clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MA(q) process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 100;
theta = [-0.6 0.0 -0.5 0.5];
mu = 1.3;
nlag = 25;
q = numel(theta);
nt = n+q;
sigma= 1.5;

%%%%%%%%
%% DGP
%%%%%%%%
% Gaussian white noise
eps = sigma*randn(nt,1);
% MA(4) process
Y = mu + eps(q+1:end)+theta(1)*eps(q:end-1)+...
    theta(2)*eps(q-1:end-2)+theta(3)*eps(q-2:end-3)+...
    theta(4)*eps(q-3:end-4);
% Metodo generale per q qualsiasi
Y = mu + eps(q+1:end);
for j=1:q
    Y = Y+ theta(j)*eps((q-j+1):(end-j));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Traccia la serie nella sua
%% evoluzione temporale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
plot(Y)
axis tight
title(['Realizzazione di un processo MA(',num2str(q),')'])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Il comando acf stima 
%% i coeff. di autocorrelazione.
%% Li si visualizza con bar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
acf1 = acf(Y,nlag);
hold on
bar(1:acf1.k,acf1.ac)
line(1:acf1.k,acf1.lowb,'LineStyle',':','Color','r')
line(1:acf1.k,acf1.topb,'LineStyle',':','Color','r')
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Il comando sacf stima e visualizza
%% i coeff. di autocorrelazione
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
axis tight
acf1b = sacf(Y,nlag,0);
% Il terzo parametro di sacf e' gflag:
% gflag = 0, flag for graphing, 1 = no graph

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Il comando spacf stima e visualizza
%% i coeff. di autocorrelazione parziale
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
axis tight
pacf1 = spacf(Y,nlag);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima parametri modello ARMA
%% mediante comando armaxfilter dell'Econometrics 
%% Toolbox di LeSage (non funziona su Matlab 7.4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[parms, errors, LLF , SEregression, stderrors, robustSE, scores, likelihoods]=armaxfilter(Y,1,0,4);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima parametri modello ARMA
%% mediante comando garchfit del GARCH
%% Toolbox di Mathworks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definisci specifiche del modello
arma_spec = garchset('Distribution','Gaussian','R',0,'M',4,'C',1,...
    'VarianceModel','Constant','Display','On');
% Stima i paramentri del modello
[parms,errors,LLF,resid,stderrors,summary] = garchfit(arma_spec, Y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Confrontiamo valori veri con valori stimati
%% Significativita' statistica dei parametri
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=0.05;
k = numel(theta)+1;
fprintf(1,'\n%s\n','************************************************');
fprintf(1,'%s\n','Valori veri');
fprintf(1,'%4.5f ',[mu theta]);
fprintf(1,'\n%s\n','Valori stimati');
pars = [parms.C parms.MA];
fprintf(1,'%4.5f ',pars);
fprintf(1,'\n%s\n','Std.Errors');
se_pars = [errors.C errors.MA];
fprintf(1,'%4.5f ',se_pars);
fprintf(1,'\n%s\n','t-statistics');
pars_tstat = pars./se_pars;
fprintf(1,'%4.5f ',pars_tstat);
fprintf(1,'\n%s\n','p-values');
tstat_pvals = 2*(1-tcdf(abs(pars_tstat),n-k));
fprintf(1,'%4.5f ',tstat_pvals);
fprintf(1,'\n%s\n','sigma');
sigma = sqrt((resid'*resid)/n);
fprintf(1,'%4.5f',sigma);
fprintf(1,'\n%s\n','************************************************');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Verifichiamo la struttura
%% di autocorrelazione dei
%% residui del modello
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
acf2 = sacf(resid,nlag);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% AR(p) process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 300;
phi = [0.8 -0.9 0.4];
mu = -0.9;
sigma=1.5;
nlag = 25;
p = length(phi);
nt = n+p;

% Radici polinomio caratteristico
figure
ro = roots([fliplr(-phi) 1]);
hold on;
plot(ro,'x','MarkerSize',12,'Color','r')
plot(ro,'.','MarkerSize',18,'Color','r')
c1=rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],...
    'EdgeColor',[.3 .3 1],'LineWidth',1.5)
axis equal
xlim = get(gca,'xlim'); ylim = get(gca,'ylim');
set(gca,'xlim',[xlim(1)-.1 xlim(2)+.1]);
set(gca,'ylim',[ylim(1)-.1 ylim(2)+.1]);
box on
title(['Processo AR(',num2str(p),...
    ') - Radici polinomio caratteristico']);
hold off;

%%%%%%%%
%% DGP
%%%%%%%%
% Gaussian white noise
eps = sigma*randn(nt,1);
% AR(q)
Y = mu+zeros(nt,1);
Y(1:p) = randn(p,1);
for t=p+1:n
    for j=1:p
        Y(t) = Y(t)+ phi(j)*Y(t-j);
    end
    Y(t)=Y(t)+eps(t);
end

% Traccia la serie
figure
plot(Y)
axis tight
title(['Realizzazione di un processo AR(',num2str(p),')'])

% ACF e PACF
figure
acf1 = acf(Y,nlag);
hold on
bar(1:acf1.k,acf1.ac)
line(1:acf1.k,acf1.lowb,'LineStyle',':','Color','r')
line(1:acf1.k,acf1.topb,'LineStyle',':','Color','r')
hold off
figure
axis tight
acf1b = sacf(Y,nlag);
figure
axis tight
pacf1 = spacf(Y,nlag);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima parametri modello ARMA
%% mediante comando armaxfilter dell'Econometrics 
%% Toolbox di LeSage (non funziona su Matlab 7.4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[parms, errors, LLF , SEregression, stderrors, robustSE, scores, likelihoods]=armaxfilter(Y,1,p,0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima parametri modello ARMA
%% mediante comando garchfit del GARCH
%% Toolbox di Mathworks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 3;
q = 3;
% Definisci specifiche del modello
arma_spec = garchset('Distribution','Gaussian','R',p,'M',q,'C',1,...
    'VarianceModel','Constant','Display','On');
% Stima i paramentri del modello
[parms,errors,LLF,resid,stderrors,summary] = garchfit(arma_spec, Y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Confrontiamo valori veri con valori stimati
%% Significativita' statistica dei parametri
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=0.05;
k = numel(theta)+1;
fprintf(1,'\n%s\n','************************************************');
fprintf(1,'%s\n','Valori veri');
fprintf(1,'%4.5f ',[mu phi zeros(1,q)]);
fprintf(1,'\n%s\n','Valori stimati');
pars = [parms.C parms.AR parms.MA];
fprintf(1,'%4.5f ',pars);
fprintf(1,'\n%s\n','Std.Errors');
se_pars = [errors.C errors.AR errors.MA];
fprintf(1,'%4.5f ',se_pars);
fprintf(1,'\n%s\n','t-statistics');
pars_tstat = pars./se_pars;
fprintf(1,'%4.5f ',pars_tstat);
fprintf(1,'\n%s\n','p-values');
tstat_pvals = 2*(1-tcdf(abs(pars_tstat),n-k));
fprintf(1,'%4.5f ',tstat_pvals);
fprintf(1,'\n%s\n','sigma');
sigma = sqrt((resid'*resid)/n);
fprintf(1,'%4.5f',sigma);
fprintf(1,'\n%s\n','************************************************');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Verifichiamo la struttura
%% di autocorrelazione dei
%% residui del modello
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
acf2 = sacf(resid,nlag);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ARMA(p,q) process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 500;
phi = [0.8 -0.9 0.4];
theta = [-0.6 0.3 -0.5 0.5];
mu = -0.9;
sigma = 1.2;
nlag = 25;
p = length(phi);
q = length(theta);
pq=max(p,q);
% Transiente: quanti dati iniziali 
% generare e poi buttare
trans = 100; 
% Quanti dati generare per la
% successiva fase di previsione
nprev = 100;
% Lunghezza totale della serie
nt = n+trans+nprev;

% Radici del polinomio caratteristico
figure
ro = roots([fliplr(-phi) 1]);
hold on;
plot(ro,'x','MarkerSize',12,'Color','r')
plot(ro,'.','MarkerSize',18,'Color','r')
c1=rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],...
    'EdgeColor',[.3 .3 1],'LineWidth',1.5)
axis equal
xlim = get(gca,'xlim'); ylim = get(gca,'ylim');
set(gca,'xlim',[xlim(1)-.1 xlim(2)+.1]);
set(gca,'ylim',[ylim(1)-.1 ylim(2)+.1]);
box on
title(['Processo ARMA(',num2str(p),',',num2str(q),...
    ') - Radici polinomio caratteristico']);
hold off;

%%%%%%%%
%% DGP
%%%%%%%%
% Gaussian white noise
eps = sigma*randn(nt,1);
% AR(q)
Y = mu+zeros(nt,1);
Y(1:pq) = randn(pq,1);
for t=(pq+1):nt
    for j=1:p
        Y(t) = Y(t) + phi(j)*Y(t-j);
    end
    for j=1:q
        Y(t) = Y(t) + theta(j)*eps(t-j);
    end
    Y(t)=Y(t)+eps(t);
end
Youtsample = Y(end-nprev+1:end);
Y = Y(trans+1:trans+n);

% Traccia la serie
figure
plot(Y)
axis tight
title(['Realizzazione di un processo ARMA(',...
    num2str(p),',',num2str(q),')'])

% ACF e PACF
figure
acf1 = acf(Y,nlag);
hold on
bar(1:acf1.k,acf1.ac)
line(1:acf1.k,acf1.lowb,'LineStyle',':','Color','r')
line(1:acf1.k,acf1.topb,'LineStyle',':','Color','r')
hold off
figure
axis tight
acf1b = sacf(Y,nlag);
figure
axis tight
pacf1 = spacf(Y,nlag);
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima parametri modello ARMA
%% mediante comando armaxfilter dell'Econometrics 
%% Toolbox di LeSage (non funziona su Matlab 7.4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[parms, errors, LLF , SEregression, stderrors, robustSE, scores, likelihoods]=armaxfilter(Y,1,p,q);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima parametri modello ARMA
%% mediante comando garchfit del GARCH
%% Toolbox di Mathworks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 4;
q = 4;
% Definisci specifiche del modello
arma_spec = garchset('Distribution','Gaussian','R',p,'M',q,'C',1,...
    'VarianceModel','Constant','Display','On');
% Stima i paramentri del modello
[parms,errors,LLF1,resid,stderrors,summary] = garchfit(arma_spec, Y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Confrontiamo valori veri con valori stimati
%% Significativita' statistica dei parametri
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=0.05;
k = numel(theta)+1;
fprintf(1,'\n%s\n','************************************************');
fprintf(1,'%s\n','Valori veri');
fprintf(1,'%4.5f ',[mu phi zeros(1,p-numel(phi)) theta zeros(1,q-numel(theta))]);
fprintf(1,'\n%s\n','Valori stimati');
pars = [parms.C parms.AR parms.MA];
fprintf(1,'%4.5f ',pars);
fprintf(1,'\n%s\n','Std.Errors');
se_pars = [errors.C errors.AR errors.MA];
fprintf(1,'%4.5f ',se_pars);
fprintf(1,'\n%s\n','t-statistics');
pars_tstat = pars./se_pars;
fprintf(1,'%4.5f ',pars_tstat);
fprintf(1,'\n%s\n','p-values');
tstat_pvals = 2*(1-tcdf(abs(pars_tstat),n-k));
fprintf(1,'%4.5f ',tstat_pvals);
fprintf(1,'\n%s\n','sigma');
sigma = sqrt((resid'*resid)/n);
fprintf(1,'%4.5f',sigma);
fprintf(1,'\n%s\n','************************************************');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Verifichiamo la struttura
%% di autocorrelazione dei
%% residui del modello
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
acf2 = sacf(resid,nlag);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LR-test ARMA(4,4) vs. ARMA(3,4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 3;
q = 4;
% Definisci specifiche del modello
arma_spec = garchset('Distribution','Gaussian','R',p,'M',q,'C',1,...
    'VarianceModel','Constant','Display','Off');
% Stima i paramentri del modello
[parms2,errors2,LLF2,resid2,stderrors2,summary2] = garchfit(arma_spec, Y);

LRstat = 2*(LLF1-LLF2);
chi2tab = chi2inv(0.95,1);
LRprob = 1-chi2cdf(LRstat,1);

fprintf(1,'\n%s\n','Likelihood Ratio Test');
fprintf(1,'%s\n','***************************************************************');
fprintf(1,'LR = = %4.2f    ',LRstat);
fprintf(1,'Chi-sq. = %4.2f    ',chi2tab);
fprintf(1,'Prob. = %4.6f\n',LRprob);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LR-test ARMA(3,4) vs. ARMA(3,3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 3;
q = 3;
% Definisci specifiche del modello
arma_spec = garchset('Distribution','Gaussian','R',p,'M',q,'C',1,...
    'VarianceModel','Constant','Display','Off');
% Stima i paramentri del modello
[parms3,errors3,LLF3,resid3,stderrors3,summary3] = garchfit(arma_spec, Y);

LRstat = 2*(LLF2-LLF3);
chi2tab = chi2inv(0.95,1);
LRprob = 1-chi2cdf(LRstat,1);

fprintf(1,'\n%s\n','Likelihood Ratio Test');
fprintf(1,'%s\n','***************************************************************');
fprintf(1,'LR = = %4.2f    ',LRstat);
fprintf(1,'Chi-sq. = %4.2f    ',chi2tab);
fprintf(1,'Prob. = %4.6f\n',LRprob);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Previsione usando modello ARMA stimato
%% Comando garchpred del GARCH
%% Toolbox di Mathworks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Previsione ad 1 step avanti nel tempo
oriztemp = 1;
% Definisci le stesse specifiche del modello ARMA(p,q)
% usate per garchfit
p = numel(phi);
q = numel(theta);
arma_spec = garchset('Distribution','Gaussian','R',p,'M',q,...
    'C',parms2.C,'AR',parms2.AR,'MA',parms2.MA,'K',NaN,...
    'VarianceModel','Constant','Display','On');
Yprev = zeros(nprev,1);
pq = max(p,q);
for t=pq:nprev
    [SigmaForecast, Yprev(t)] = ...
        garchpred(arma_spec, Youtsample(t-pq+1:t), oriztemp);
end
Yve = Youtsample(pq+1:nprev);
Ypr = Yprev(pq:nprev-1);
eprev = (Yve-Ypr);
Theil = sqrt(eprev'*eprev)/sqrt(sum((Yve-mean(Y)).^2));

% Traccia valore vero vs. valore previsto
figure
plot(Yve,Ypr,'.')
mn = min([Yve; Ypr]);
mx = max([Yve; Ypr]);
line([mn mx],[mn mx],'Color','red','LineWidth',2,'LineStyle',':')
set(gca,'xlim',[mn mx],'ylim',[mn mx],'dataaspectratio',[1 1 1])
title(['Previsione di un ARMA(',num2str(p),',',num2str(q),...
    ') - Theil index = ',num2str(round(Theil*100)/100)])
xlabel('Valore vero')
ylabel('Valore previsto')