options crt;  ? 6e9 -- a difficult exercise
in kopcke;
freq q;
frml e637a sigma = sumphic/(1-sumom);  ? (6.37), but note change in sign
frml e637b eta = sumphiy/(1-sumom);    ? for sumom because we est. -om* .
const phic1-phic12; const phiy1-phiy12;  ? make sure all phis are defined
DOT 1-11;   ? equations for part (d) -- PDL w/ AR(2) error
  FRML Ephic. phic. = C1*A.1 + C2*A.2 + C3*A.3 + C4*A.4;
  ? FRML Ephic1 phic1 = C1*A11 + C2*A12 + C3*A13 + C4*A14;
  ? ... ETC. ...
  ? FRML Ephic11 phic11 = C1*A111 + C2*A112 + C3*A113 + C4*A114;
  FRML Ephiy. phiy. = D1*A.1 + D2*A.2 + D3*A.3 + D4*A.4;
  CONST A.1 A.2 A.3 A.4;  ? fill in values via AMAT in eqsetup below
ENDDOT;
eqsetup;  ? see PROC eqsetup below -- sets up equations for ANALYZ
? equations to copy om* (so all params are together in @VCOVA for ANALYZ)
DOT 1-3;  FRML Eom. om. = om. ;  ENDDOT;
LIST AEQS Ephic1-Ephic11 Ephiy1-Ephiy11 Eom1;
PARAM A0 C1-C4 D1-D4;   ? NOTE:  C1-C4, D1-D4 ARE THE "SCRAMBLED" COEFFICIENTS
?
set equip=2;  ? T/F flag for Equipment vs. Structures
dot e s;
set equip = equip-1;   ? will be 1, then 0
smpl 52:1,86:3;
lk. = log(k.lag(+1));  ? note lead on lagged K to make K.LAG current --
                       ? this means we can only estimate up to 86:3, vs. 86:4
ly = log(y);
lpc. = -log(c.);
smpl 52:2,86:3;  ? first difference variables for part (b)
dlk. = lk. - lk.(-1);
dly  = ly - ly(-1);
dlpc. = lpc. - lpc.(-1);
smpl 56:1,86:3;
?
?  set up summary coefficient and confidence interval files
write(file='coefa.6e9',format=
"('   m   s    sigma     eta [ siglo , sighi][ etalo , etahi]   SBIC')");
write(file='coefb.6e9',format=
"('   m   s    sigma     eta [ siglo , sighi][ etalo , etahi]   SBIC')");
write(file='coefc.6e9',format=
"('   m   s    sigma     eta [ siglo , sighi][ etalo , etahi]  rho')");
?
?  PDLs for different lag lengths
do m=9,12;
 do s=1,3;
  set sneg = -s;
  ? If output is not suppressed with SILENT, 6e9.out will be 500K bytes.
  olsq(silent) lk. c lpc.(4,m,none) ly(4,m,none)
             lk.(-1)-lk.(sneg);  ? note special "variable lag list"
  do637 1;  ? (a) -- This is a PROC listed below which uses ANALYZ to
            ? compute standard errors for sigma and eta
  olsq(silent) dlk. dlpc.(4,m,none) dly(4,m,none) dlk.(-1)-dlk.(sneg);
  do637 2;  ? (b)
  ? to constrain rho >=0 , use METHOD=HILU, RMIN=0
  ar1(silent,method=hilu,rmin=0,rstep=.05)
   lk. c lpc.(4,m,none) ly(4,m,none) lk.(-1)-lk.(sneg);
  do637 3;  ? (c)
 enddo;
