*
* Chapter 4. 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
graph(footer="Figure 4.1 Log UK drivers KSI with time lines for years",$
  grid=january)
# logksi
*
* Create the component models for the level and seasonal.
*
@LocalDLM(type=level,a=al,c=cl,f=fl)
@SeasonalDLM(type=additive,a=as,c=cs,f=fs)
*
* Glue them together to make the full system matrices
*
compute a=al~\as,f=fl~\fs,c=cl~~cs
*
@LocalDLMInit(deseasonalize,irreg=sigsqeps) logksi
*
* First model, deterministic level and seasonal pegs the variances for
* xi and omega to zero
*
nonlin sigsqeps sigsqxi=0.0 sigsqomega=0.0
*
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi,$
  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,$
  type=smoothed) / xstates
set level    = xstates(t)(1)
set irreg    = logksi-level
set seasonal = xstates(t)(2)
set fitted   = %dot(c,xstates)
*
graph(footer="Figure 4.2 Combined deterministic level and seasonal",$
 key=upright,klabels=||"log UK drivers KSI","deterministic level plus seasonal"||) 2
# logksi
# fitted
*
graph(footer="Figure 4.3 Deterministic level and seasonal",$
 key=upright,klabels=||"log UK drivers KSI","deterministic level"||) 2
# logksi
# level
graph(footer="Figure 4.4 Deterministic seasonal",$
 key=upright,klabels=||"deterministic seasonal"||)
# seasonal
graph(footer="Figure 4.5 Irregular component",$
 key=upright,klabels=||"irregular"||)
# irreg
*
* Same model with variances freed
*
nonlin sigsqeps sigsqxi sigsqomega
*
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi,$
  method=bfgs,vhat=vhat,svhat=svhat)
*
* The variance of the seasonal comes in negative, so we re-estimate with
* that set to zero. This is the same as the model in section 4.3.
*
nonlin sigsqeps sigsqxi sigsqomega=0.00
dlm(a=a,c=c,sv=sigsqeps,f=f,sw=%diag(||sigsqxi,sigsqomega||),exact,y=logksi,$
  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,$
  type=smooth) / xstates
set level    = xstates(t)(1)
set irreg    = logksi-level
set seasonal = xstates(t)(2)
set fitted   = %dot(c,xstates)
*
graph(footer="Figure 4.6 Stochastic level",$
 key=upright,klabels=||"log UK drivers KSI","stochastic level"||) 2
# logksi
# level
graph(footer="Figure 4.7 Stochastic seasonal",$
 key=upright,klabels=||"stochastic seasonal"||)
# seasonal
graph(footer="Figure 4.8 Stochastic seasonal for the year 1969",$
 key=upleft,klabels=||"seasonal 1969"||)
# seasonal 1969:1 1969:12
graph(footer="Figure 4.9 Irregular component",$
 key=upright,klabels=||"irregular"||)
# irreg