new;

/*
**            PROGRAM FOR GARCH(p,q)
**
**                 WRITTEN BY
**               SEUNG CHAN AHN
**           DEPARTMENT OF ECONOMICS
**             COLLEGE OF BUSINESS
**               TEMPE, AZ 85287
**
*/

@   OPEN OUTPUT FILE  @

/* ---------------------------------------- */

    output file = ahngarch.out reset ;

/* ---------------------------------------- */

@   DATA LOADING AND TRANSFORMATION  @

/* ------------------------------------------------------- */

    load dat[301,1] = exdmdo.db;
    y = 100*ln(dat[2:rows(dat)]./dat[1:rows(dat)-1]) ;

/* ------------------------------------------------------- */

@   DEPENDENT VARIABLE AND EXOGENOUS REGRESSORS @
@   Z: VARIABLES WHICH MAY APPEAR IN THE GARCH TERM @

/* ------------------------------------------------------- */

    yy = y ;
    xx = ones(rows(y),1)   ;
    zz = ones(rows(y),1)~y ;

/* ------------------------------------------------------- */

@   CONDITIONAL MEAN SPECIFICATION  @

proc res(mb)  ;
local ms      ;

/* ------------------------------------------------------- */

    ms = yy - xx*mb ;

/* ------------------------------------------------------- */

retp(ms)  ;
endp      ;

@ DEFINE K, P AND Q FOR GARCH @

/* ------------------ */

    k = cols(xx) ;
    p = 1 ;
    q = 1 ;

/* ------------------ */

@ ARCH-M @
@ IP1 = 0 IF NO ARCH-M ; IP1 = 1 IF ARCH-M @
@ IP2 = 0 IF ARCH-M WITH DEV; IP2 = 1 IF WITH VAR @

/* ------------------------------- */

    ip1 = 0 ;
    ip2 = 0 ;

/* ------------------------------- */

@ SPECIFY GARCH TERM  @

proc hterm(a,h,u2,tar,z) ;
local hc,pp,qq ;

    pp = p+1 ; qq = q+1 ;

/* ------------------------------------------------------------- */

   hc = a[1] + a[2]*u2[qq-1] + a[3]*h[pp-1] ;

   @ TARCH(0,1)  hc = a[1] + a[2]*u2[qq-1]
                       + a[3]*tar[qq-1]*u2[qq-1] ; @

   @ IGARCH(1,1) hc = a[1] + a[2]*u2[qq-1]
                       + (1-a[2])*h[pp-1] ;  @

   @ GARCH(1,1)  hc = a[1] + a[2]*u2[qq-1] + a[3]*h[pp-1]
                       + a[4]*h[pp-1]*zz[1] ;      @

/* ------------------------------------------------------------- */

retp(hc) ;
endp    ;

@ INITIAL VALUES FOR PARAMETER IN YOUR MODEL @

/* ----------------------------------------------------------- */

   inib1 = invpd(xx'xx)*(xx'yy) ; @ for mean @
   inib2 = {0.1, 0.2, 0.7 }     ; @ for conditional variances  @
   inib3 = {0.1}                ; @ for ARCH-M term (optimal) @
   inib = inib1|inib2 ;
   @ if ARCH-M, use inib = inib1|inib2|inib3 ; @

/* ------------------------------------------------------------ */

@   CONTROL ALGORITHM  @

/* -------------------------------------------------------- */

    algo = 2 ; @ 2=BFGS; 3=DFP; 4=NEWTON; 5=BHHH @
    lser = 2 ; @1=one; 2=STEPBT; 4=BRENT; 5=BHHHSTEP @
    gtol = 1e-3 ;

/* -------------------------------------------------------- */

@   GIVE TITLE FOR YOUR OUTPUT  @

/* ------------------------------ */

    qqt = "GARCH(1,1)" ;

/* ------------------------------ */

@   Max. Lag for LM tests  @

/* ------------------------------- */

    maxlag = 10 ;

/* ------------------------------- */



                     /**************************
                      *  Do not change below   *
                      **************************/

