*
*  Illustration from section 2.6.1
*  Simulation and conditional simulation
*
open data nile.dat
calendar 1871
data(format=free,org=columns,skips=1) 1871:1 1970:1 nile
*
dlm(type=smooth,y=nile,a=1.0,c=1.0,sx0=1.e+7,sv=15099.0,sw=1469.1,$
  vhat=vhat,what=what) / xstates
dlm(type=simulate,yhat=nilehat,a=1.0,c=1.0,sx0=1.e+7,sv=15099.0,sw=1469.1) / xstates_s
dlm(type=csimulate,y=nile,a=1.0,c=1.0,sx0=1.e+7,sv=15099.0,sw=1469.1,what=wdraws,vhat=vdraws) / xstates_c
*
* Because the initial state for the unconditional sample is random with a huge
* variance, it wouldn't work to put the direct draw for the states on the same
* graph as the smoothed states. Instead, we peg the sampled states to start at the
* same value as the smoothed states (by subtracting off the initial difference).
*
set smoothed = xstates(t)(1)
set simulate = xstates_s(t)(1)-xstates_s(1871:1)(1)+smoothed(1871:1)
set csimulate = xstates_c(t)(1)
spgraph(footer="Figure 2.4 Simulations",vfields=2,hfields=2)
graph(overlay=dots,ovsame,header="Smoothed State and Unconditional Sample") 2
# smoothed
# simulate
graph(overlay=dots,ovsame,header="Smoothed State and Conditional Sample") 2
# smoothed
# csimulate
*
set smoothedobsv   = vhat(t)(1)
set sampleobsv     = vdraws(t)(1)
graph(overlay=dots,ovsame,header="Smoothed Observation Error and Sample") 2
# smoothedobsv
# sampleobsv
*
set smoothedstatew = what(t)(1)
set samplestatew   = wdraws(t)(1)
graph(overlay=dots,ovsame,header="Smoothed State Error and Sample") 2
# smoothedstatew
# samplestatew
spgraph(done)