options crt;  ? 10e10
? Note:  the grid search parts of this exercise take a lot of time to run
? (On a Pentium/166Mhz, 25 seconds for the first grid of 49 points, NAR=0).
in taylor;
freq q;
smpl 53:1 75:4;
ydev = log(gnp/potgnp);
rm1 = log(m1/gnpdef);
smpl 53:1 75:3;
inf = log(gnpdef(+1)/gnpdef);

smpl 53:3 75:3;
frml eq1 ydev =     (b0 + b1*ydev(-1) + b2*ydev(-2) + b3*rm1 + b4*rm1(-1)
                    + b5*inf(-1) + b6*time + b5*g0)/(1-b5*g1);
frml eq2 inf  = (g1*(b0 + b1*ydev(-1) + b2*ydev(-2) + b3*rm1 + b4*rm1(-1)
                                 + b6*time)
                 + inf(-1) + g0)/(1-b5*g1);
param b0-b6 g0 g1;
lsq eq1 eq2;  ? restricted model
copy @logl loglr;
var ydev inf | c ydev(-1) ydev(-2) rm1 rm1(-1) inf(-1) time;  ? unrestricted
set LR = 2*(@logl - loglr);
cdf(chisq,df=5) LR;  ? (a)

smpl 53:2 75:3;
e1 = 0;  ? for (10.54) presample zero value
e2 = 0;  ? e1 = eta, e2 = epsilon
smpl 53:3 75:3;
ydev1 = ydev(-1);
ydev2 = ydev(-2);
rm11 = rm1(-1);
inf1 = inf(-1);
list rhs c ydev1 ydev2 rm1 rm11 inf1 time;
frml eq1s ydevs1 = (b0*cs1 + b1*ydev1s1 + b2*ydev2s1 + b3*rm1s1 + b4*rm11s1
                    + b5*inf1s1 + b6*times1 + b5*g0*cs1)/(1-b5*g1);
frml eq2s infs2  = (g1*(b0*cs2 + b1*ydev1s2 + b2*ydev2s2 + b3*rm1s2 + b4*rm11s2
                                 + b6*times2)
                 + inf1s2 + g0*cs2)/(1-b5*g1);
dot ydev inf rhs;
  .s1 = .;  ? ydevs1 = ydev,  ydev1s1 = ydev1,  ...   times1 = time
  .s2 = .;  ? for first observation of *1 (s1) and *2 (s2) transformations
enddot;

? initialize central points in grid, +/- width, and best function value
const phi1best 0 phi2best 0 dphi .99 depebest 1e30;
supres smpl coef logl nob;  ? turn off output
regopt(nopvcalc) dw;  ? speed up LSQ, by skipping DW P-value
?  Note: we can compute from Taylor's results that his optimal value
?  of phi1 was about -.67 and the optimal value of phi2 was about .66 :
const b5tay, -.447 g1tay,.018 thet1tay,-.38 thet2tay,.67;
set phi1tay = (b5tay*thet2tay+thet1tay)/(1-g1tay*b5tay);
set phi2tay = (g1tay*thet1tay+thet2tay)/(1-g1tay*b5tay);
print phi1tay phi2tay;
?  Since those value are interior to the grid (-1,1) that we use
?  here, we should find that optimum.  We can't grid search directly
?  on theta1,theta2, because we would have to know the values of
?  b5,g1 to compute phi1,phi2 and do the data transformations.
?  Using the grid search below, we find phi1=.0543 and phi2=.3413 ,
?  +/- .00015, with depebest = 6.89305d-6 .

