% Sample size
n=4000;
n
% ng must be greater than the sample size
ng=18000;

% parameters

alpha_0=1;
beta_1=0.3;
beta_2=0.1;
alpha_1=0.2;
alpha_2=0.1;
v=5;
mu=2;

a=ones(ng,1)*mu;
h=(alpha_0/(1-alpha_1-alpha_2-beta_1-beta_2))*ones(ng,1);

for i=3:ng
    h(i) = alpha_0+alpha_1*(a(i-1)-mu).^2+alpha_2*(a(i-2)-mu).^2+beta_1*h(i-1)+beta_2*h(i-2);
    a(i) = mu+sqrt(h(i))*trnd(v);
end

t=(1:1:ng-n);
a=a(n+1:ng);

% Estimation of GARCH model 1

dist = struct('Name','Gaussian');

ModelToEst = garch('Offset',NaN,'Constant',NaN,'ARCHLags',[1,2],'GARCHLags',[1 2],'Distribution',dist);
[EstMdl, EstParamCov, LogL, info] = estimate(ModelToEst,a);

EstMdl
EstParamCov

% Estimation of GARCH model 2

dist = struct('Name','t','DoF',NaN);

ModelToEst = garch('Offset',NaN,'Constant',NaN,'ARCHLags',[1,2],'GARCHLags',[1 2],'Distribution',dist);
[EstMdl, EstParamCov, LogL, info] = estimate(ModelToEst,a);

EstMdl
EstParamCov


figure(1)
plot(t,a)
box off
grid off
xlabel('Time Index')
ylabel('Value of GARCH(2,1) process')
title('GARCH(2,1)')


figure(2)
normplot(a)
title('GARCH(2,1)')



