options crt nodate double;
name agarch 'Asymmetric GARCH(1,1) with MA(1) residual';
?  Following Ding, Granger, and Engle, Journal of Empirical Finance,
?  1993 pp.83-106  (here we do not estimate the power delta; leave it
?  fixed at 2).

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 gamma1,.2 beta1,.3 theta,.5;
h = a0 + a1*h(-1)*[abs(e(-1)) - gamma1*e(-1)]**2  ? asymmetric part with gamma
       + beta1*h(-1);  ?  AGARCH(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 gamma1,.2 beta1,.3 theta,0;
? new ML PROC method
ML(maxit=50) AGARCH b0 b1 a0 a1 gamma1 beta1 theta;   ? 19 iter, 419 feval

Proc AGARCH;

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

  smpl 2,nob;
  u = y - b0 - b1*x;      ? residual
  h = a0/h0den;  ? instead of estimating h0, we use the expected value
                 ? a0/(1-(a1*(1+gamma1**2)+beta1)) .
  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*[abs(eps(-1)) - gamma1*eps(-1)]**2  ? asymmetric part with gamma
       + 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;