// simulate a GARCH model and analyze the data

#include <oxstd.h>
#include <oxfloat.h>

#import <maximize>


main()
{
     decl d,x,r,rdata,h,n,st,i,eps,q,acfs,acfd,file;

     d=loadmat("vwretd1.mat");
     
     rdata=100*d[5000:][3];

// model parameters 
// r_t = mu + rho*r_{t-1} + esp_t
// eps_t = sqrt(h_t)*z_t, z_t~N(0,1)
// h_t = omega + alpha*eps_{t-1}^2 + beta*h_{t-1}

// parameter order: mu,rho,omega,alpha,beta

x=<     0.053896;
      0.17838;
    0.0063616;
     0.089330;
      0.90533>;

//      x[4]=.91;

      print("model parameters ",x);
      print("Persistence alpha+beta= ",x[3]+x[4]);
      
      
// total simulations of the model
      n=100000;

// start up values that will be discarded
      st=100;

// initialize vectors that will hold simulated data

      r=zeros(n+st,1);
      h=zeros(n+st,1);
      eps=zeros(n+st,1);
      
// set startup conditions
      r[0]=0.d0;
      h[0]=1.2d0;
      eps[0]=sqrt(h[0])*rann(1,1);
      
// simulate model
	   
      for(i=1;i<n+st;i++)
      {
	   h[i]=x[2]+x[3]*eps[i-1]^2+x[4]*h[i-1];
	   eps[i]=sqrt(h[i])*rann(1,1);
	   r[i]=x[0]+x[1]*r[i-1]+eps[i];
      }

      // compare moments from simulated data and real data
      print("\nobs, mean, stdev, skew ,kurt\n");
      print("\nSimulated GARCH(1,1) data\n",moments(r[st:],4)');

      print("\n Data\n",moments(rdata,4)',"\n");
     
      q=<.01,.05,.10,.5,.90,.95,.99>;
      print("      Quantiles    sim data     CRSP      normal distr\n");
      
      print(q'~quantilec(r[st:],q)~quantilec(rdata,q)~quann(q)');

// compare ACF of abs value of returns
      acfs=acf(fabs(r[st:]),1000);
      acfd=acf(fabs(rdata),1000);
 

      file=fopen("sacf.dat","w");
      fprint(file,acfs[1:]~acfd[1:]);



     
}