options crt;
?  name armax 'Regression with MA(1) error term' ;
?
?  ARMAX(1,1) estimation using the method in Andrew Harvey, 
?  The Econometric Analysis of Time Series (1990), p.273.
?  TSP 4.3A (May 1995)
?  Author: Clint Cummins
? 
set n=2000;
?
?  Simulated data generation.
?
param y0,1 theta,.5 alpha,.4 beta,2;
smpl 1,n;
e = 0; y = y0;
deda = 0; dedb = 0; dedt = 0;
smpl 2,n;
random(seedin=1498742) e x;
u = e + theta*e(-1);
y = alpha*y(-1) + beta*x + u;
?
?  Estimation
?  1. initial values from OLS residuals and BJEST.
?
ols y y(-1) x;
unmake @coef alpha beta;
u = @res;
bjest(nma=1,nocumplo,noplot) u;
set theta = -@coef(1);       ? sign of coef is reversed in BJEST
list bnames alpha beta theta;
mmake b bnames;
?
?  2. iteration
?
?  Equation for residual:
?
frml eqe e = y - alpha*y(-1) - beta*x - theta*e(-1);
?
?  Derivatives.
?
frml d1 deda = -y(-1) - theta*deda(-1);
frml d2 dedb = -x     - theta*dedb(-1);
frml d3 dedt = -e(-1) - theta*dedt(-1);
?
genr eqe;
mat ssrold = e'e;
?
mat novar = nrow(b);
supres @logl @nob @coef;
const maxit,20 maxsqz,10 tol,.001;
do iter=1,maxit;
  genr(silent) d1; genr(silent) d2; genr(silent) d3;
  ? Gauss-Newton iteration:  regress e on derivatives
  ols(silent) e deda dedb dedt;
  rename @coef delt;
  ? test for convergence
  mat worst = max(abs(delt)/(abs(b) + 10*tol));
  if worst<tol; then; goto 20;
  do isqz=1,maxsqz;
    mat testb = b - delt/(4**(isqz-1));
    unmake testb bnames;
    genr(silent) eqe;
    mat ssrnew = e'e;
    ?write ssrold ssrnew isqz;
    if ssrnew<ssrold; then; goto 15;
  enddo;
  write(form='(" Failure to improve after ",F4.0," squeezes.")') isqz;
  write b delt;
  goto 18;
  15 write iter,ssrnew;
  set ssrold = ssrnew;
  copy testb b;
enddo;

18 write(form='(" Failure to converge after ",F4.0," iterations.")') iter;
goto 30;

20 write(form='(" Convergence after ",F4.0," iterations.")') iter;

?  Print results of estimation.
?
30 unmake b bnames;
write ssrold bnames;
nosupres @coef;
tstats(names=bnames) b @vcov;