enddo;
?
? (d) -- doing a higher-order AR with PDLs is feasible, but rather messy,
? because a very complex nonlinear equation must be constructed for LSQ.
? Let's list the complications:
? 1. The PDLs have to be written explicitly, including the associated
?    Lagrangian interpolation polynomials.
? 2. Writing the AR pseudo-differencing will involve duplicating the terms
?    in (1) and lagging all the variables.  This is actually fairly simple
?    with the "lagged eqsub" feature as of 6/1996.
? 3. The sum of phics and phiys (for ANALYZ) will now involve a
?    wild bunch of the explicit terms in (1).
? The program above must be run first to find the "preferred specification".
? At least the do loops don't have to be handled anymore, and the test for
? AR(p) vs. AR(1) is straightforward.
?
?  There is an independent PROC on the TSP Examples web page called PDLAR
?  which will do a single PDL, with or without restrictions, for many
?  different AR orders.  We can't use it here, because there are 2 PDLs.
?
? Preferred models, based on min SBIC for the OLS models (less easy to
? compute SBIC for AR1 -- see 6e7 and 6e8, and rho is close to one anyway):
? Equipment:   m=11, s=3 (SBIC= -12.7663).
? Structures:  m=11, s=1 (SBIC= -14.3586).
?
PARAM om1-om3;
if equip; then;
  set s=3;
else;
  const s 1 om2,0 om3,0;  ? zero extra lagged lks for Structures
set sneg=-s;
set m=11;
smpl 55:4,86:3;  ? adjust SMPL because CORC drops an obs.
?  Tighten the convergence tolerance, to make sure that AR1 and LSQ
?  agree very closely on the estimates for the AR(1) case, before
?  extending to the AR(2) model, where no such check is available.
ar1(rstart=.7,method=corc,TOL=.001)
  lk. c lpc.(4,m,none) ly(4,m,none) lk.(-1)-lk.(sneg);
smpl 56:1,86:3;
?
? Setup AR(2) estimation with PDLs in LSQ (nonlinear least squares)
?
? AR(2) error:  Y = X*B + U;  U = E + P1*U(-1) + P2*U(-2);
?       so      Y = X*B + P1*[Y(-1)-X(-1)*B] + P2*[Y(-2)-X(-2)*B] + E;
?
?  First write the equation for the residual U
?
FRML U. lk. - (A0 +
               phic1*lpc.      + phic2*lpc.(-1)  + phic3*lpc.(-2) +
               phic4*lpc.(-3)  + phic5*lpc.(-4)  + phic6*lpc.(-5) +
               phic7*lpc.(-6)  + phic8*lpc.(-7)  + phic9*lpc.(-8) +
               phic10*lpc.(-9) + phic11*lpc.(-10) +
               phiy1*ly      + phiy2*ly(-1)  + phiy3*ly(-2) +
               phiy4*ly(-3)  + phiy5*ly(-4)  + phiy6*ly(-5) +
               phiy7*ly(-6)  + phiy8*ly(-7)  + phiy9*ly(-8) +
               phiy10*ly(-9) + phiy11*ly(-10) +
               om1*lk.(-1) + om2*lk.(-2) + om3*lk.(-3) );
?
?  Now the AR2 equation for the pseudo differenced residual
?  E = U - P1*U(-1) - P2*U(-2)
?
FRML AR2. U. - P1*U.(-1) - P2*U.(-2);
?
EQSUB AR2. U. Ephic1-Ephic11 Ephiy1-Ephiy11;   ? substitute in scrambled coefs
CONST P1 0 P2 0;  ? start with plain OLS model to generate good starting
                  ? values for everything but RHO (P1)
LSQ(silent) AR2.;
? AR(1)
CONST P2 0;
PARAM P1 .9;  ? ESTIMATE FIRST AR COEFFICIENT ONLY
LSQ(silent,TOL=.001) AR2.;  copy @logl loglar1;
ar2o;
? AR(2)
PARAM P1 .6 P2 .3;  ? ESTIMATE BOTH AR COEFFICIENTS
LSQ(silent,TOL=.001,maxit=50) AR2.;
ar2o;
title 'Tests for AR(1) vs. AR(2)';
frml wald p2wald = p2;  ANALYZ WALD;
set lr = 2*(@logl-loglar1);  cdf(chisq,df=1) lr;
?
enddot;
?
Proc ar2o;  ? output for AR(2) just above
  if equip; then;
    analyz AEQS Eom2 Eom3;  ? consolidate output for phi*, om1-om3
  else;
    analyz AEQS;  ? Structures, s=1:  phi*, om1 only
  analyz(noprint,vcov=@VCOVA,names=@RNMSA,COEF=@COEFA) e637a e637b;
  lohi sl sh el eh;
  mmake v m s sigma eta sl sh el eh p1 p2;
  write(file='coefc.6e9',format='(f5.0,f4.0,6f8.4,2f5.2)') v;
