*
* Chapter 8. UK data.
*
open data ukdriversksi.txt
calendar(m) 1969
data(format=free,org=columns,skips=1) 1969:01 1984:12 ksi
set logksi = log(ksi)
open data logukpetrolprice.txt
data(format=free,org=columns,skips=1) 1969:01 1984:12 logukpetrol
set seatbelt = t>=1983:2
*
* Fixed coefficients regression done with linreg
*
seasonal seasons
linreg logksi * 1983:1
# constant logukpetrol seasons{0 to -10}
*
compute sigsqeps=%sigmasq,beta=%beta(2),lambda=%beta(3)
*
* Create the standard structural model for local level + seasonal
*
@LocalDLM(type=level,a=al,c=cl,f=fl)
@SeasonalDLM(type=additive,a=as,c=cs,f=fs)
compute a=al~\as,f=fl~\fs,c=cl~~cs
*
frml explan = beta*logukpetrol
*
* Estimates through Jan 1983 (prior to seatbelt law) using stochastic
* level, deterministic seasonal and the log price of petrol as
* explanatory variable.
*
compute sigsqxi=.01*sigsqeps
nonlin sigsqeps sigsqxi sigsqomega=0.00 beta
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  method=bfgs) * 1983:1
*
* Using those estimates, forecast the last 23 periods of data, holding
* back the actual values of logksi, but using the petrol price.
*
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,$
  y=%if(t<=1983:1,logksi-explan,%na),yhat=yhat,svhat=svhat) / xstates vstates
*
set forecast  1970:1 * = explan+%if(t<=1983:1,%scalar(yhat),%dot(c,xstates))
set fvariance 1970:1 * = %if(t<=1983:1,%scalar(svhat),%qform(vstates,c))
*
set upper = forecast+1.64*sqrt(fvariance)
set lower = forecast-1.64*sqrt(fvariance)
*
graph(footer="Figure 8.15 Forecasts out of sample") 3
# forecast 1983:2 *
# upper    1983:2 * 2
# lower    1983:2 * 2
*
* Redo estimation with all data, including the intervention variable
*
frml explan = beta*logukpetrol+lambda*seatbelt
nonlin sigsqeps sigsqxi sigsqomega=0.00 beta lambda
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,$
  method=bfgs,y=logksi-explan,yhat=yhat,svhat=svhat) / xstates vstates
*
set cforecast 1970:1 * = %scalar(yhat)+explan
*
graph(footer="Figure 8.16 Last four years of data with forecasts obtained using\\"+$
 "data through Jan 1983 and with intervention") 3
# logksi    1982:1 *
# forecast  1982:1 *
# cforecast 1982:1 *