options crt;  ? 7e6
?
?  The AR1 command was changed in TSP 4.4 (June 1998) so that it
?  combines grid search and iterations for the GLS objective function.
?  This means it will automatically note the local optima in sections
?  (d)-(f).  Most of the commands here are intended for earlier versions
?  of TSP, which do not automatically find multiple local optima.
?  So that this exercise will run in earlier versions of
?  TSP or later versions, enter your version number below.  (Use 4.5
?  as the version number if your copy of TSP is dated 6/1998 or later).
set version = 4.5;
in nerc;
freq a;
smpl 51,84;
lnkwh = log(kwh);
lnpe = log(pelec);
lngnp = log(gnp);
smpl 52,84;
regopt(pvprint,stars) dhalt;  ? This P-value is on by default in TSP 4.4
regopt(dwpval=exact) ;  ? turn on exact Durbin-Watson P-value
?
?  Note that the P-value for Durbin-Watson is lower than that for
?  Durbin's h or Durbin's m.  So the D-W is still fairly powerful
?  (and conservative) with lagged dependent variables, even though
?  the DW is biased towards 2 and the P-value is biased away from zero.
ols lnkwh c lnkwh(-1) lnpe lngnp;
regopt(dwpval=approx) ;
copy @dhalt dha;  ? Save dhalt value so it is not overwritten by later OLS
copy %dhalt %dha;
e = @res;  ? save residuals for (c)
frml ephi phi = 1-lnkwh(-1);
frml eb1 b1 = lnpe/(1-lnkwh(-1));
frml eb2 b2 = lngnp/(1-lnkwh(-1));
analyz ephi eb1 eb2;  ? phi and long-run elasticities  (a)
?
? do D-W test even though it's biased towards 2 due to lagged dependent var.
? However, it often has better power than Durbin's h and h-alternative.
?
? early versions of TSP 4.4 will refuse to calculate the
? D-W when there is a lagged dependent variable,
? but it can be computed from the residuals.
?
if (version=4.4); then;
olsq(silent) e c;
?
if (version >= 4.4); then;
  set dwtest = (%dw < .05);
else;
  set dwtest = (@dw < 1.26);
if (dwtest); then;  title 'Reject null: rho=0';
             else;  title 'Accept null: rho=0';
set rho = 1 - @dw/2;
print rho;  ? (b)
?
title "Compute Durbin's m manually";
smpl 53,84;
ols e c lnpe lngnp lnkwh(-1) e(-1);
set mstat = @t(5);
cdf(inv,t,df=27) .05 tcrit;
print dha %dha mstat tcrit;  ? (c) - compare manual and automatic Durbin's m
?
smpl 52,84;  ? AR1 will adjust the SMPL automatically
if (version >= 4.5); then;
 ?
 ? The 2 local optima are:  rho=.7001147 (local) and rho=.9431796 (global).
 ar1 lnkwh c lnkwh(-1) lnpe lngnp;  ? combines grid search and iteration
 ?
 ? if you want to repeat the pure grid search, use maxit=0:
 ?ar1(maxit=0,rmin=0,rmax=1,rstep=.05) lnkwh c lnkwh(-1) lnpe lngnp;
 ?ar1(maxit=0,rmin=.65,rmax=.75,rstep=.01) lnkwh c lnkwh(-1) lnpe lngnp;
 ?ar1(maxit=0,rmin=.9,rmax=1,rstep=.01) lnkwh c lnkwh(-1) lnpe lngnp;
else; do;
 ar1(method=hilu,rmin=0,rmax=1,rstep=.05) lnkwh c lnkwh(-1) lnpe lngnp;
 ar1(method=hilu,rmin=.65,rmax=.75,rstep=.01) lnkwh c lnkwh(-1) lnpe lngnp;
 ar1(method=hilu,rmin=.9,rmax=1,rstep=.01) lnkwh c lnkwh(-1) lnpe lngnp;
 ?
 ? The last equation yields an (apparent) global optimum at rho=.94
 ? Note that there is no guarantee that this is even within .01 of the
 ? true global optimum.  There might be a true global optimum at, say
 ? rho = -.5 or even rho = .02 that was not picked up by the
 ? initial grid search.  All we can say for sure is that .94 is the
 ? optimum of the rho values we tried.  But the log likelihood appears
 ? to be fairly smooth as a function of rho, so it is quite likely that
 ? the global optimum is .94 +/- .01 .
enddo;
?
analyz ephi eb1 eb2;  ? phi and long-run elasticities  (d)
?
if (version >= 4.5); then; do;
 title 'one-step Cochrane-Orcutt estimator';
 ar1(method=corc,maxit=1,rstart=0) lnkwh c lnkwh(-1) lnpe lngnp;  ? (e)
 title 'global GLS estimator estimator (same as (d))';
 ar1 lnkwh c lnkwh(-1) lnpe lngnp;  ? combines grid search and iteration
