*
* Chapter 7. 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
# constant logukpetrol seatbelt 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
*
* Model with both explanatory variables with fixed coefficients,
* deterministic level and seasonal. It will be easier to work with this
* if we create a FRML for the explanatory variable component.
*
frml explan = beta*logukpetrol+lambda*seatbelt
*
nonlin sigsqeps sigsqxi=0.0 sigsqomega=0.0 beta lambda
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  method=bfgs,vhat=vhat,svhat=svhat)
*
set resids = %scalar(vhat)/sqrt(%scalar(svhat))
@STAMPDiags(ncorrs=15) resids
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  type=smooth) / xstates
set level     = xstates(t)(1)
set levelplus = level+explan
graph(footer="Figure 7.1 Deterministic level plus explanatory variables",$
  key=upright,klabels=||"log UK drivers KSI","deterministic level+beta*log(petrol)+lambda*seatbelt"||) 2
# logksi
# levelplus
*
* Same model with variances on the level and seasonal freed.
*
nonlin sigsqeps sigsqxi sigsqomega beta lambda
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  method=bfgs,vhat=vhat,svhat=svhat)
*
* Since that doesn't work well, same with deterministic seasonal
*
nonlin sigsqeps sigsqxi sigsqomega=0.00 lambda beta
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  method=bfgs,vhat=vhat,svhat=svhat)
*
set resids = %scalar(vhat)/sqrt(%scalar(svhat))
@STAMPDiags(ncorrs=15) resids
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi-explan,$
  type=smooth) / xstates
set level     = xstates(t)(1)
set levelplus = level+explan
graph(footer="Figure 7.2 Stochastic level plus explanatory variables",$
  key=upright,klabels=||"log UK drivers KSI","stochastic level+beta*log(petrol)+lambda*seatbelt"||) 2
# logksi
# levelplus