* PS9.9, using DATA9-3, for Application Section 9.7  
read (data9-3) period reskwh nocust price cpi incm cdd hdd pop
* create dummy variables for four periods
genr d74 = dum(period-1974)
genr d79 = dum(period-1979)
genr d83 = dum(period-1983.2)
* generate log variables  
genr kwh = reskwh/nocust
genr lkwh = log(kwh)
genr lprice = log(100*price/cpi)
genr ly = log(100*incm/(cpi*pop))
* check correlations for possible sign of multicollinearity  
stat ly lprice cdd hdd  / cor 
* estimate basic model by ols  
ols lkwh ly lprice cdd hdd / resid=ut 
* generate forecast after correcting for log model bias
genr f1 = exp(lkwh-ut+($sse/($df*2)))
* forecast error
genr abserr1 = abs(kwh - f1)
* generate lagged values for four periods  
genr ut_1=ut(-1)
genr ut_2=ut(-2)
genr ut_3=ut(-3)
genr ut_4=ut(-4)
* suppress first four observations  
sample 5 87
* estimate auxiliary regression  
ols ut ut_1 ut_2 ut_3 ut_4 ly lprice cdd hdd  
* compute nRsquared statistic  
gen1 LM = $n*$r2
* compute pvalue
distrib LM / type=chi df=4  
* reset sample range back to beginning  
sample 1 87
* estimate AR(4) model by generalized CORC  
auto lkwh ly lprice cdd hdd / order=4 resid=uhat drop
* again obtain forecasts and errors
genr f2 = exp(lkwh-uhat+($sse/($df*2)))
genr abserr2 = abs(kwh - f2)
* generate interaction terms to test for structural change  
genr lyd74=ly*d74
genr lyd79=ly*d79
genr lyd83=ly*d83
genr lprd74=lprice*d74
genr lprd79=lprice*d79
genr lprd83=lprice*d83
genr d74cdd=d74*cdd
genr d79cdd=d79*cdd
genr d83cdd=d83*cdd
genr d74hdd=d74*hdd
genr d79hdd=d79*hdd
genr d83hdd=d83*hdd
genr lprcdd=lprice*cdd
genr lprhdd=lprice*hdd
* estimate full model by ols and note that there is no autocorrelation  
ols lkwh d74 d79 d83 ly lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d74cdd d79cdd d83cdd lprcdd hdd d74hdd d79hdd d83hdd lprhdd  
* omit variable with highest pvalue, one at a time starting with d83 
ols lkwh d74 d79 ly lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d74cdd d79cdd d83cdd lprcdd hdd d74hdd d79hdd d83hdd lprhdd  
* omit d74  
ols lkwh d79 ly lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d74cdd d79cdd d83cdd lprcdd hdd d74hdd d79hdd d83hdd lprhdd  
* omit d79cdd  
ols lkwh d79 ly lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d74cdd d83cdd lprcdd hdd d74hdd d79hdd d83hdd lprhdd  
* omit d79hdd  
ols lkwh d79 ly lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d74cdd d83cdd lprcdd hdd d74hdd d83hdd lprhdd  
* omit ly  
ols lkwh d79 lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d74cdd d83cdd lprcdd hdd d74hdd d83hdd lprhdd  
* omit d74cdd  
ols lkwh d79 lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd d83cdd lprcdd hdd d74hdd d83hdd lprhdd  
* omit d83cdd  
ols lkwh d79 lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd lprcdd hdd d74hdd d83hdd lprhdd  
* omit lprhdd  
ols lkwh d79 lyd74 lyd79 lyd83 lprice lprd74 lprd79 &
 lprd83 cdd lprcdd hdd d74hdd d83hdd / resid=uhat
* generate final model forecasts and errors*
genr f3 = exp(lkwh-uhat+($sse/($df*2)))
genr abserr3 = abs(kwh - f3)
* exclude four initial observations and compute percent errors, mean 
* absolute percent error, and mean squared error
sample 5 87
genr pcterr1 = 100*abserr1/kwh
genr pcterr2 = 100*abserr2/kwh
genr pcterr3 = 100*abserr3/kwh
genr mape1=sum(pcterr1)/83
genr mape2=sum(pcterr2)/83
genr mape3=sum(pcterr3)/83
genr mse1=sum(pcterr1*pcterr1)/83
genr mse2=sum(pcterr2*pcterr2)/83
genr mse3=sum(pcterr3*pcterr3)/83
print period pcterr1 pcterr2 pcterr3 
sample 87 87
print mape1 mape2 mape3 mse1 mse2 mse3 
stop