% La zona C6:C25 del foglio primapar contiene l'ammontare delle vendite di
% un bene nel periodo primo-trimestre-2017  quarto-trimestre-2021.
%
% 1) Importare i dati in MATLAB e rappresentare l'andamento della serie
% storica nel periodo considerato inserendo sull'asse x del grafico le date
% nel formato QQ-yyyy (trimestre ed anno a 4 cifre)
% https://www.mathworks.com/help/matlab/ref/datetick.html#btpmlwj-1-dateFormat
%
% 2) Stimare i parametri di un modello di regressione lineare multiplo che
% spieghi l'andamento delle vendite con un trend lineare e un effetto
% dovuto alla stagionalità
%
% 3) Stimare i coefficienti del modello di regressione quando nella matrice
% X si esclude l'intercetta (prima parametrizzazione) oppure il livello
% associato alla stagione invernale (seconda parametrizzazione)
%
% 4) Confrontare ed interpretare le stime dei coefficienti di regressione
% per entrambe le parametrizzazioni
%
% 5) Costruire e commentare la serie dei valori stimati per entrambe le
% parametrizzazioni
%
% 6) Imporre il vincolo di somma a zero dei coefficienti stagionali.
% Adottando questo vincolo calcolare: la stima della componente stagionale,
% la serie destagionalizzata, la stima del trend, la serie detrendizzata e
% la serie destagionalizzata e detrendizzata.
% Rappresentare graficamente:
% a) la serie storica originaria e la stima del trend.
% b) la serie detrendizzata
% c) la serie destagionalizzata
% e) la series destagionalizzata e detrendizzata (ossia la componente
% irregolare)
%
% 7) Creare un array denominato YY con dimensione Tx4 con le seguenti
% caratteristiche
% Prima colonna: y (serie storica originaria)
% Seconda colonna: la stima del trend
% Terza colonna: stima della componente stagionale
% Quarta colonna: stima della componente irregolare
% Verificare che somma delle colonne 2:4 sia esattamente uguale alla colonna 1

% 8) Creare la timetable (denominateYYt) corrispondente all'array YY nei 4 seguenti diversi
% modi:
% Modo 1: specificare inizialmente nella sintassi di timetable la
% dimensione e le caratteristiche delle variabili
% Fare il grafico delle 4 colonne della timetable appena creata tramite la
% funzione stackedplot
% Modo 2: fornire direttamente le 4 variabili che compongono le colonne
% dell'array YY
% Modo 3: creare in anticipo la sequenza temporale delle date
% Modo 4: creare prima la table dall'array YY e poi passare alla timetable
% Modo 5: creare prima la table con solo la colonna di YY(:,1), chiamare la
% funzione addvars per aggiungere le nuove variabili alla table e poi
% convertirla in timetable
%
% 9) Prevedere l'andamento della serie storica nei periodi T+1:T+4 insieme
% ad un intervallo di confidenza al 99 per cento dei valori previsti
% Rappresentare graficamente la serie storica osservata, la serie storica
% delle previsione ed il relativo intervallo di confidenza
%
% 10) Utilizzando la funzione trenddecomp scomporre la serie storica in trend
% stagionalità e componente irregolare. 
% Rappresentare nel grafico tutte e
% 3 le componenti e confrontare il risultato con quello ottenuto in
% precedenza.
% 10a) usare l'algoritmo algoritmo STL
% 10b) usare l'algoritmo SSA (vado a stimare NumSeasonal)
% 10c) usare l'algoritmo SSA con un numero ottimo di componenti stagionali (NumSeasonal)


