options crt nodate double;
name garchm 'GARCH(1,1)-M -- ARCH and ML PROC commands';
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)  variance
eps = e*sqrt(h);
param theta,.5 lambda,.5;
y = 1.2 + 1.4*x + theta*h**lambda + eps;  ? M term in regression

? check against arch command, to reproduce its results
? default internal starting values
arch(nar=1,nma=1,mean,maxit=50) y c x;                  ? 25 iter, 98 feval

? common starting values
param b0,0 b1,0 theta,.1 a0,1 a1,.4 beta1,.3 h0,1;
mmake @start b0 b1 theta a0 a1 beta1 h0;
arch(nar=1,nma=1,mean,maxit=50) y c x;                  ? 5 iter, 51 feval
                                                        ? (fail to improve)

? new hiter option  (although using analytic derivatives)
arch(nar=1,nma=1,mean,maxit=50,hiter=f,hcov=f) y c x;   ? 29 iter, 85 feval

supres smpl;
param b0,0 b1,0 theta,.1 a0,1 a1,.4 beta1,.3 h0,1;

?  Set up for GARCH-M eval PROC by writing recursive FRMLs that
?  can be evaluated by SOLVE (first h, then eps, for each time period)
?  This contrasts with simple GARCH, where h and eps can be evaluated
?  independently with plain GENR commands.  In GARCH-M, h and eps are
?  related recursively.
frml eqh h = a0 + a1*eps(-1)**2 + beta1*h(-1);   ? variance of residual
frml eqe eps = y - (b0 + b1*x + theta*h**lambda);      ? residual
model eqh eqe mgm;

? new ML PROC method
ML(maxit=50) GARCHM b0 b1 theta a0 a1 beta1 h0;   ? 23 iter, 548 feval

Proc GARCHM;

  ? 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;
  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
  genr eqe;
  smpl 3,nob;
  solve(silent,noprnsim,tag=none) mgm;

  ? we could check @ifconv with the 2 commands below, but in this
  ? particular model (lambda=.5), it is sufficient to check for min(h) <= 0 .
  ?mat nconv = sum(@ifconv);
  ?if (nconv < @nob); then; goto 10;  ? solve might exit without
            ? converging in all period if h(t)<0, due to negative
            ? square root (or similar, when lambda is not an integer)

  ? 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;