endproc;
?
?----------------------------------------------
?
?  Here is the PROC which sets up the equations
?  (it could have been left at the top of the program, but
?  it's cleaner to hide it down here)
?
Proc eqsetup;
?
?  equations to set up for ANALYZ in PROC  do637  below
?
?  note 0/1 expressions like (m>9) to "turn on" terms in the sum based on m
frml sumphic phic1+phic2+phic3+phic4+phic5+phic6+phic7+phic8+phic9+
             (m>9)*phic10 + (m>10)*phic11 + (m>11)*phic12;
frml sumphiy phiy1+phiy2+phiy3+phiy4+phiy5+phiy6+phiy7+phiy8+phiy9+
             (m>9)*phiy10 + (m>10)*phiy11 + (m>11)*phiy12;
frml sumom om1 + (s>1)*om2 + (s>2)*om3;
eqsub e637a sumphic sumom;
eqsub e637b sumphiy sumom;
?
?  set up the PDL scrambling matrices for (d)
?
AMAT 4,11,A;
?PRINT A;
DOT(index=I) 1-11;
  DOT(index=J) 1-4;
    SET A.. = A(I,J);  ? SETS A02 = A(1,2), ETC. for FRML EPHIC1, etc.
  ENDDOT;
ENDDOT;
endproc;
?
?  Here is the PROC which handles the messy lists for ANALYZ
?
Proc do637 case;
? the argument case handles part (a), (b), or (c)
? set up the phi lists
?
  LIST(LAST=m,prefix=phic) phics;  ? same as  List phics phic1-phic9; (if m=9)
  LIST(LAST=m,prefix=phiy) phiys;
? set up list for omegas based on s
  LIST(LAST=s,prefix=om) oms;
  const om1-om3;  ? make sure all oms are defined
? list of coefficient names
  if case=1; then; list bs a phics phiys oms;      ? (a) rho=0
  if case=2; then; list bs   phics phiys oms;      ? (b) rho=1 w/ no C
  if case=3; then; list bs a phics phiys oms rho;  ? (c) AR1 w/ lagged dep.
  param bs;
  unmake @coef bs;
  analyz(noprint,names=bs) e637a e637b;
  lohi sl sh el eh;
  ? put values in vector for write
  if case<3; then; mmake v m s sigma eta sl sh el eh @sbic;
             else; mmake v m s sigma eta sl sh el eh rho;
  if case=1; then; write(file='coefa.6e9',format='(f5.0,f4.0,6f8.4,f9.4)') v;
  if case=2; then; write(file='coefb.6e9',format='(f5.0,f4.0,6f8.4,f9.4)') v;
  if case=3; then; write(file='coefc.6e9',
                         format='(f5.0,f4.0,6f8.4,f6.2)') v;
endproc;
?
Proc lohi sl sh el eh;
  set sl=sigma-1.96*@sesa(1); set sh=sigma+1.96*@sesa(1);
  set el=  eta-1.96*@sesa(2); set eh=  eta+1.96*@sesa(2);
endproc;
?
PROC AMAT N,P,A;   ? for PDL with LSQ (d)
?  THIS MAKES THE LAGRANGIAN INTERPOLATION POLYNOMIALS
?  USING COOPER'S FORMULA.  BASED ON TSP SUBROUTINE AMAT.
?  DOES NOT ALLOW ZERO RESTRICTIONS IN THIS VERSION, FOR SIMPLICITY.
?  There are more general versions of this PROC in the Examples on the
?  TSP web page.
LOCAL INC,J,DENOM,I,T,PHINUM;
MFORM(NROW=P,NCOL=N) A=0;
IF (N.GT.1); THEN; SET INC = (P-1)/(N-1);
DO J=1,N;
  SET DENOM=1;
  DO I=1,N;
    IF (I.NE.J); THEN; SET DENOM = DENOM*INC*(J-I);
  ENDDO;
  DO T=1,P;
    SET PHINUM=1;
    DO I=1,N;
      IF (I.NE.J); THEN; SET PHINUM = PHINUM*(T - 1 - INC*(I-1));
    ENDDO;
    SET A(T,J) = PHINUM/DENOM;
  ENDDO;
ENDDO;
ENDPROC;