@   GARCH Specification   @
@   The GARCH-M term is the last entry of vb   @

    if p == 0   ; wp = 1 ; else ; wp = p ; endif ;

proc (3) = garch(mb,vb,p,q) ;
local uncm,n,u,uu,hh,i,t,th,j,c,du,dh,tarch ;

    uncm  = meanc(res(mb)^2) ;
    n     = rows(yy)         ;

    u          = sqrt(uncm)*ones(n+q,1) ;
    u[q+1:q+n] = res(mb)                ;
    hh         = uncm*ones(n+p,1)       ;

    i  = 1 ; do while i <= n ;

        dh = hh[i:i+wp-1] ;
        du = u[i:i+q-1]   ;
        tarch = dummydn(du,0,2) ;
        hh[p+i] = hterm(vb,dh,du^2,tarch,zz[i,.]) ;

        th = hh[p+i]      ;
        c  = vb[rows(vb)] ;
        u[q+i] = u[q+i] - ip1*(ip2*c*th + (1-ip2)*c*sqrt(th)) ;

    i = i + 1 ; endo ;

    hh = hh[p+1:rows(hh)]  ;
    uu = u[q+1:rows(u)]^2  ;

retp(hh,uu,u[q+1:rows(u)]) ;
endp        ;

library maxlik;
#include maxlik.ext ;
maxset ;

@   Define log-likelihood function   @

cc = rows(inib2) ;
proc func(b,yy)  ;
local hh,uu,u,logl ;

  b[k+1:k+cc] = abs(b[k+1:k+cc]) ;

  {hh,uu,u} = garch(b[1:k],b[k+1:k+cc],p,q)    ;
  logl      = -0.5*(ln(hh) + uu./hh)-0.9189385 ;

retp(logl) ;
endp       ;

    b0 = inib ;

@ Calculate the maximum likelihood estimates @

_max_CovPar = 3        ;
_max_algorithm = algo  ;
_max_LineSearch = lser ;
_max_GradTol = gtol    ;
__output = 10          ;

{b,f1,f2,f3,retcode }= maxlik(yy,0,&func,b0) ;

b[k+1:k+cc] = abs(b[k+1:k+cc]) ;
n  = rows(yy) ;
rc = (n^2)*(_max_finalHess)*(_max_XprodCov)*(_max_finalHess) ;
rc = invpd(rc)      ;
uc = _max_XprodCov  ;

ww = rows(b) ;
{hh,uu,u} = garch(b[1:k],b[k+1:ww],p,q) ;

@  Computing QMLE using approximated Hessian @

proc dmu(bb)  ;
local mm      ;

    mm = res(bb[1:k]) ;

retp ( mm ) ;
endp ;

proc dhh(bb)      ;
local mm1,mm2,mm3 ;

  {mm1,mm2,mm3} = garch(bb[1:k],bb[k+1:k+cc],p,q)  ;

retp (mm1) ;
endp       ;

dm = gradp(&dmu,b)./sqrt(hh)  ;
dh = gradp(&dhh,b)./hh        ;
apphess = dm'dm + 0.5*(dh'dh) ;
rdd = (apphess)*(_max_XprodCov)*(apphess) ;
rdd =invpd(rdd)                      ;


/*
**  Output
*/

format /rd 8,3;
print " ";
print " ";
print "S. C. AHN ---   " qqt "";
print " " ;
print "=========================================";
print "log-likelihood =" n*f1 ;
print "" ;
print "=========================================";

st1=sqrt(diag(uc))  ;
st2=sqrt(diag(rc))  ;
st3=sqrt(diag(rdd)) ;

print "Results with Usual Covariance Matrix" ;
print "" ;
print "For mean Specification" ;
print " Estimate std. err. t-stat" ;
print b[1:k]~st1[1:k]~(b[1:k]./st1[1:k]) ;

print "";

print "For GARCH Specification" ;
print " Estimate std. err. t-stat" ;
print b[k+1:ww]~st1[k+1:ww]~(b[k+1:ww]./st1[k+1:ww]) ;

