*
* Chapter 8. UK data. Missing values.
*
open data ukdriversksi.txt
calendar(m) 1969
data(format=free,org=columns,skips=1) 1969:01 1984:12 ksi
set logksi = log(ksi)
*
@LocalDLM(type=level,a=al,c=cl,f=fl) r
@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
*
* Stochastic level, deterministic seasonal.
*
* Patch over series with the range that we want to treat as missing
*
set withmiss = %if(t>=48.and.t<=62.or.t>=120.and.t<=140,%na,logksi)
*
nonlin sigsqeps sigsqxi sigsqomega=0.00
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=withmiss,$
  method=bfgs,type=smooth) / xstates vstates
*
set levelvar = vstates(t)(1,1)
graph(footer="Figure 8.17 Stochastic level estimation error variance\\"+$
 "with missing values")
# levelvar
*
set level = xstates(t)(1)
set upper = level+1.64*sqrt(levelvar)
set lower = level-1.64*sqrt(levelvar)
graph(footer="Figure 8.18 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.20 Deterministic seasonal and 90% confidence interval") 3
# seasonal 1971:1 1975:12
# upper    1971:1 1975:12 2
# lower    1971:1 1975:12 2
*
* The sum of the level and seasonal is the c vector dotted with the
* state vector. The irregular will be the difference between the data
* and that value. Since the data are missing in certain ranges, the
* irregular will be NA there as well.
*
set combined = %dot(c,xstates)
set irreg    = withmiss-combined
graph(footer="Figure 8.21 Irregular component")
# irreg