options crt nodate;
name ar1lag 'AR(1) with lagged dependent variable';
?  Demonstrates:
?  1. Testing for serial correlation with lagged dependent variable
?     (Durbin's h, Durbin's m, and Durbin-Watson).
?  2. Multiple local optima for AR(1) model (conditional ML) with
?     lagged dependent variable.
?        The exact ML estimator does not usually have this problem, but
?     it is often not used, because it is biased in small samples due
?     to the correlation between the lagged dependent variable and
?     the residual in the transformed initial observation.
?  3. Calculation of standard errors of coefficients for AR(1) with
?     lagged dependent variable (information matrix is no longer block
?     diagonal in rho).
?  Follows Berndt, E.R. "The Practice of Econometrics," (1990),
?  Chapter 7, Exercise 6.
load;
lnkwh = log(kwh);
lnpe = log(pelec);
lngnp = log(gnp);
smpl 52,84;

title 'OLS - testing for serial correlation with lagged dependent var.';
?    Type of Test      TestStat P-value
?      Durbin-Watson = 1.69145 [.074]   (exact P-value)
?         Durbin's h = .454320 [.650]   computed using rho as regr. coef.
?                                       of e on e(-1), not rho = 1-DW/2.
?                                       (both ways are valid)
?    Durbin's h alt. = .740731 [.459]   a.k.a. "Durbin's m"
?    Log Likelihood  = 83.2166
? Note that even though the Durbin-Watson is biased towards 2,
? it often dominates the other tests in power.  I.e. its P-value is
? biased upwards, but is still the lowest of the 3!
regopt(dwpval=exact);
ols lnkwh c lnkwh(-1) lnpe lngnp;
regopt(dwpval=bounds);

title 'AR(1) (conditional ML) with default starting values';
?              TSP      Berndt
? Const.   -1.97242
? (s.e.)    (.351050)          Note:  consistent standard errors.
? LNKWH(-1)  .433996
?           (.070534)
? LNPE      -.142029
?           (.050908)
? LNGNP      .883004
?           (.105590)
? rho        .700115     near .70  (.67 is premature convergence)
?           (.106674)
? LogL       87.7452          (apparent local optimum)
ar1(method=corc,tol=1e-7,maxit=100) lnkwh c lnkwh(-1) lnpe lngnp;

title 'Repeat conditional ML AR(1) via nonlinear least squares';
? (to verify that the standard errors are computed correctly).
form(param,coefpref=b) ar1;
smpl 53,84;
lsq(tol=1e-7) ar1;

title 'AR(1) (conditional ML) with special starting value for rho';
?               TSP      Berndt
? Const.     1.80567
? (s.e.)    (1.78843)
? LNKWH(-1)   .137712
?            (.111433)
? LNPE       -.178422
?            (.090199)
? LNGNP       .704903
?            (.157024)
? rho         .943180     .94
?            (.021153)
? LogL       88.6444          (apparent global optimum)
param rhoar1,.95;   ? special starting value to locate better optimum
lsq(tol=1e-7) ar1;

? Can verify that AR(1) with grid search locates best optimum
?  (to fineness of grid used, .05 here)
smpl 52,84;
ar1(method=hilu,rmin=0,rmax=1,rstep=.05) lnkwh c lnkwh(-1) lnpe lngnp;
? Can also show that rho=.66 (mentioned in Berndt book as .67)
? results from tol=.01 (old looser convergence tolerance)
?ar1(method=corc,tol=.01,print) lnkwh c lnkwh(-1) lnpe lngnp;
end;
noprint;
freq a;
smpl 51,84;  ? data from Berndt (1990) Chapter 7 (NERC dataset)
load
  KWH   PELEC   GNP;
 330.00 3.12  579.40
 356.00 3.09  600.80
 396.00 3.01  623.60
 424.00 2.97  616.10
 497.00 2.75  657.50
 546.00 2.61  671.60
 576.00 2.57  683.80
 588.00 2.59  680.90
 647.00 2.50  721.70
 689.00 2.65  737.20
 723.00 2.63  756.60
 778.00 2.55  800.30
 833.00 2.47  832.50
 896.00 2.38  876.40
 954.00 2.29  929.30
1035.00 2.16  984.80
1099.00 2.09 1011.40
1203.00 1.97 1058.10
1314.00 1.88 1087.60
1392.00 1.83 1085.60
1470.00 1.84 1122.40
1595.00 1.86 1185.90
1712.91 1.85 1254.30
1705.92 2.16 1246.30
1747.09 2.32 1231.60
1855.25 2.33 1298.20
1948.36 2.44 1369.70
2017.92 2.45 1438.60
2071.10 2.44 1479.40
2094.45 2.65 1474.00
2147.10 2.79 1502.60
2086.44 2.96 1479.98
2150.96 2.92 1534.68
2278.37 2.92 1639.32