options double crt nodate;
name bjacls 'ARMA(1,1) Conditional L.S. on Box-Jenkins series A';
?
?   These 2 CLS benchmarks may be useful for comparison with
?   the exact ML benchmark.  (I prefer exact ML benchmarks in
?   general, because they avoid arbitrary choices for how the
?   presample u and e values are formed).
?   TSP code by Clint Cummins, 6/2000
?   CLS-1 benchmark coefficients agree to 7 digits with code
?     by Bruce McCullough, and to 5 digits with RATS and EViews
title 'ARMA(1,1)  with constant';


title 'CLS coded by hand';
set pi = 4*atan(1);
e = 0;
u = 0;
?  Benchmarks     CLS-1               CLS-2
?             TSP Proc ML         TSP Proc ML
?   mu        17.0937523895       17.0387257391
?  (s.e.)      (.106236540817)     (.098083165019)
?   phi         .906587036003       .903692939150
?              (.054937896394)     (.058144772525)
?   theta       .568809102811       .558037113139
?              (.118965622194)     (.121855335748)
?   LogL     -50.7889277349      -50.8698692896
?  Num.Obs.       196                 197
? 1. Standard errors are computed from numeric second derivatives
?    of the LogL.
? 2. Model is parameterized as:
?    (1 - phi*B)*(y - mu) = (1 - theta*B)*e    , e ~ N(0,sigma2)

title 'CLS-1:  drop one obs for u(-1), e(-1)=0';
param mu,17 phi,.9 theta,.6;
supres smpl;
options nwidth=20 signif=11;
ml(tol=1e-11,grad=c4,hiter=u) cls1 mu phi theta;
proc cls1;
  local u e sig2;
  smpl 1,197;
    u = a-mu;
    e = 0;
  smpl 2,197;
    genr(silent) e = u - phi*u(-1) + theta*e(-1);
    mat sig2 = (e'e)/@nob;
    set @logl = -(@nob/2)*(log(sig2) + 1 + log(2*pi) );
endproc;
smpl 1,197;

title 'CLS-2:  full sample, u(-1)=0, e(-1)=0';
ml(tol=1e-11,grad=c4,hiter=u) cls2 mu phi theta;
proc cls2;
  local u e sig2;
  smpl 1,197;
    u = a-mu;
    e = u;     ? e = u - phi*0     + theta*0    for first obs.
  smpl 2,197;
    genr(silent) e = u - phi*u(-1) + theta*e(-1);
  smpl 1,197;
    mat sig2 = (e'e)/@nob;
    set @logl = -(@nob/2)*(log(sig2) + 1 + log(2*pi) );
endproc;

title 'BJEST(NBACK=0) is not standard Conditional Least Squares';
?     BJEST with NBACK=0 does not exactly reproduce the above.
?  The reason is that BJEST still back-forecasts epsilon for
?  the first observation, instead of setting it to zero.
?  BJEST will always back-forecast NLAG epsilons, where NLAG=
?  max(NAR,NMA) for an ARMA model.  NBACK can used to do further
?  forecasts before these.
?     BJEST also parametrizes the intercept differently, with the
?  MA terms instead of using (y-mu):
?    (1 - phi*B)*y = alpha + (1 - theta*B)*e    , e ~ N(0,sigma2)
?   so  alpha = (1-phi)*mu  .
bjest(const,nar=1,nma=1,nback=0,nocumpl,noplot,nlag=0,tol=1e-7) a;
unmake @coef phi theta alpha;
smpl 1,5;
print a @res;
title 'verify calculation of second residual';
smpl 2,2;
res = a - phi*a(-1) - alpha + theta*@res(-1);
print res @res;

end;
noprint;
smpl 1,197;
read a;  ? p.525
17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 17.0
16.7 17.4 17.2 17.4 17.4 17.0 17.3 17.2 17.4 16.8
17.1 17.4 17.4 17.5 17.4 17.6 17.4 17.3 17.0 17.8
17.5 18.1 17.5 17.4 17.4 17.1 17.6 17.7 17.4 17.8
17.6 17.5 16.5 17.8 17.3 17.3 17.1 17.4 16.9 17.3
17.6 16.9 16.7 16.8 16.8 17.2 16.8 17.6 17.2 16.6
17.1 16.9 16.6 18.0 17.2 17.3 17.0 16.9 17.3 16.8
17.3 17.4 17.7 16.8 16.9 17.0 16.9 17.0 16.6 16.7
16.8 16.7 16.4 16.5 16.4 16.6 16.5 16.7 16.4 16.4
16.2 16.4 16.3 16.4 17.0 16.9 17.1 17.1 16.7 16.9
16.5 17.2 16.4 17.0 17.0 16.7 16.2 16.6 16.9 16.5
16.6 16.6 17.0 17.1 17.1 16.7 16.8 16.3 16.6 16.8
16.9 17.1 16.8 17.0 17.2 17.3 17.2 17.3 17.2 17.2
17.5 16.9 16.9 16.9 17.0 16.5 16.7 16.8 16.7 16.7
16.6 16.5 17.0 16.7 16.7 16.9 17.4 17.1 17.0 16.8
17.2 17.2 17.4 17.2 16.9 16.8 17.0 17.4 17.2 17.2
17.1 17.1 17.1 17.4 17.2 16.9 16.9 17.0 16.7 16.9
17.3 17.8 17.8 17.6 17.5 17.0 16.9 17.1 17.2 17.4
17.5 17.9 17.0 17.0 17.0 17.2 17.3 17.4 17.4 17.0
18.0 18.2 17.6 17.8 17.7 17.2 17.4 ;