clear all;
close all;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GARCH(p,q)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parametri DGP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10000;
alpha=0.05;

% Scegli il tipo di DGP
% Opzioni: 'gaussian', 'arch', 
% 'garch1', 'garch2', 'garcht'
% 'egarch'
dgp = 'egarch';
% Parametri del modello GARCH
K = 0.05;
coef_arch=[0.31 0.26];
coef_garch=[0.22 0.15];
lev = [0.15 0.5];
% Se disturbi di tipo t allora 
% specifica i gradi di liberta'
dof=5;
% Primi elementi della serie
% da scartare
trans = 100;
% Coeff. modello lineare
beta_teo = [3 -1 2.5 0 -1.5]';
p = length(beta_teo)-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parametri dell'algoritmo di stima
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Qest = 2;  % Ordine ARCH
Pest = 2;  % Ordine GARCH
varmod_est = 'egarch';
distrib_est = 'gaussian';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Data Generating Process (DGP)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xobs = rand(N,p);
X = [ones(N,1) Xobs];
switch dgp
    case 'gaussian'
        Q = 0; P = 0; L = 0;
        distrib = 'Gaussian';                   
        % I.i.d. gaussian disturbances
        u = randn(N,1);
        y = X*beta_teo+sigma*u;
    case 'arch'
        % Genera i dati usando la definizione di ARCH
        Q = numel(coef_arch); P = 0; L = 0;
        modello = ['ARCH(',num2str(Q),')'];            
        u = randn(N+trans,1);
        h = u;
        for t=(Q+1):(N+trans)
            h(t) = K;
            for j=1:Q
                h(t)=h(t)+coef_arch(j)*(u(t-j))^2;
            end
            u(t)=randn(1)*sqrt(h(t));
        end
        u=u(trans+1:end);
        y = X*beta_teo+u;
    case 'garch1'
        % Genera i dati usando la definizione di GARCH
        P = numel(coef_garch); Q = numel(coef_arch); L = 0;        
        modello = ['GARCH(',num2str(P),',',num2str(Q),')'];                
        ord = max(Q,P);
        u = randn(N+trans,1);
        h = u;
        for t=(ord+1):(N+trans)
            h(t) = K;
            for j=1:Q
                h(t)=h(t)+coef_arch(j)*(u(t-j))^2;
            end
            for j=1:P
                h(t)=h(t)+coef_garch(j)*h(t-j);
            end
            u(t)=randn(1)*sqrt(h(t));
        end
        u=u(trans+1:end);
        y = X*beta_teo+u;
    case 'garch2'
        % Genera i dati con il comando garchsim
        P = numel(coef_garch); Q = numel(coef_arch); L = 0;        
        modello = ['GARCH(',num2str(P),',',num2str(Q),')'];                
        distrib = 'Gaussian';
        varmod = 'Garch';        
        spec = garchset('C',0,'R',0,'AR',[],'M',0,'MA',[],...
            'VarianceModel',varmod,'P',P,'Q',Q,'K',K, ...
            'GARCH',coef_garch,'ARCH',coef_arch,'Distribution',distrib);
        [u, s, y] = garchsim(spec, N, 1);
    case 'garcht'
        % Genera i dati con il comando garchsim
        P = numel(coef_garch); Q = numel(coef_arch); L = 0;        
        modello = ['GARCH-T(',num2str(P),',',num2str(Q),')'];                
        distrib = 'T';
        varmod = 'Garch';        
        spec = garchset('C',0,'R',0,'AR',[],'M',0,'MA',[],...
            'VarianceModel',varmod,'P',P,'Q',Q,'K',K, ...
            'GARCH',coef_garch,'ARCH',coef_arch,...
            'Distribution',distrib,'DoF',dof);
        [u, s, y] = garchsim(spec, N, 1);        
    case 'egarch'
        % Genera i dati con il comando garchsim
        P = numel(coef_garch); Q = numel(coef_arch); L = numel(lev);            
        modello = ['EGARCH(',num2str(P),',',num2str(Q),')'];                
        distrib = 'Gaussian';
        varmod = 'EGarch';           
        spec = garchset('C',0,'R',0,'AR',[],'M',0,'MA',[],...
            'VarianceModel',varmod,'P',P,'Q',Q,'K',K, ...
            'GARCH',coef_garch,'ARCH',coef_arch,...
            'Distribution',distrib,'Leverage',lev);
        [u, s, y] = garchsim(spec, N, 1);         
