*
* 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)
set january = %period(t)==1
*
* Create the structural model
*
@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
*
@LocalDLMInit(deseasonalize,irreg=sigsqeps) logksi
compute sigsqxi=sigsqeps*.01
*
* Settings for stochastic level/deterministic seasonal
*
nonlin sigsqeps sigsqxi sigsqomega=0.0
*
* If you don't need the filtered information, you can just do
* type=smooth on the DLM that does the estimation. The estimation is
* done using Kalman filtering passes, then a smoothing is done at the
* final values (only).
*
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi,$
  method=bfgs,type=smooth) / xstates vstates
*
set levelvar = vstates(t)(1,1)
graph(footer="Figure 8.1 Level estimation error variance for stochastic level\\"+$
  "and deterministic seasonal model")
# levelvar
*
set level = xstates(t)(1)
set upper = level+1.64*sqrt(levelvar)
set lower = level-1.64*sqrt(levelvar)
graph(footer="Figure 8.2 Stochastic level and 90% confidence interval") 4
# level
# upper  / 2
# lower  / 2
# logksi / 3
*
* The current seasonal is state 2
*
set seasonal = xstates(t)(2)
set upper    = seasonal+1.64*sqrt(vstates(t)(2,2))
set lower    = seasonal-1.64*sqrt(vstates(t)(2,2))
*
* This is graphed under a shorter range, since it repeats exactly, and
* over the full range, the lines run together because there's so much
* rapid vertical movement.
*
graph(footer="Figure 8.3 Deterministic seasonal and 90% confidence interval") 3
# seasonal 1981:1 1984:12
# upper    1981:1 1984:12 2
# lower    1981:1 1984:12 2
*
* The sum of the level and seasonal is the c vector dotted with the
* state vector. So the variance is the quadratic form of the variance
* matrix with c.
*
set combined = %dot(c,xstates)
set upper    = combined+1.64*sqrt(%qform(vstates,c))
set lower    = combined-1.64*sqrt(%qform(vstates,c))
graph(footer="Figure 8.4 Stochastic level + deterministic seasonal") 3
# combined 1981:1 1984:12
# upper    1981:1 1984:12  2
# lower    1981:1 1984:12  2