%% Caricamento dati in memoria
ytable=readtable('destag.xlsx','Sheet','Prima_par','Range','B5:C25','ReadRowNames',true);
y=ytable{:,1};
T=size(ytable,1);
% mo =vettore colonna che contiene la sequenza 1,2,3,4 ripetuta 5 volte
mo = repmat((1:4)',5,1);
% Xseaso = matrice che contiene le variabili stagionali dummies
Xseaso = dummyvar(mo);
% Xtrend vettore colonna che contiene la sequenza 1, 2, ...., T
Xtrend =(1:T)';
% one = vettore colonna con tutti gli elementi uguali a 1 di lunghezza T
one=ones(T,1);

% Costruire la sequenza di date delle serie storica in formato datetime
% Osservazione: il testo dell'esercizio afferma che la data di inizio è 
% il primo-trimestre-2017. Di conseguenza qualsiasi data all'interno del
% primo trimestre 2017 (ad esempio primo gennaio) va bene.

% Osservazione le variabili datainizio, cadenzaRilevazione e datafine
% verranno utilizzate diverse volte all'interno del file.
datainizio=datetime('2017-01-1');
cadenzaRilevazione=calquarters(1);
datafine=datetime('2021-10-1');

% rowTimes è un vettore riga di lunghezza T che contiene le date
rowTimes=datainizio:cadenzaRilevazione:datafine;


%% Rappresentazione grafica della serie storica
plot(rowTimes,y,'b-')
xlabel('t')
ylabel('Valori effettivi')
datetick('x','QQ-yyyy')


%% Costruzione matrice X prima parametrizzazione
figure
Xprimapar=[Xtrend Xseaso];
bprimapar=Xprimapar\y;
% yhat primapar è il vettore che contiene i valori teorici (trend +
% stagionalità) in base alla prima parametrizzazione
yhatprimapar=Xprimapar*bprimapar;
% Grafico valori effettivi e valori teorici
plot(rowTimes,y,'b-',rowTimes,yhatprimapar,'r--')
xlabel('t')
ylabel('Valori effettivi e valori teorici (trend+ stagionalità)')
datetick('x','QQ-yyyy')

%% Costruzione matrice X seconda parametrizzazione
Xsecondapar=[one Xtrend Xseaso(:,1:end-1)];
bsecondapar=Xsecondapar\y;
yhatsecondapar=Xsecondapar*bsecondapar;

% I valori stimati utilizzando la prima e la seconda parametrizzazione non
% cambiano
disp('Confronto tra i valori stimati ottenuti con le due parametrizzazioni')
disp([yhatprimapar yhatsecondapar])

%% Imposizione vincolo somma a zero dei coefficienti stagionali
% 6) Imporre il vincolo di somma a zero dei coefficienti stagionali.
% Adottando questo vincolo calcolare: la stima della componente stagionale,
% la serie destagionalizzata, la stima del trend, la serie detrendizzata e

btrend=bprimapar(1);
bseaso=bprimapar(2:end);
% La media dei coefficienti stagionali è l'intercetta
intercept=mean(bseaso);
% bseasosommazero = vettore degli scostamenti dalla media
% La media di questo vettore è zero. bseasosommazero è interpretabile come
% gli effetti stagionali rispetto alla media.
bseasosommazero=bseaso-intercept;
bterzapar=[intercept;btrend;bseasosommazero];
rownam=["Intercetta"; "trend"; "Q"+(1:4)']; 
bterzapartable=array2table(bterzapar,'RowNames',rownam);
disp('Coefficienti stagionali a somma zero')
disp('Si noti che in questo caso la somma dei coefficienti stagionali è 0')
disp(bterzapartable)

%% plot y and trend
% a) la serie storica originaria e la stima del trend.

figure
plot(rowTimes,y,'b-')
hold('on')
yhattrend=intercept+btrend*Xtrend;
plot(rowTimes,yhattrend,'k--')
datetick('x','QQ-yyyy')
title('Serie storica originale e stima del trend lineare')

%% Stima serie detrendizzata
% b) la serie detrendizzata

figure
plot(rowTimes,y-yhattrend,'k--')
title('Serie detrendizzata=stagionalità+irregolare')
datetick('x','QQ-yyyy')

%% Stima della componente stagionale
figure
yhatseaso=Xprimapar(:,2:end)*bseasosommazero;
plot(rowTimes,yhatseaso,'k--')
title('Stima della componente stagionale')
datetick('x','QQ-yyyy')


%% Stima della serie  destagionalizzata
% c) la serie destagionalizzata

figure
yhatdestag=y-yhatseaso;
plot(rowTimes,yhatdestag,'k--')
title('Stima della serie destagionalizzata=trend+irregolare')
datetick('x','QQ-yyyy')

%% Stima della serie  detrendizzata e destagionalizzata
% e) la series destagionalizzata e detrendizzata (ossia la componente
% irregolare)

