options crt nodate double;
name garchma 'GARCH(1,1) with MA(1) error term';
const nob,100;
smpl 1,nob;

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

supres smpl;
param b0,0 b1,0 a0,1 a1,.4 beta1,.3 h0,1 theta,0;
? new ML PROC method
ML(maxit=50) GMA b0 b1 a0 a1 beta1 theta h0;   ? 21 iter, 489 feval

Proc GMA;

  ? check parameter constraints
  if ( abs(beta1) > exp(6/@nob) ) .or.    ? so that beta1**@nob <= exp(6)
                                          ? to prevent overflow in h
     ( h0 <= 0 ) .or.      ? bad h0
     ( abs(theta) > 1 ) ;  ? invertibility for MA
     then; goto 10;

  smpl 2,nob;
  u = 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.
  eps = 0;  ? set initial residual to zero (could estimate eps0)
  smpl 3,nob;
  eps = u - theta*eps(-1);   ? recursion for eps in terms of u for MA(1)
  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;