options crt nodate double;
name aparch 'Asymmetric Power GARCH(1,1) with MA(1) residual';
?  Following Ding, Granger, and Engle, Journal of Empirical Finance,
?  1993 pp.83-106  (here we *do* estimate the power delta, starting
?  from a true value of 1.5).

const nob,100;
smpl 1,nob;

? generate artificial data from true model
random(seedin=189741) x e;
s = 1;
eps = 0;
smpl 2,nob;
param a0,1 a1,.4 gamma1,.2 beta1,.3 theta,.5 delta,1.5;
s = ( a0 + a1*s(-1)*[abs(e(-1)) - gamma1*e(-1)]**delta
       + beta1*s(-1)**delta ) ** (1/delta) ;  ?  APARCH(1,1) std.dev.
eps = e*s;
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 delta,1.6;
? new ML PROC method
ML(maxit=50) APARCH b0 b1 a0 a1 gamma1 beta1 theta delta;  ? 30 iter, 749 feval

Proc APARCH;

  ? 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
     ( delta  <=  0 ) .or.      ? bad h0
     ( abs(theta) > 1 ) ;  ? invertibility for MA
     then; goto 10;

  smpl 2,nob;
  u = y - b0 - b1*x;      ? residual
  s = (a0/h0den)**(1/delta);  ? instead of estimating h0,
                              ? we use the expected value
                              ? a0/(1-(a1*(1+gamma1**2)+beta1)) .
                              ? (not sure if this formula really carries
                              ?  over from AGARCH model...)
  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)
  s = ( a0 + a1*[abs(eps(-1)) - gamma1*eps(-1)]**delta
           + beta1*s(-1)**delta ) ** (1/delta);  ?  std.dev. of residual

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

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

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