figure
yhatirregular=y-yhatseaso-yhattrend;
plot(rowTimes,yhatirregular,'k--')
title('Stima della serie destagionalizzata e detrendizzata=irregolare')
datetick('x','QQ-yyyy')



%% Creazione dell'array YY
% YY ha dimensione Tx4 e contiene rispettivamente
% Prima colonna: y (serie storica originaria)
% Seconda colonna: la stima del trend
% Terza colonna: stima della componente stagionale
% Quarta colonna: stima della componente irregolare
% La somma delle colonne 2:4 è esattamente uguale alla colonna 1
YY=[y yhattrend yhatseaso yhatirregular];
% Osservazione: YY è un array di conseguenza non contiene nomi delle righe
% e/o delle colonne

%% Verificare che somma delle colonne 2:4 sia esattamente uguale alla colonna 1
controllougu=isequal(YY(:,1),sum(YY(:,2:4),2));
% Osservazione: il messaggio di errore che segue appare solo se la
% variabile controllougu è un booleano di tipo false
assert(controllougu,'Errore: la scomposizione della serie è errata')

%% 8) Creare la timetable corrispondendente all'array YY nei 4 seguenti diversi
% modi:
% Modo 1: specificare inizialmente nella sintassi di timetable la
% dimensione e le caratteristiche delle variabili
% Modo 2: forniscono direttamente le 4 variabili che compongono le colonne
% dell'array YY
% Modo 3: crea in anticipo la sequenza temporale delle date
% Modo 4: creare prima la table e poi passare alla timetable
%

%% Modo 1 per creare la timetable (si specifica la dimensione e le caratteristiche delle variabili)
% Dimensione  delle timetable Tx4
dimYY=[T 4];
% Nomi della timetable 
nomiYY={'y','trend','componente stagionale', 'componente irregolare'};
% Tipo di variabili nella timetable. Le nostre 4 variabili sono tutte di tipo double
vartypeYY=repelem("double",4,1);
% datainizio della timetable
YYt=timetable('Size',dimYY,'VariableTypes',vartypeYY,'VariableNames',nomiYY,...
    'StartTime',datainizio,'TimeStep',cadenzaRilevazione);
disp('Table appena creata con date e nomi e delle variabili ma con valori pari a 0')
disp(YYt)

% Nell'istruzione di seguito vado a popolare YYt con YY
disp('Table popolata con i valori delle 4 variabili')
YYt{:,:}=YY;
disp(YYt)
stackedplot(YYt)
title('Scomposizione della serie storica in termini di trend, stagionalità e componente irregolare')

%% Modo 2 per creare la timetable (si forniscono direttamente le variabili)

nomiYY={'y','trend','componente stagionale', 'componente irregolare'};

YYt=timetable(y,yhattrend,yhatseaso,yhatirregular,'VariableNames',nomiYY,...
    'StartTime',datainizio,'TimeStep',cadenzaRilevazione);
disp(YYt)
% Nuova figura per ospitare lo stackedplot di YYt
figure
stackedplot(YYt)
title('Scomposizione della serie storica in termini di trend lineare, stagionalità e componente irregolare')

%% Modo 3 per creare la timetable (si crea in anticipo la sequenza temporale delle date)
% In questo caso si crea in anticipo quello che dovrà essere inserito nella prima
% colonna della timetable

YYtchk=timetable(y,yhattrend,yhatseaso,yhatirregular,'VariableNames',nomiYY,...
   'RowTimes',rowTimes);

%% Modo 4 per creare la timetable (si passa attraverso la table)
% Creare prima la table dall'array YY e poi passare alla timetable

% Si crea prima la table
YYtable=array2table(YY,'VariableNames',nomiYY);
% si converte la time in timetable
YYtchk1=table2timetable(YYtable,'RowTimes',rowTimes);

%% Modo 5 per creare la timetable (si passa attraverso la table e si chiama addvars)
% creare prima la table con solo la colonna di YY(:,1), chiamare la
% funzione addvars per aggiungere le nuove variabili alla table e poi
% convertirla in timetable
ytable=table(YY(:,1));
YYtable=addvars(ytable,yhattrend,yhatseaso,yhatirregular,'NewVariableNames',nomiYY(2:end));
% si converte la time in timetable
YYt=table2timetable(YYtable,'RowTimes',rowTimes);

