*
* Chapter 9
*
open data ukfrontrearseatksi.txt
calendar(m) 1969
data(format=free,org=columns,skips=1) 1969:01 1984:12 driver front rear kilos petrol
set seatbelt = t>=1983:2
*
graph(footer="Figure 9.1. Front and rear seat passengers killed and seriously injured\\"+$
 "in the UK in the period 1969-1984",key=attached) 2
# front
# rear
*
* Generate the basic seasonal and local level component models.
*
@SeasonalDLM(type=additive,a=as,c=cs,f=fs)
@LocalDLM(a=al,c=cl,f=fl)
compute a1=(al~\as)
compute f1=(fl~\fs)
compute c1=(cl~~cs)
*
* Because the two series will have separately estimated level and
* seasonal components, we have to "double" all the component matrices
* up: the A matrix will be 2N x 2N, the F matrix will be 2N x 2M and the
* C matrix will be 2N x 2. There are several ways to organize the
* combined states. This one is most convenient, because it keeps the
* contemporaneously correlated shocks in consecutive positions. What
* this does is to interleave the states and the shocks.
*
dec rect a(%rows(a1)*2,%rows(a1)*2)
ewise a(i,j)=%if(%clock(i,2)==%clock(j,2),a1((i+1)/2,(j+1)/2),0.0)
dec rect f(%rows(f1)*2,%cols(f1)*2)
ewise f(i,j)=%if(%clock(i,2)==%clock(j,2),f1((i+1)/2,(j+1)/2),0.0)
dec rect c(%rows(a1)*2,2)
ewise c(i,j)=%if(%clock(i,2)==j,c1((i+1)/2,1),0.0)
*
* The component variances will now be 2 x 2 symmetric matrices
*
dec symm sigmaeps(2,2) sigmaxi(2,2) sigmaomega(2,2)
dec symm sw
*
* sigmaxi and sigmaomega are modeled in packed lower triangular form
* since it's quite possible for one of them to go non-positive definite.
*
dec packed pxi(2,2) pomega(2,2)
*
* This function takes the current packed forms and converts them into
* standard covariance matrices, then concatenates them diagonally to
* form the overall covariance matrix for the shocks.
*
function %%DLMSetupSW
type symmetric %%DLMSetupSW
compute sigmaxi=%ltouterxx(pxi),sigmaomega=%ltouterxx(pomega)
compute %%DLMSetupSW=sigmaxi~\sigmaomega
end %%DLMSetupSW
*
compute sigmaeps=%mscalar(.005),pxi=%mscalar(.0003),pomega=%mscalar(0.0)
*
dec vect bpetrol(2) bkilos(2) blaw(2)
dec vect[frml] explan(2)
frml explan(1) = bpetrol(1)*petrol+bkilos(1)*kilos+blaw(1)*seatbelt
frml explan(2) = bpetrol(2)*petrol+bkilos(2)*kilos+blaw(2)*seatbelt
*
* Estimate the unconstrained model
*
nonlin sigmaeps pxi pomega bpetrol bkilos blaw
dlm(start=(sw=%%DLMSetupSW()),sw=sw,a=a,c=c,f=f,sv=sigmaeps,$
  y=||front-explan(1),rear-explan(2)||,exact,$
  method=bfgs) / xstates
*
* Re-estimate with the seasonal variance zeroed out
*
compute pomega=%zeros(2,2)
nonlin sigmaeps pxi bpetrol bkilos blaw
dlm(start=(sw=%%DLMSetupSW()),sw=sw,a=a,c=c,f=f,sv=sigmaeps,$
  y=||front-explan(1),rear-explan(2)||,exact,$
  method=bfgs,vhat=vhat,svhat=svhat)
*
* This standardizes the residuals using the inverse Choleski factor (a
* matrix generalization of the reciprocal square root). We can apply
* standard diagnostics to each component of this.
*
dec vect stdres
set r1 = stdres=inv(%decomp(svhat))*vhat,stdres(1)
set r2 = stdres=inv(%decomp(svhat))*vhat,stdres(2)
@STAMPDiags r1
@STAMPDiags r2
*
disp "H=" #.####### sigmaeps
disp "Q=" #.####### sigmaxi
*
dlm(start=(sw=%%DLMSetupSW()),sw=sw,a=a,c=c,f=f,sv=sigmaeps,$
  y=||front-explan(1),rear-explan(2)||,exact,$
  type=smoothed,what=what) / xstates
*
set lshock1 = what(t)(1)
set lshock2 = what(t)(2)
linreg lshock1
# constant lshock2
scatter(footer="Figure 9.2 Level disturbances for rear vs front seat",$
  line=%beta)
# lshock2 lshock1
*
set level1 = xstates(t)(1)
set level2 = xstates(t)(2)
spgraph(vfields=2,samesize,$
 footer="Figure 9.3 Levels of treatment and control series in the SUR model")
graph(key=upright,klabels=||"level front"||)
# level1
graph(key=upright,klabels=||"level rear"||)
# level2
spgraph(done)
*
scatter(footer="Figure 9.4 Level of treatment vs level of control")
# level2 level1
*
* Re-estimation, restricting level disturbances to rank one and
* eliminating coefficient on the seatbelt term.
*
nonlin sigmaeps pxi bpetrol bkilos blaw blaw(2)=0.0 pxi(2,2)=0.0
dlm(start=(sw=%%DLMSetupSW()),sw=sw,a=a,c=c,f=f,sv=sigmaeps,$
  y=||front-explan(1),rear-explan(2)||,exact,$
  method=bfgs)
dlm(start=(sw=%%DLMSetupSW()),sw=sw,a=a,c=c,f=f,sv=sigmaeps,$
  y=||front-explan(1),rear-explan(2)||,exact,$
  type=smoothed,what=what) / xstates
*
set level1 = xstates(t)(1)
set level2 = xstates(t)(2)
scatter(footer="Figure 9.6 Level of treatment vs level of control")
# level2 level1
*
spgraph(vfields=2,samesize,$
 footer="Figure 9.7 Levels of treatment and control series in the rank one model")
graph(key=upright,klabels=||"level front"||)
# level1
graph(key=upright,klabels=||"level rear"||)
# level2
spgraph(done)
*
set levelplus1 = level1+blaw(1)*seatbelt
spgraph(vfields=2,samesize,$
 footer="Figure 9.8 Levels of treatment plus intervention\\and control series in the rank one model")
graph(key=upright,klabels=||"level front"||)
# levelplus1
graph(key=upright,klabels=||"level rear"||)
# level2
spgraph(done)
*
set seas1 = xstates(t)(3)
set seas2 = xstates(t)(4)
spgraph(vfields=2,samesize,$
  footer="Figure 9.9 Deterministic seasonal of treatment and control series, rank one model")
graph(key=upleft,klabels=||"seasonal front"||)
# seas1
graph(key=upleft,klabels=||"seasonal rear"||)
# seas2
spgraph(done)