options crt nodate double;
name archml 'various estimation methods on artificial data';
const nob,100;
smpl 1,nob;

random(seedin=189741) x e;
h = 1;
smpl 2,nob;
param a0,1 a1,.4;
h = a0 + a1*h(-1)*e(-1)**2 ;  ?  ARCH(1)
eps = e*sqrt(h);
y = 1.2 + 1.4*x + eps;
arch(nar=1) y c x;  ? 11 iter, 26 feval

frml arch1 logl = -log(h)/2 + lnorm(eps/sqrt(h));
frml eqh h = a0 + a1*eps(-1)**2;
frml eqe eps = y - b0 - b1*x;  param b0,b1;
eqsub arch1 eqh eqe;
smpl 3,nob;
ml arch1;  ? 10 iter, 36 feval

? compare w/ hiter=n,hcov=n, using same starting values
param a0,1 a1,.4 b0,0 b1,0;
ml(hiter=n,hcov=n) arch1;  ? 8 iter, 25 feval

? compare w/ hiter=f,hcov=f, using same starting values
param a0,1 a1,.4 b0,0 b1,0;
ml(hiter=f,hcov=f) arch1;  ? 15 iter, 51 feval

? compare w/ hiter=f,hcov=f,grad=c2 using same starting values
param a0,1 a1,.4 b0,0 b1,0;
ml(hiter=f,hcov=f,grad=c2) arch1;  ? 10 iter, 148 feval (to bad LogL)
? again from "converged" values
ml(hiter=f,hcov=f,grad=c2) arch1;  ? 8 iter, 123 feval to right LogL

? PROC version
proc af;
  smpl 2,nob;
  eps = y - b0 - b1*x;      ? residual
  smpl 3,nob;
  h = a0 + a1*eps(-1)**2;   ? variance of residual
  ? check positive variance constraint
  mat @min = min(h);   ? don't use msd; it can overflow
  if miss(@min) .or. (@min <= 0); then; do;
    set @logl = @miss;
    goto 10;
  enddo;
  logl = -log(h)/2 + lnorm(eps/sqrt(h));
  mat @logl = sum(logl);
  10 nodbug;
endproc;

supres smpl;
param a0,1 a1,.4 b0,0 b1,0;
ml af a0 a1 b0 b1;   ? 12 iter, 197 feval to right LogL