/* This procedure simulates a GARCH(1,1) procedure
   Written by Christos S. Savva */

new;cls;
library maxlik,pgraph;
#include maxlik.ext;
maxset;
__title="GARCH Simulation Process";

_max_ParNames = "c" | "w" | "a" | "b" ;


_max_GradTol=1.e-5;                
_max_Algoritm = 2;
_max_CovPar = 1;
_max_MaxIters = 150;
__output=1;     


/* Monte Carlo simulation*/

rep=5;
j=1;
n=1100;
c0=1;
w0=1;
a0=0.15;
b0=0.85;
let estimates={};

do while j<=rep; 
   {yt}=datagen(n,c0,w0,a0,b0);
   T=rows(yt); /* sample size*/


/*combined vector of starting values of the parameter*/
  
   bstart=c0|w0|a0|b0;

  {bhat,mllf1,grad1,cov1,rc1}=maxprt(maxlik(yt,0,&garchll,bstart));

estimates=estimates|bhat';

j=j+1;
endo;

print "no of successful convergence";
rows(estimates);

print "mean estimates";
meanc(estimates);



/*********************************************************************************/

/* Data generation Proc*/

proc(1)= datagen(n,c0,w0,a0,b0); 
  local ini,ksi,y,h,e_t,i, hup, hlow;  

ini=100;

ksi=rndn(n,1);
y=zeros(n,1);
h=zeros(n,1);
e_t=zeros(n,1);
h[1]=1;

y[1]=0;

       for i(2,n,1);
          h[i]= w0 + (a0*ksi[i-1]^2 + b0)*h[i-1];
          hup=1e6; hlow=1e-6;
          if h[i]>hup;
             h[i]=hup;
          endif;
          if h[i]<hlow;
             h[i]=hlow;
          endif;
          e_t[i]=ksi[i]*sqrt(h[i]);
          y[i] = c0 + e_t[i];
       endfor;


y=y[ini+1:n];
h=h[ini+1:n];
e_t=e_t[ini+1:n];


retp(y);
endp;



/****************************************************************************/



/*Likelihood proc*/

proc garchll(b,z);

  local res,res_sq,ht,i,lik,_volup,_volmin;

  _volup=1e+6; _volmin=1e-6;
  
  res = yt - b[1];
  res_sq = res.^2;
   
  
ht=zeros(T,1);
ht[1]=stdc(res)^2;

for i(2,T,1);
  ht[i]= b[2] + b[3]*res_sq[i-1] + b[4]*ht[i-1];
    if ht[i]>_volup;
      ht[i]=_volup;
    endif;
    if ht[i]<_volmin;
       ht[i]=_volmin;
    endif; 
endfor;

  
lik = -1/2*ln(2*pi) - 1/2*sumc(ln(ht)) - 1/2*sumc(res_sq./ht);

  retp(lik);
endp;