print "";

print "=========================================";
print "Results with Robust Covariance Matrix" ;
print "" ;
print "For mean Specification" ;
print " Estimate std. err. t-stat" ;
print b[1:k]~st2[1:k]~(b[1:k]./st2[1:k]) ;

print "";

print "For GARCH Specification" ;
print " Estimate std. err. t-stat" ;
print b[k+1:ww]~st2[k+1:ww]~(b[k+1:ww]./st2[k+1:ww]) ;

print "";

print "=========================================";
print "Results with Robust Covariance Matrix by Approx. Hessian" ;
print "" ;
print "For mean Specification" ;
print " Estimate std. err. t-stat" ;
print b[1:k]~st3[1:k]~(b[1:k]./st3[1:k]) ;

print "";

print "For GARCH Specification" ;
print " Estimate std. err. t-stat" ;
print b[k+1:ww]~st3[k+1:ww]~(b[k+1:ww]./st3[k+1:ww]) ;

print "";


/*
**  LM Test for ARCH effects
*/

proc lmt(r) ;
local y,x,r2,b,i,t,f,lm ;

    t = rows(uu) - r ;
    y = uu[r+1:r+t]  ;
    x = zeros(rows(y),r) ;

    i = 1 ; do while i <= r ;

        x[.,i] = uu[i:i+t-1] ;

    i = i + 1 ; endo ;

    x  = ones(t,1)~x      ;
    b  = invpd(x'x)*(x'y) ;
    f  = x*b              ;
    lm = t*(f'f)/(y'y)    ;

retp(lm) ;
endp     ;

print "=========================================";
print "LM TESTS";
print "" ;

    i = 1 ; do while i <= maxlag ;


    "LM for lag=" i ", P =" lmt(i) cdfchic(lmt(i),i) ;
    "" ;

    i = i + 1 ; endo ;

/*
**  Tests Based on Standardized Residuals
*/

uh  = u./sqrt(hh);
u2h = uu./hh     ;
m1  = meanc(uh)  ;
m2  = meanc(u2h) ;
m3  = meanc(uh.*u2h)  ;
m4  = meanc(u2h.*u2h) ;
b3  = m3/(m2^1.5) ;
b4  = m4/(m2*m2)  ;
bera=n*( (b3^2)/6 + ((b4-3)^2)/24 ) ;
print " ";
print "=====================================" ;
print " Skewness =" b3;
print " Kurtosis =" b4;
print " Bera-Jarque Normality (df=2), p =" bera cdfchic(bera,2) ;
"";

print "=====================================" ;
ncor=20;
cr=ones(ncor,1);
q=0.0;
uh1=uh-m1;
il=1;
do while il<=ncor;
     uh1=0.0|uh1[1:n-1];
     cr[il,1]=meanc((uh-m1).*uh1)/meanc((uh-m1).*(uh-m1));
     q=q+cr[il,1]*cr[il,1]/(n-il);
     il=il+1;
endo;
q=q*(n+2)*n;
print "  Q(20) for the levels (df=20), p =" q cdfchic(q,20) ;
print "  Lag 1,...,20:";
print cr[1:5]'    ;
print cr[6:10]'   ;
print cr[11:15]'  ;
print cr[16:20]'  ;
print "";

q=0.0;
uh1=u2h-m2;
il=1;
do while il<=ncor;
     uh1=0.0|uh1[1:n-1];
     cr[il,1]=meanc((u2h-m2).*uh1)/meanc((u2h-m2).*(u2h-m2));
     q=q+cr[il,1]*cr[il,1]/(n-il);
     il=il+1;
endo;
q=q*(n+2)*n;
print "  Q(20) for the squares (df=20), p =" q cdfchic(q,20) ;
print "  Lag 1,...,20:";
print cr[1:5]'   ;
print cr[6:10]'  ;
print cr[11:15]' ;
print cr[16:20]' ;
print "";
print "";

output off ;