new;
library maxlik;
#include maxlik.ext ;
maxset ;


/*   GARCH (1,1)  TIM'S DATA SET
                                                            */

load p[3006,1]=exdmdo.txt;
n=rows(p);
p=100.0*ln(p);
y=p[2:n,1]-p[1:n-1,1];
n=rows(y);


/*                                                             */
/*   This function sets up the likelihood function.            */
/*                                                             */

proc func(b,y);
     local e,h,h1,y1,unc,e2,e2m1,logl;
     e=y - b[1,1]; @ error terms @
     e2=e.*e;
     unc=meanc(e2);
     e2m1=unc|e2[1:n-1,1] ;
     b[2,1]=abs(b[2,1]);
     b[3,1]=abs(b[3,1]);
     b[4,1]=abs(b[4,1]);
     h1=unc;
     h=recserar(b[2,1]+b[3,1]*e2m1,h1,b[4,1]);
     logl=-0.5*(ln(h) + e2./h)-0.9189385;
     retp(logl);
endp;

/*                                                             */
/*   Initialize the parameters                                 */
/*                                                             */

let b0[4,1] = -0.004  0.03  0.07 0.89 ;

/*   Calculate the maximum likelihood estimates                */
/*                                                             */

_max_CovPar = 3        ; 
    @ 1 = cov by hess; 2 = cov by bhhh; 3 = robust cov @
_max_algorithm = 2     ; 
    @ 2=BFGS; 3=DFP; 4=NEWTON; 5=BHHH @
_max_LineSearch = 2    ; 
    @1=one; 2=STEPBT; 4=BRENT; 5=BHHHSTEP @
_max_GradTol = 1e-3    ;
__output = 10          ;

{b,f1,f2,f3, retcode }=maxlik(y,0,&func,b0);
@ b = est. b; f1 = max. logl/# of obs; f2 = grad; f3 = cov @


format /rd 13,5;
st1=sqrt( diag(f3) ) ; @ s.e. from robust cov @
st2=sqrt( diag(_max_HessCov) ) ; @ s.e. from hess @
st3=sqrt( diag(_max_XprodCov) ) ; @ s.e. from bhhh @

/*                                                             */
/*   Output                                                    */
/*                                                             */

output file = timgarch_1.out  reset;        /*Output*/

format /rd 11,4;
print " ";
print " ";
print " ";
print "        TIM  ---     GARCH(1,1) Model" ;

tsthess=b./st2; 
tstrob=b./st1;
print "" ;
print "     Estimate   SE (ROB.)   t-St.(Rob.) SE(Hess.)  t-St.(Hess.)" ;
print b~st1~tstrob~st2~tsthess ;
print " ";
print " ";
print "Obs.:" n;

e=y - b[1,1]; @ error terms @
e2=e.*e;
unc=meanc(e2);
e2m1=unc|e2[1:n-1,1] ;
b[2,1]=abs(b[2,1]);
b[3,1]=abs(b[3,1]);
b[4,1]=abs(b[4,1]);
h1=unc;
h=recserar(b[2,1]+b[3,1]*e2m1,h1,b[4,1]);
logl=-0.5*(ln(h) + e2./h)-0.9189385;
PRINT " ";
PRINT "logl:"  sumc(logl);

eh=e./sqrt(h);
e2h=e2./h;
m1=meanc(eh);
m2=meanc(e2h);
m3=meanc(eh.*e2h);
m4=meanc(e2h.*e2h);
b3=m3/(m2^1.5);
b4=m4/(m2*m2);
print " ";
print " Skewness :" b3;
print " Kurtosis :" b4;

ncor=20;
cr=ones(ncor,1);
q=0.0;
eh1=eh-m1;
il=1;
do while il<=ncor;
     eh1=0.0|eh1[1:n-1];
     cr[il,1]=meanc((eh-m1).*eh1)/meanc((eh-m1).*(eh-m1));
     q=q+cr[il,1]*cr[il,1]/(n-il);
     il=il+1;
endo;
q=q*(n+2)*n;
print "";
print "";
print "  Q(20) for the levels: " q;
print "Lag 1,...,10:";
print cr[1:10]';
print "";

q=0.0;
eh1=e2h-m2;
il=1;
do while il<=ncor;
     eh1=0.0|eh1[1:n-1];
     cr[il,1]=meanc((e2h-m2).*eh1)/meanc((e2h-m2).*(e2h-m2));
     q=q+cr[il,1]*cr[il,1]/(n-il);
     il=il+1;
endo;
q=q*(n+2)*n;
print "";
print " Q(20) for the squares:" q;
print " Lag 1,...,10:";
print cr[1:10]';
print "";
print "";

end;