? GARCH(1,1) with unconditional presample h0 and e0**2
? 1st and 2nd derivatives
?  mostly following O. Linton, "Econometric Theory", 8/97
?  and Fiorentini, Calzolari, and Panattoni, JAE 96
?  h0 = e0**2 = sig0/(1-alpha1-beta1)   or   SSR/T
? by Clint Cummins 1/98

? Code below uses Procs  MLBHHH/MLNEWT to estimate, and GRADCHEC to check
? analytic derivatives.
? Define Procs  DOLOGL  and  DERIVS , and define some global variables
? (shared by these Procs).

options double;  ? make sure recursive variables are in double precision

supres @smpl;
copy @smpl smpall;   ? includes initial obs
copy @smpl smpw;     ? smpl for LogL
copy @smpl smp2;     ? smpl for constructing e and derivatives recursively
set smp2(1) = smp2(1)+1;
smpl 1,1;
copy @smpl smp0;  ? first obs

list bnames sig0 beta1 alpha1 gamma0;
list dnames dlds dldb dlda dldg;
const h0;  ? h0 is a global variable

? compute first and second derivatives at true parameter values,
? but without knowing htrue and etrue
list res h u e rh;
list dh1 dhds dhdb dhda dhdg;
list dh2 d2hdss d2hdsb d2hdsa d2hdsg   ? d2hdsb = d2h/(ds*db), etc.
                d2hdbb d2hdba d2hdbg
                       d2hdaa d2hdag
                              d2hdgg;
list dl2 d2ldss d2ldsb d2ldsa d2ldsg
                d2ldbb d2ldba d2ldbg
                       d2ldaa d2ldag
                              d2ldgg;
smpl smpall;
? init all variables to full length zeros, to make them global variables
dot res dh1 dnames;
   . = 0;
enddot;
if ifhess; then; do;
dot dh2 dl2;
   . = 0;
enddot;
enddo;
set h0 = 0;  ? make h0 global to pass from dologl to derivs
set sig2uols = 0;  ? make global
? initialize derivs w.r.t. h0 for all INITs
dot dh0dg d2h0dgg dh0ds dh0db d2h0dsb d2h0dbb;
  set . = 0;
enddot;
?------------
Proc initd;  ? (re)initialize derivatives
  dot dh0dg d2h0dgg dh0ds dh0db d2h0dsb d2h0dbb;
    set . = 0;
  enddot;
  ? run OLS if INIT=1
  if (INIT=1); then; do;
    olsq(silent) y c;
    set sig2uols = @SSR/@NOB;
    print sig2uols;
  enddo;
  smpl smpall;
  ? init all variables to full length zeros
  dot res dh1 dnames;
     . = 0;
  enddot;
  if ifhess; then; do;
    dot dh2 dl2;
       . = 0;
    enddot;
  enddo;
endproc;

?--------------------------------------------------------------
Proc DOLOGL logl;

? check constraints
if (alpha1+beta1)>=1 .or.
   sig0 <= 0 .or.
   alpha1 <= 0 .or.
   beta1 <= 0
   ; then; goto 800;

smpl smpall;
u = y - gamma0;

if (INIT=1); then;
  set h0 = sig2uols;
