*
*  Illustration from section 6.6.3.
*  Smoothing Spline model
*
open data nile.dat
calendar 1871
data(format=free,org=columns,skips=1) 1871:1 1970:1 nile
*
* Set the 1890-1900 and 1950-1960 ranges to NA
*
set nile 1890:1 1900:1 = %na
set nile 1950:1 1960:1 = %na
*
* The local linear trend model with the timing properties used by RATS is:
*
*   y(t)=mu(t)+eps(t)
*   mu(t)=mu(t-1)+tau(t-1)
*   tau(t)=tau(t-1)+eta(t)
*
* The states are mu(t) and tau(t).
*
*  y(t)=[1,0][mu(t)|tau(t)]+eps(t)
*
*  [mu(t) ]    [1  1] [mu(t-1)]      [0]
*  [tau(t)] =  [0  1] [tau(t-1)]  +  [eta(t)]
*
* lambda=2500 would be accomplished by pegging the variance of eta at
* 1/2500=.0004. Because the graphs are done with the (incorrect) .004, we will
* use that as well. The only "free" parameter is the variance scale factor, which
* is concentrated out. The "SCALE" option is used to indicate that that should be
* done. Because the states are fully non-stationary, we can just use the EXACT
* option to handle the initial conditions.
*
* Do filter first
*
dlm(a=||1.0,1.0|0.0,1.0||,c=||1.0,0.0||,sv=1.0,sw=||0.0|0.0,.004||,y=nile,$
     exact,scale,type=filter) / xstates vstates
set fit   = xstates(t)(1)
set vfit  = %variance*vstates(t)(1,1)
set upper = fit + sqrt(vfit)*2.0
set lower = fit - sqrt(vfit)*2.0
spgraph(footer="Figure 6.1. Spline Estimates",vfields=2)
graph(header="Filtered Estimate with 95% Confidence Interval") 4
# nile
# fit
# upper / 3
# lower / 3
*
* Now smooth
*
dlm(a=||1.0,1.0|0.0,1.0||,c=||1.0,0.0||,sv=1.0,sw=||0.0|0.0,.004||,y=nile,$
   exact,scale,type=smooth) / xstates vstates
set fit   = xstates(t)(1)
set vfit  = %variance*vstates(t)(1,1)
set upper = fit + sqrt(vfit)*2.0
set lower = fit - sqrt(vfit)*2.0
graph(header="Smoothed Estimate with 95% Confidence Interval") 4
# nile
# fit
# upper / 3
# lower / 3
spgraph(done)