enddo;
else; do;
 title 'look at first iteration for one-step Cochrane-Orcutt estimator';
 ar1(method=corc,print) lnkwh c lnkwh(-1) lnpe lngnp;  ? (e)
enddo;
?
? AR1 will normally refuse to do ML because the lagged dependent variable
?  will be correlated with the residual in the first (non rho-differenced)
?  observation.  This correlation does not have much effect if there is
?  a large number of observations (i.e. it is asymptotically ignorable).
?  We can defeat this check by renaming the lagged dependent  variable.
?  It would then be necessary to use Hatanaka's procedure to obtain
?  correct standard errors (i.e. including RHO).
?  The ML estimate of RHO is .658564 (accurate to 6 digits).
?
lnkwhl = lnkwh(-1);
ar1(maxit=100) lnkwh c lnkwhl lnpe lngnp;  ? METHOD=ML is the default  (f)
?
?  In TSP 4.4 and earlier, don't do ANALYZ here because the standard
?  errors are wrong (although the point estimates of the elasticities
?  would be OK)
if (version >= 4.5); then; do;
  frml ephi phi = 1-lnkwhL;
  frml eb1 b1 = lnpe/(1-lnkwhL);
  frml eb2 b2 = lngnp/(1-lnkwhL);
  analyz ephi eb1 eb2;
enddo;
?
? Additional estimates to check for local optima.  (TSP 4.4 and earlier)
? Note:  of the 3 datasets I've seen with local optima for CORC/HILU,
?  none of them appears to have local optima for ML.
?
if (version <= 4.4); then; do;
  ar1(method=mlgrid) lnkwh c lnkwhl lnpe lngnp;  ? this should normally
                              ? be sufficient to check for local optima
  ar1(rstart=.5) lnkwh c lnkwhl lnpe lngnp;  ? here are more checks
  ar1(rstart=.7) lnkwh c lnkwhl lnpe lngnp;  ? as suggested in the text
  ar1(rstart=.9) lnkwh c lnkwhl lnpe lngnp;
  ar1(rstart=.95) lnkwh c lnkwhl lnpe lngnp;
enddo;
?
? The global optimum depends not only on the first residual, but on
? the Jacobian term log(1-rho**2) that helps bound rho away from +/- 1.
? To assess the contributions of all these terms, plot the logLikelihood
? for HILU, "PRWI" (i.e. Prais-Winston GLS procedure with the first
? observation retained but without the Jacobian term), and ML as a
? function of rho, say in the range rho= [.01,.99] .
? This could be a part (h) of the exercise.
? Let's assume we're using a graphics version of TSP.
?    Judging from the graph, the ML and Prais-Winsten objective functions
? are very close to each other and unimodal, while the HILU objective
? function is bimodal.  So the inclusion of the first observation is more
? important than the Jacobian term for unimodality in this particular
? example.
?
mform(nrow=99,ncol=1) rhov;
copy rhov hiluv;
copy rhov prwiv;
copy rhov mlv;
supres @coef @logl @nob;
noplot;
do i=1,99;
?do i=10,10;
  set rho = i/100;
  set rhov(i) = rho;
?  ar1(method=corc,rstart=rho,maxit=0,silent) lnkwh c lnkwhl lnpe lngnp;
  if (version >= 4.5); then;
    ar1(method=hilu,rstart=rho,silent,maxit=0) lnkwh c lnkwhl lnpe lngnp;
  else;
    ar1(method=hilu,rmin=rho,rmax=rho,silent) lnkwh c lnkwhl lnpe lngnp;
  set hiluv(i) = @logl;
?  ar1(method=ml,rstart=rho,maxit=0,silent) lnkwh c lnkwhl lnpe lngnp;
  if (version >= 4.5); then;
    ar1(method=mlgrid,rstart=rho,maxit=0,silent) lnkwh c lnkwhl lnpe lngnp;
  else;
    ar1(method=mlgrid,rmin=rho,rmax=rho,silent) lnkwh c lnkwhl lnpe lngnp;
  set mlv(i) = @logl;
  set prwiv(i) = @logl - log(1-rho**2)/2;
enddo;
freq n;
smpl 1,99;
dot rho hilu prwi ml;
  unmake .v .;
enddot;
options display=vga;
graph(line,preview,title='Logl(rho)') rho hilu prwi ml;
write(file='7e6.dat') rho hilu prwi ml;
?smpl 10,10;
?print rho hilu gls ml;