options crt;  ? 7e8
in nerc;

title 'Exercise 8, part (a) - stationarity check';
freq a;
smpl 52,84;
lnk1 = log(kwh/kwh(-1));
lnp1 = log(pelec/pelec(-1));
lng1 = log(gnp/gnp(-1));
smpl 52,73;
dot k p g;
  bjident(nlag=5,nlagp=5,ndiff=2,noplot) ln.1;
enddot;

? create a few global variables for use after bcd exits
const alpha,beta1,beta2; lnp1f=0; lng1f=0;
bcd 73 1;
proc bcd last ifp;
?
?   Do parts (b), (c) and (d) for an arbitrary last year.
?
if last=73 & ifp; then;
title 'Exercise 8, part (b) - AR(1) regression';
? Actually no observations are lost for AR(1) ML -- the only
? observation lost is 1951, due to differencing.
?smpl 54,73;
smpl 52,last;
if (ifp); then; do;
  ar1(silent) lnk1 c lnp1 lng1;
  unmake @coef alpha beta1 beta2 rho;  ? assuming TSP 4.5
  enddo;
else;
  ar1(silent) lnk1 c      lng1;
copy @coef bar1;
? The above are short-run elasticities.
title 'long-run elasticities';
? Without reference to a specific formula giving the
? expectations process, it is not clear how to derive the
? long-run elasticities.  However, one possibility is that
? phi = rho.
set rho = @rho;  ? for TSP 4.4 and earlier
frml ephi phi = rho;
frml eb1 b1lr = lnp1/rho;
frml eb2 b2lr = lng1/rho;
if (ifp); then;
  analyz ephi eb1 eb2;  ? phi and long-run elasticities
else;
  analyz ephi     eb2;  ? phi and long-run elasticities

set fy1 = last+1;
set fy10 = last+10;
if last=73 & ifp; then;
title 'Exercise 8, part (c) - ARIMA(0,1,1) forecasts of RHS vars';
dot p g;
  smpl 53,last;
  bjest(nar=0,ndiff=1,nma=1,const,exactml,grad=c4,silent) ln.1;
  bjfrcst(nhoriz=10,noprint) ln.1;
  smpl fy1,fy10;
  ln.1f = @fit;
enddot;

if last=73 & ifp; then;
title 'Exercise 8, part (d) - growth rate forecasts';
smpl last,last;  ? initialize RHS vars in presample period for AR(1) forecast
lnp1f = lnp1;
lng1f = lng1;
smpl fy1,fy10;
if (ifp); then;
  forcst(print,rho=rho,coef=bar1,depvar=lnk1) lnkf c lnp1f lng1f;
else;
  forcst(print,rho=rho,coef=bar1,depvar=lnk1) lnkf c       lng1f;
msd lnkf;
endproc;

title 'Exercise 8, part (e) - growth rates of RHS vars';
smpl 73,73;
dot p g;
  aagr.  = exp[(ln.1(+10) - ln.1)/10] - 1;
  aagr.f = exp[(ln.1f(+10) - ln.1)/10] - 1;
enddot;
print aagrp aagrpf aagrg aagrgf;
title 'forecasted demand growth using actual RHS var growth';
?  (and ignoring rho)
lnkfa = alpha + beta1*aagrp + beta2*aagrg;
print lnkfa;

title 'Exercise 8, part (f) - varying the last year';
do last=74,84;
  print last;
  bcd last 1;
enddo;

title 'Exercise 8, part (g) - without price, varying the last year';
do last=73,84;
  print last;
  bcd last 0;
enddo;