? MA(1) with 1st and 2nd derivatives
? Follows Davidson and MacKinnon, p.354-356.
? by Clint Cummins 1/2000

? Code below uses Procs  MLBHHH/MLNEWT to estimate, and GRADCHEC to check
? analytic derivatives.
? We have to transform the data to a simple FREQ/SMPL, define the
? Procs  DOLOGL  and  DERIVS , and define some global variables
? (shared by these Procs).

options double;  ? make sure recursive variables are in double precision

supres @smpl;
mmake my y;  ? copy dep var to a matrix
mat n = nrow(my);
copy @freq freqsav;  ? convert to freq n
copy @smpl smplsav;  ? (no plans to convert back, at present...)
freq n;
set np1 = n+1;
smpl 1,np1;
copy @smpl smpall;   ? includes presample obs
copy @smpl smpw;   ? includes presample obs
copy @smpl smp2;     ? smpl for true residuals
set smp2(1) = smp2(1)+1;
smpl smp2;
unmake my yy;
smpl 1,1;
copy @smpl smp0;  ? first obs

list bnames theta1 sigi;
list dnames dldt dlds;

? global scalar variables, shared between DoLogL and Derives
const e0 Jac;

? compute first and second derivatives at true parameter values,
? but without knowing htrue and etrue
list res e ys dysdt dzsdt dedt;
list de2 d2edtt d2ysdtt d2zsdtt;
list dl2 d2ldtt d2ldts
                d2ldss;

? init all variables to full length zeros, to make them global variables
?  (Same code as in Proc initd, but run it outside the Proc, so the
?   variables are made global).
  smpl smpall;
  ? init all variables to full length zeros
  dot res dnames;
     . = 0;
  enddot;
  zs = -1;
  if ifhess; then; do;
    dot de2 dl2;
       . = 0;
    enddot;
  enddo;

? a few more global scalar variables
set pi = 4*atan(1);
set loglcons =  -n*log(sqrt(2*pi));
?regopt(nopvcalc) dw;  ? to speed up OLS

?------------
Proc initd;  ? (re)initialize derivatives
  smpl smpall;
  ? init all variables to full length zeros
  dot res dnames;
     . = 0;
  enddot;
  zs = -1;
  if ifhess; then; do;
    dot de2 dl2;
       . = 0;
    enddot;
  enddo;
endproc;


?--------------------------------------------------------------
Proc DOLOGL logl;

? check constraints
if sigi <= 0 .or.
   abs(theta1) > 1   ? or we could set theta1=1 ...
   ; then; goto 800;

smpl smp2;
ys = yy + theta1*ys(-1);
zs = theta1*zs(-1);
smpl smpall;
?ols(silent) ys zs;
?e = @res;
?set e0 = @coef;
?  zs'zs = Jac = (1-theta1**(2*np1))/(1-theta1**2)  or  n+1
set Jac = (abs(theta1)<1)*[ (1-theta1**(2*np1))/(1-theta1**2) ]
        + (abs(theta1)=1)*[ np1 ] ;
mat zspys = zs'ys;
?mat e0 = (zs'zs)"(zs'ys);
set e0 = zspys/Jac;
e = ys - e0*zs;
mat @ssr = e'e;
set logl = n*log(sigi) - log(Jac)/2 + -@ssr*sigi**2/2 + loglcons;
goto 900;

800  set logl = @miss;

900 nodbug;
endproc;

?--------------------------------------------------------------
Proc DERIVS;
?  derivatives of LogL

? first derivatives
set sigi2 = sigi**2;

?  Do the later observations first, in order to calculate de0/dt
smpl smp2;
dysdt = ys(-1) + theta1*dysdt(-1);
dzsdt = zs(-1) + theta1*dzsdt(-1);
?mat zspys = zs'ys;
mat dzspysdt = zs'dysdt + dzsdt'ys;
set dljdt = (abs(theta1)<1)*[ -theta1/(1-theta1**2) +
                             np1*theta1**(2*np1-1)/(1-theta1**(2*np1))] ;
set dJacdt = -2*Jac*dljdt;  ? dljdt = (-1/2)*dlog(Jac)/dt
                            ?       = (-1/2)*(1/Jac)*dJacdt
?set e0 = zspys/Jac;
set de0dt = dzspysdt/Jac - dJacdt*e0/Jac;  ? e0/Jac = zspys/Jac**2
?e = ys - e0*zs;
?  Note:  we could leave out the  -de0dt*zs term in dedt for first
?  derivatives, as long as we do it in both smp2 and smp0.  These
?  terms add to zero, because e'zs = 0, due to orthogonality of zs
?  with the regression residuals e from regressing ys on zs.
?  However, the term is nonzero in the second derivatives (when
?  dedt is squared in each observation), so it cannot be left out
?  from usage there.
dedt = dysdt - e0*dzsdt - de0dt*zs;
?logl = log(sigi) - sigi2*e**2/2 - log(sqrt(2*pi));
dldt = -sigi2*e*dedt;
dlds = 1/sigi - sigi*e**2;

if ifhess; then; do;
d2ysdtt = 2*dysdt(-1) + theta1*d2ysdtt(-1);
d2zsdtt = 2*dzsdt(-1) + theta1*d2zsdtt(-1);
mat d2zydtt = d2zsdtt'ys + 2*(dzsdt'dysdt) + zs'd2ysdtt;
set d2ljdtt = (abs(theta1)<1)*[
        -1/(1-theta1**2) +
        -2*theta1**2/(1-theta1**2)**2 +
        (2*np1-1)*np1*theta1**(2*np1-2)/(1-theta1**(2*np1)) +
        2*np1**2*theta1**(4*np1-2)/(1-theta1**(2*np1))**2
        ] ;
set d2Jacdtt = -2*Jac*(d2ljdtt - 2*dljdt**2);
set d2e0dtt = d2zydtt/Jac - dJacdt*dzspysdt/Jac**2
      - de0dt*dJacdt/Jac + e0*(dJacdt/Jac)**2 - e0*d2Jacdtt/Jac ;
?d2edtt = d2ysdtt - e0*d2zsdtt;   ? too simple
d2edtt = d2ysdtt - e0*d2zsdtt - 2*de0dt*dzsdt - zs*d2e0dtt;  ? correct
d2ldtt = sigi2*(dedt**2 + e*d2edtt);
d2ldts = 2*sigi*e*dedt;
d2ldss = 1/sigi2 + e**2;
enddo;

?  Now do the first observation, since de0/dt has been calculated above.
smpl smp0;
?logl = -log(Jac)/2 - sigi2*e0**2/2;
dldt = dljdt - sigi2*e0*de0dt;
dlds = -sigi*e0**2;

if ifhess; then; do;
d2ldtt = -d2ljdtt + sigi2*(de0dt**2 + e0*d2e0dtt);
d2ldts = 2*sigi*e0*de0dt;
d2ldss = e0**2;
enddo;
smpl smpall;

endproc;