%% 9) Prevedere l'andamento della serie storica nei periodi T+1:T+4 
% insieme ad un intervallo di confidenza al 99 per cento dei valori
% previsti

% Stimo il modello con la funzione fitlm
% Dentro fitlm l'intercetta viene aggiunta di default perciò devo escludere
% la prima colonna che contiene il vettore di uno
out=fitlm(Xsecondapar(:,2:end),y);

nuovit=(T+1:T+4)';
% Coordinate della matrice X se si adotta la seconda parametrizzazione in
% corrispondenza delle 4 stagioni del 2022
Xnew=[nuovit,eye(4,3)]; 
[ypred,yci]=predict(out,Xnew,'Alpha',0.01,'Prediction','observation');

%% Grafico valori effettivi + valori previsti per il 2022 ed intervallo di confidenza

% Creazione della sequenza in formato datetime per il 2022
% Costruire la sequenza di date delle serie storica in formato datetime
datainizioP=datetime('2022-01-1');
datafineP=datetime('2022-10-1');
Xtempi4futuri=datainizioP:cadenzaRilevazione:datafineP;
% Grafico serie originaria e trend nel periodo di rilevazione
plot(rowTimes,YY(:,1),'b-',rowTimes,YY(:,2),'k--')
hold('on')
plot(Xtempi4futuri,[ypred yci],'r-.','LineWidth',2)
datetick('x','QQ-yyyy')
title('Serie originale, stima del trend ed intervallo di confidenza al 99% per il 2022')


%% 10a) Scomposizione della serie storica tramite trenddecomp usando algoritmo STL
[yhattrend1,yhatseaso1,yhatirregular1] = trenddecomp(y,"stl",4);
YYt1=timetable(y,yhattrend1,yhatseaso1,yhatirregular1,'VariableNames',nomiYY,...
    'StartTime',datainizio,'TimeStep',cadenzaRilevazione);
disp(YYt1)


% Nuova figura per ospitare lo stackedplot di YYt1
figure
stackedplot(YYt1)
title('Scomposizione della serie storica tramite algoritmo STL')

% YYt conteneva la vecchia scomposizione parametrica
figure
stackedplot(YYt)
title('Scomposizione della serie storica in termini di trend lineare, stagionalità e componente irregolare')

%% 10b) Scomposizione della serie storica tramite algoritmo SSA (vado a stimare NumSeasonal)
[yhattrend1,yhatseaso1,yhatirregular1] = trenddecomp(y,NumSeasonal=4);
% Dato che non so quanti cicli stagionali (semestrali, quadrimestrali, trimestrali, mensili...)
% ci possono essere inserisco NumSeasonal=4
YYt1=timetable(y,yhattrend1,yhatseaso1,yhatirregular1,'VariableNames',nomiYY,...
    'StartTime',datainizio,'TimeStep',cadenzaRilevazione);
stackedplot(YYt1)
title('Scomposizione della serie storica tramite algoritmo SSA (4 sottocomponenti stagionali)')

% Dall'esame delle componente stagionali nello stackedplot mi accorgo che
% le componenti stagionali 3 e 4 sono molto vicine a zero. 
% La linea blu (componente stagionale 1) sembra essere collegata ad una
% componente semestrale La linea rossa riflette il fatto che le vendite nei
% trimestri 3 e 4 sembrano essere più alte di quelle dei trimestri 1 e 2.
% Nella section di seguito inserisco NumSeasonal=2 e calcolo la componente
% stagionale come somma delle due sottocomponenti stagionali che sono state
% trovate.


%% 10c) Scomposizione della serie storica tramite algoritmo SSA (fisso NumSeasonal=2)
[yhattrend1,yhatseaso1,yhatirregular1] = trenddecomp(y,NumSeasonal=2);

YYt1=timetable(y,yhattrend1,sum(yhatseaso1,2),yhatirregular1,'VariableNames',nomiYY,...
    'StartTime',datainizio,'TimeStep',cadenzaRilevazione);
figure
stackedplot(YYt1)
title('Scomposizione della serie storica tramite algoritmo SSA')