? Do Grid Search on phi1,phi2  (b)
const narmax 8;
date time1;
do nar=0,narmax;
  if nar=narmax; then; do;   ? set up printing for final result
    set dphi=0;
    nosupres coef logl nob;  ? turn on output
  enddo;
  ? narrow the grid, around previous best answers
  set phi1L = phi1best - dphi;
  set phi1H = phi1best + dphi;
  set phi2L = phi2best - dphi;
  set phi2H = phi2best + dphi;
  set dphi = dphi/3;
  ? Grid search over phi1,phi2, given their bounds
  ? phi1L <= phi1 <= phi1H
  ? phi2L <= phi2 <= phi2H
  ? and number of values nphi
  do phi1=phi1L,phi1H,dphi;  ? grid on phis, not thetas
    smpl 53:4 75:3;  ? observations 2 to T
    dot ydev rhs;
      .s1 = . + phi1*.s1(-1);  ? recursive transformations
    enddot;
    do phi2=phi2L,phi2H,dphi;
      smpl 53:4 75:3;  ? observations 2 to T
      dot inf rhs;
        .s2 = . + phi2*.s2(-1);
      enddot;
      smpl 53:3 75:3;   ? all observations for estimation
      if nar=0; then; do;
        param b0,0 b1,0 b2,0 b3,0 b4,0 b5,0 b6,0 g0,0 g1,0;
        lsq(silent) eq1s eq2s;
      enddo;  else;
        lsq(silent,wname=own) eq1s eq2s;
      ? compute depe = det(e'e), equivalent to e'Omega"e
      unmake @res w1 w2;
      e2 = w2 + phi2*e2(-1);
      e1 = w1 + phi1*e2(-1);  ? (10.56)
      mmake e e1 e2;
      mat depe = det(e'e);
      ?print depe phi1 phi2;
      if (depe < depebest); then; do;
        set depebest = depe;
        set phi1best = phi1;
        set phi2best = phi2;
        set theta1 = phi1 - b5*phi2;
        set theta2 = phi2 - g1*phi1;
      enddo;
    enddo;
  enddo;
  date time2;  set seconds = time2-time1;  set time1 = time2;
  print nar theta1 theta2 depebest phi1best phi2best dphi seconds;
enddo;
set loglr = @logl;

? Do the same type of grid search on the unrestricted model
form(param,varpref=y_) eq1u ydevs1
      cs1 ydev1s1 ydev2s1 rm1s1 rm11s1 inf1s1 times1;
form(param,varpref=i_) eq2u infs2
      cs2 ydev1s2 ydev2s2 rm1s2 rm11s2 inf1s2 times2;

? initialize central points in grid, +/- width, and best function value
const phi1best 0 phi2best 0 dphi .99 depebest 1e30;
supres smpl coef logl nob;  ? turn off output

? Do Grid Search on phi1,phi2  (c)
const narmax 8;
date time1;
do nar=0,narmax;
  if nar=narmax; then; do;   ? set up printing for final result
    set dphi=0;
    nosupres coef logl nob;  ? turn on output
  enddo;
  ? narrow the grid, around previous best answers
  set phi1L = phi1best - dphi;
  set phi1H = phi1best + dphi;
  set phi2L = phi2best - dphi;
  set phi2H = phi2best + dphi;
  set dphi = dphi/3;
  ? Grid search over phi1,phi2, given their bounds
  ? phi1L <= phi1 <= phi1H
  ? phi2L <= phi2 <= phi2H
  ? and number of values nphi
  do phi1=phi1L,phi1H,dphi;  ? grid on phis, not thetas
    smpl 53:4 75:3;  ? observations 2 to T
    dot ydev rhs;
      .s1 = . + phi1*.s1(-1);  ? recursive transformations
    enddot;
    do phi2=phi2L,phi2H,dphi;
      smpl 53:4 75:3;  ? observations 2 to T
      dot inf rhs;
        .s2 = . + phi2*.s2(-1);
      enddot;
      smpl 53:3 75:3;   ? all observations for estimation
      ? Since the unrestricted model is linear, it will be solved
      ? in two iterations regardless of starting values, so there is
      ? no need to use wname=own.
      lsq(silent) eq1u eq2u;
      ? compute depe = det(e'e), equivalent to e'Omega"e
      unmake @res w1 w2;
      e2 = w2 + phi2*e2(-1);
      e1 = w1 + phi1*e2(-1);  ? (10.56)
      mmake e e1 e2;
      mat depe = det(e'e);
      ?print depe phi1 phi2;
      if (depe < depebest); then; do;
        set depebest = depe;
        set phi1best = phi1;
        set phi2best = phi2;
      enddo;
    enddo;
  enddo;
  date time2;  set seconds = time2-time1;  set time1 = time2;
  print nar depebest phi1best phi2best dphi seconds;
enddo;
set LR = 2*(@logl - loglr);
cdf(chisq,df=5) LR;  ? (c)