// Estimate an ARCH(1) model by MLE and report diagnostic tests

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

// declare global variables that can be used in all subprograms 

decl r,h,npar,n,istart,iend,nobs,e;

// define arch loglikelihood function

ARCH1(x,f,score,hessian)
{
decl q,i,lgl,uncond;

 uncond=x[1]/(1-x[2]);
 
// initialize lgl function to 0
lgl=0.;

// construct lgl observation by observation, start at istart.
 for(i=istart;i<=iend;i++)
 {

// construct the conditional variance

      if(i==istart)
      {
	   // deal with startup issues (set presample values of var = unconditional values)
	   h[i]=x[1]+x[2]*uncond;
      }
      else
      {
	   h[i]=x[1]+x[2]*e[i-1]^2;
      }
      
// construct the conditional mean
      e[i]=r[i]-x[0];
      
      lgl=lgl -.5*log(M_2PI)-.5*log(h[i]) - .5*e[i]^2/h[i];
 }
 
// return the function (lgl) value
 f[0]=lgl;

 return 1;
}

// define Ljung-Box test for serial correlation Q(p), for data vector z

LB(z,p)
{
     decl rho,i,T,test,pvalue;
     
     T=rows(z);     
     rho=acf(z,p);
     test=0.;
     
     for(i=1;i<=p;i++)
     {
	  test=test+rho[i]^2/(T-i);
     }
     test=T*(T+2)*test;
     // get p-value
     pvalue=1.-probchi(test,p);
     print("\n",p," ",test," ",pvalue,"\n");
}



main()
{
     decl d,x,f,ir,NH,Hess,var,se,z;

     d=loadmat("vwretd1.mat");
     
     r=100*d[][3];
     print("total nobs = ",rows(r),"\n");
     
// start estimation at observation istart.     
     istart=5000;
// end estimation at iend
     iend=7000;
// number of parameters
     npar=3;
// total number of data observations
     n=rows(r);
// number of obs used in estimation
     nobs=iend-istart+1;
// initialize the cond var, and error vector
     h=zeros(n,1);
     e=zeros(n,1);
// set startup values for optimization
     x=<0.0; .512 ; .5 >;

// set some defaults on the optimizor
     MaxControl(100,10);
     MaxControlEps(.0001, .001);
     NH=unit(rows(x))*.001;
// call BFGS to find max of lgl
     ir=MaxBFGS(ARCH1, &x,&f,NH,1);
     MaxConvergenceMsg(ir);
     print(x,"lgl = ",f,"\n");

     
// calculate numerial hessian at the optimum.
     ir=Num2Derivative(ARCH1,x,&Hess);
// calc Var-Cov matrix estimate of the parameters
     var=invert(-Hess);
     print("\n\n"," Var-Cov = \n",var);
     print("\n obs = ",nobs,"\n\n");
     print(" lgl = ",f,"\n\n");

// obtain standard errors   
     se=sqrt(diagonal(var)');

// report estimates, standard errors and t-statistics
     print("       Estimates   SE           T-stats","\n");
     print(x~se~(x./se),"\n");

//construct standardized residuals
z=e[istart:iend]./sqrt(h[istart:iend]);

// calc moments of standardized residuals
print("\n"," Moments of raw data (returns)");
print("\n","obs, mean, stdev, skew ,kurt\n",moments(r[istart:iend],4)',"\n");
print("\n"," Moments of standadized residuals");
print("\n","obs, mean, stdev, skew ,kurt\n",moments(z,4)',"\n");


print(" LB test of stardardized residual\n\n");
print("k   Q(k)   pvalue\n");
LB(z,5);     
LB(z,20);

z=z.^2;
print("\n\n LB test of squared stardardized residual\n\n");
print("k   Q(k)   pvalue\n");
LB(z,5);     
LB(z,20);

// LB test on squared returns
print("\n\n LB test on squared returns\n\n");
print("k   Q(k)   pvalue\n");
LB(r[istart:iend].^2,5);     
LB(r[istart:iend].^2,20);


     
     
     
}