end

%%%%%%%%%%%%%%
%% Specifica la struttura del modello da stimare
%%%%%%%%%%%%%%
garch_spec = garchset('Distribution',distrib_est,'R',0,'M',0,'C',0,...
    'VarianceModel',varmod_est,'Display','On','P',Pest,'Q',Qest,'K',1);
% Stima i parametri del modello
[parms,errors,LLF,resid,stderrors,summary] = garchfit(garch_spec,u);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Confrontiamo valori veri con valori stimati
%% Significativita' statistica dei parametri
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fase preparatoria di alcuni oggetti 
% che servono per costruire la tabella
if Pest==0 
    parms.GARCH = [];
    errors.GARCH = [];
    parms.Leverage = [];
    errors.Leverage = [];  
elseif Pest~=0 && ~strcmp(dgp,'egarch') 
    parms.Leverage = [];
    errors.Leverage = [];
end

valori_veri=K;
if Q>0
    valori_veri = [valori_veri coef_arch];
end
if Qest>Q
    valori_veri = [valori_veri zeros(1,Qest-Q)];
end
if Pest>0 && ~strcmp(dgp,'arch')
    valori_veri = [valori_veri coef_garch];
end
if Pest>P
    valori_veri = [valori_veri zeros(1,Pest-P)];
end
if L>0
    valori_veri = [valori_veri lev];
end
if numel(parms.Leverage)>L
    valori_veri = [valori_veri zeros(1,Lest-L)];
end

%%%%%%%%%%%%%%
%% Alcuni calcoli
%%%%%%%%%%%%%%
k = numel(parms.ARCH)+numel(parms.GARCH)+numel(parms.Leverage)+1;
pars = [parms.K parms.ARCH parms.GARCH parms.Leverage]';
se_pars = [errors.K errors.ARCH errors.GARCH errors.Leverage]';
pars_tstat = pars./se_pars;
tstat_pvals = 2*(1-tcdf(abs(pars_tstat),N-k));
sigma = sqrt((resid'*resid)/N);

% Fase preparatoria di altri oggetti 
% che servono per costruire la tabella
mtx = [valori_veri' pars se_pars pars_tstat tstat_pvals ...
       pars+tinv(alpha/2,N-k)*se_pars pars+tinv(1-alpha/2,N-k)*se_pars];
info.cnames = strvcat('Veri','Stimati','S.E.','t-stat','p-value',...
    '[95% Conf.','Interval]');
info.rnames = strvcat(' ','K');
for j=1:Qest
    info.rnames = strvcat(info.rnames,['ARCH',num2str(j)]);
end
for j=1:Pest
    info.rnames = strvcat(info.rnames,['GARCH',num2str(j)]);
end
for j=1:numel(parms.Leverage)
    info.rnames = strvcat(info.rnames,['Lev',num2str(j)]);
end
info.fmt=strvcat('%8.2f','%8.5f','%8.5f','%8.2f','%8.4f','%8.4f','%8.4f');
info.fid = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Tabella dei risultati
%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'\n\n***************************************************************\n\n');
fprintf(1,['Stima parametri del modello ',modello,'\n']);
fprintf(1,'No. of obs     = %9d \n',N);
fprintf(1,'Log-Likelihood = %9.2f \n',LLF);
fprintf(1,'Sigma          = %9.2f \n\n',sigma);
mprint(mtx,info)
fprintf(1,'***************************************************************\n');

%%%%%%%%%%
%% Plot
%%%%%%%%%%
figure(1);
subplot(2,1,1);
plot(u)
title(modello)
subplot(2,2,3);
hold on;
[dens,xi]=ksdensity(u);
normdens=normpdf(xi,mean(u),std(u));
plot(xi,dens,'-b',xi,normdens,'-r');
axis([-7.5 7.5 get(gca,'ylim')]);
box on;
title('Histogram')
hold off;
subplot(2,2,4);
qqplot(u);