options double crt;
name bg11 'GARCH(1,1) benchmark from Bollerslev-Ghysels JBES 4/96';
?  (using MLNEWT and computing White SEs)
?  by Clint Cummins 2/1/98
?
?  Below we reproduce p.146 Table 2 column 4 - GARCH(1,1) results
?          B-G     TSP MLNEWT1  TSP MLNEWT2  TSP MLNEWT3
?  mu    -.006     -.00617319   -.00619041   -.00626932
?        (.009)    (.00920255)  (.00918935)  (.00923712)
?  omega  .264      .263167      .263165      .241437
?        (.075)    (.074486)    (.074483)    (.053175)
?  alpha0 .0108     .010761      .010761      .010983
?                  (.00649317)  (.00649319)  (.00657626)
?  alpha1 .153      .153132      .153134      .148700
?        (.054)    (.053533)    (.053532)    (.054033)
?  beta1  .806      .805977      .805974      .805809
?        (.073)    (.072462)    (.072461)    (.073853)
?  LogL  -1106.6  -1106.60665  -1106.60788  -1106.94851
?--------------------
?  Notes:
?  Estimators vary according to handling of initial conditions:
?  MLNEWT1:  h(0) = e(0)**2 = SSR(OLS)/T
?  MLNEWT2:  h(0) = e(0)**2 = SSR(mu_current_iteration)/T
?  MLNEWT3:  h(0) = e(0)**2 = omega (see definition below)
?  all SEs are robust (QMLE - White)
?  SE for omega in TSP is computed with delta method (analyz)
?  omega = unconditional E(e**2)
?        = alpha0/(1-alpha1-beta1)
?  B-G LogL is computed from AIC or SIC in Table 2
?  (they agree to 5 digits)
?  Parameter names used here (g11s.tsp):
?    alpha0 = sig0
?    mu = gamma0
?
?  y = Deutschmark-British Pound daily percentage nominal returns,
?      1/3/1984 - 12/31/1991
?      Obtained on 1/30/98 from   ftp  web.amstat.org
?      jbes/View/96-2-APR/bollerslev_ghysels/bollerslev.sec41.dat
smpl 1,1974;
read(file='dmbp.dat') y monday;
olsq y c;  ? table 2 column 1 - easy to reproduce

const ifhess,1;  ? use 2nd derivatives
input 'g11s.tsp';   ? define model logl, derivatives, global variables

frml eomega omega = sig0/(1-alpha1-beta1);
set vcovw = 0;  ? global White vcov for analyz from mlnewt

do init=1,3;
  param gamma0,-.006 omega,.264 sig0 alpha1,.153 beta1,.806;
  set sig0 = omega*(1-alpha1-beta1);
  initd;
  mlnewt 30 1e-7;    ? arguments are maxit and tol
  analyz(vcov=vcovw,names=bnames) eomega;
enddo;

?input 'mlbhhh.tsp';   ? BHHH iterations
input 'mlnewt.tsp';   ? Newton iterations
input 'gradchec.tsp';   ? gradient check
input 'hesschec.tsp';   ? hessian check