options crt nodate double;  ? double is used for precise residuals in the Proc
name garchml 'various estimation methods on artificial data';
const nob,100;
smpl 1,nob;

? generate artificial data from true model
random(seedin=189741) x e;
h = 1;
smpl 2,nob;
param a0,1 a1,.4 beta1,.3;
h = a0 + a1*h(-1)*e(-1)**2 + beta1*h(-1);  ?  GARCH(1,1)
eps = e*sqrt(h);
y = 1.2 + 1.4*x + eps;

? check against arch command, to reproduce its results
arch(nar=1,nma=1,maxit=50) y c x;                  ? 30 iter, 120 feval
? new hiter option
arch(nar=1,nma=1,maxit=50,hiter=f,hcov=f) y c x;   ? 20 iter,  63 feval

supres smpl;
param b0,0 b1,0 a0,1 a1,.4 beta1,.3 h0,1;
? new ML PROC method
ML GF b0 b1 a0 a1 beta1 h0;   ? 16 iter, 345 feval

Proc GF;

  ? check for large beta1 value (can cause overflow)
  ? or bad h0
  if ( abs(beta1) > exp(6/@nob) ) .or.    ? so that beta1**@nob <= exp(6)
     ( h0 <= 0        ) ; then; goto 10;

  smpl 2,nob;
  eps = y - b0 - b1*x;      ? residual
  h = h0;  ? instead of estimating h0, we could use an expected value
           ? like  a0/(1-(a1+beta1)) .  But h0 is estimated here to
           ? reproduce results from the arch command
  smpl 3,nob;
  h = a0 + a1*eps(-1)**2 + beta1*h(-1);   ? variance of residual

  ? check positive variance constraint
  ?msd(noprint,terse) h;   ? do not use MSD; it can overflow
  mat @min = min(h);
  if miss(@min) .or. (@min <= 0); then; goto 10;

  logl = -log(h)/2 + lnorm(eps/sqrt(h));
  mat @logl = sum(logl);
  goto 20;

  10 nodbug;  ? arithmetic error
     set @logl = @miss;
  20 nodbug;
endproc;