if (INIT=2); then;
  mat h0 = (u'u)/@NOB;
if (INIT=3); then;
  set h0 = sig0/(1-alpha1-beta1);

smpl smp0;   ? first obs
h = sig0 + alpha1*h0 + beta1*h0;
smpl smp2;
h = sig0 + alpha1*u(-1)**2 + beta1*h(-1);
smpl smpall;

? check for h<=0  or  missing  (numerical error)
badh = .not.( (.not.miss(h))*(h>0) );
msd(noprint) badh;
if (@sum>0); then; do;
800  set logl = @miss;
     goto 900;
enddo;

rh = sqrt(h);
e = u/rh;
logli = -log(h)/2 + lnorm(e);

mat logl = sum(logli);

900 nodbug;
endproc;

?--------------------------------------------------------------
Proc DERIVS;
?  derivatives of LogL

?if (INIT=1); ? dh0/d(any) = 0
if (INIT=2); then; do;
  smpl smpall;
  mat dh0dg = -sum(u)*2/@nob;
  mat d2h0dgg = 2;
enddo;
if (INIT=3); then; do;
  set dh0ds = 1/(1 - alpha1 - beta1);
  set dh0db = h0*dh0ds;
  if ifhess; then; do;
    set d2h0dsb = dh0ds**2;
    set d2h0dbb = 2*h0*d2h0dsb;
  enddo;
enddo;
smpl smp0;  ? first obs
dhds = 1 + (alpha1+beta1)*dh0ds;
dhdb = h0 + (alpha1+beta1)*dh0db;
dhda = dhdb;
dhdg = (alpha1+beta1)*dh0dg;
if ifhess; then; do;
  d2hdss = 0;
  d2hdsb = dh0ds + (alpha1+beta1)*d2h0dsb;
  d2hdsa = d2hdsb;
  d2hdsg = 0;
  d2hdbb = 2*dh0db + (alpha1+beta1)*d2h0dbb;
  d2hdba = d2hdbb;
  d2hdbg = dh0dg;
  d2hdaa = d2hdbb;
  d2hdag = dh0dg;
  d2hdgg = (alpha1+beta1)*d2h0dgg;
enddo;

smpl smp2;  ? later obs - recursive calcs
dhds = beta1*dhds(-1) + 1;
dhdb = beta1*dhdb(-1) + h(-1);
dhda = beta1*dhda(-1) + u(-1)**2;
dhdg = beta1*dhdg(-1) - 2*alpha1*u(-1);
if ifhess; then; do;
?d2hdss = beta1*dd2hdss(-1) = 0;
d2hdss = 0;
d2hdsb = beta1*d2hdsb(-1) + dhds(-1);
d2hdsa = beta1*d2hdsa(-1);
?d2hdsg = beta1*d2hdsg(-1) = 0;
d2hdsg = 0;
d2hdbb = beta1*d2hdbb(-1) + 2*dhdb(-1);
d2hdba = beta1*d2hdba(-1) + dhda(-1);
d2hdbg = beta1*d2hdbg(-1) + dhdg(-1);
d2hdaa = beta1*d2hdaa(-1);
d2hdag = beta1*d2hdag(-1) - 2*u(-1);
d2hdgg = beta1*d2hdgg(-1) + 2*alpha1;
enddo;

? transform to log derivatives  (but retain same names)
?  g = log(h)
?  dgdx = dhdx/h
?  d2gdxy = d2hdxy/h - dgdx*dgdy
smpl smpall;  ? all obs
dot dh1;
  . = ./h;
enddot;
if ifhess; then; do;
d2hdss = d2hdss/h - dhds*dhds;
d2hdsb = d2hdsb/h - dhds*dhdb;
d2hdsa = d2hdsa/h - dhds*dhda;
d2hdsg = d2hdsg/h - dhds*dhdg;
d2hdbb = d2hdbb/h - dhdb*dhdb;
d2hdba = d2hdba/h - dhdb*dhda;
d2hdbg = d2hdbg/h - dhdb*dhdg;
d2hdaa = d2hdaa/h - dhda*dhda;
d2hdag = d2hdag/h - dhda*dhdg;
d2hdgg = d2hdgg/h - dhdg*dhdg;
enddo;

e2 = e**2;
e2o2 = e2/2;
e2m1o2 = (e2 - 1)/2;
eorh = e/rh;

? transform to log likelihood derivatives
dot s b a g;
  dld. = e2m1o2*dhd.;
enddot;
dldg = dldg + eorh;
if ifhess; then; do;
d2ldss = -e2m1o2*d2hdss + e2o2*dhds*dhds;
d2ldsb = -e2m1o2*d2hdsb + e2o2*dhds*dhdb;
d2ldsa = -e2m1o2*d2hdsa + e2o2*dhds*dhda;
d2ldsg = -e2m1o2*d2hdsg + e2o2*dhds*dhdg + eorh*dhds;
d2ldbb = -e2m1o2*d2hdbb + e2o2*dhdb*dhdb;
d2ldba = -e2m1o2*d2hdba + e2o2*dhdb*dhda;
d2ldbg = -e2m1o2*d2hdbg + e2o2*dhdb*dhdg + eorh*dhdb;
d2ldaa = -e2m1o2*d2hdaa + e2o2*dhda*dhda;
d2ldag = -e2m1o2*d2hdag + e2o2*dhda*dhdg + eorh*dhda;
d2ldgg = -e2m1o2*d2hdgg + e2o2*dhdg*dhdg + 2*eorh*dhdg + 1/h;
enddo;
endproc;