*
* Example from section 9.2
*
cal(m) 1969
open data seatbelt.dat
data(format=free,skips=2,org=columns) 1969:1 1984:12 driver front rear kilos petrol
*
graph(footer="Figure 9.1 Monthly numbers (logged) of drivers who were killed\\"+$
 "or seriously injured (KSI) in road accidents in cars in Great Britain")
# driver
*
* Generate the basic seasonal and local level component models
* SeasonalDLM and LocalDLM create the "A", "C" and "SW" matrices (in RATS notation)
*
@SeasonalDLM(type=fourier) as cs sws
@LocalDLM al cl swl
*
* Glue these together to make the full model
*
compute c=cl~~cs
compute a=al~\as

nonlin phil phis
*
* The success of the estimation is quite sensitive to the initial guess values.
* The following is a bit cruder than the guess values computed using the
* reference in the book, but should work in practice
*
* (a) get rid of the seasonal by regression. (This corresponds to
*     a zero variance for that component)
*
filter(remove=seasonal) driver / deseas
*
* (b) do a short span filter on the deseasonalized data to get an
*     estimate of the local trend
*
filter(type=centered,width=5) deseas / trend
*
* (c) compute the approximate irregular component and its variance
*
set irreg = deseas-trend
stats(noprint) irreg
compute vi=%variance
*
* (d) compute the first difference of the local level component and
*     its variance
*
set difft = trend-trend{1}
stats(noprint) difft
compute vl=%variance
*
* Because we're concentrating out the irregular variance, the "phil" parameter
* needs to be initialized based upon the ratio. We'll make "phis" small.
*
compute phil=log(sqrt(vl/vi))
compute phis=-4.0
*
* Because the SW matrix for the complete model depends upon the parameters being
* estimated, we couldn't "glue" it together earlier.
*
dlm(method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
   y=driver,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=filter) / xfilter
*
*  Do smoothing using the final estimated values
*
dlm(y=driver,a=a,c=c,exact,scale,sv=1.0,sw=sw,type=smooth) / xsmooth
*
set smoothed = xsmooth(t)(1)
set filtered = xfilter(t)(1)
graph(footer="Figure 9.3. Data with filtered and smoothed estimated level",overlay=dots,ovsame) 3
# smoothed
# filtered 1970:1 *
# driver
*
* Extract the seasonal. This is done by taking the dot product of the subvector of
* the states which has the 11 seasonal components with the seasonal part of the
* "c" vector.
*
set seasonal = %dot(cs,%xsubvec(xsmooth(t),2,12))
set irreg    = driver-smoothed-seasonal
spgraph(footer="Estimated components",vfields=3)
graph(header="Level") 2
# smoothed
# driver
graph(header="Seasonal")
# seasonal
graph(header="Irregular")
# irreg
spgraph(done)
*
*  Create dummy for 1983:2 on
*
set level83_2 = t>=1983:2
*
* Estimating the DLM with a regression component. There are two ways to do this
* with RATS
*
* Method I - subtract the regression function from the dependent variable.
* Estimate the coefficients as non-linear parameters
*
nonlin phil phis b1 b2
compute b1=b2=0.0
dlm(method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
   y=driver-b1*petrol-b2*level83_2,a=a,c=c,exact,scale,sv=1.0,sw=sw) $
    / xstates vstates
*
* Method II - add the coefficients as state variables fixed over time. The
* partial DLM for this takes the form
*
*  b(t) = I b[t-1] + 0
*
* with its contribution to "c(t)" vector being equal to the regressors (petrol and
* level83_2)
*
* The system A matrix needs to have the 2x2 identity put into the bottom corner.
* C needs the regressors glued on and SW needs a 2x2 block of zeros.
*
* The new "states" should be put at the end. Add the FREEPARMS option so the log
* likelihood will be computed conditional on these, rather than unconditionally.
*
dec frml[rect] cf
frml cf = c~~petrol~~level83_2

compute a=al~\as~\%identity(2)
nonlin phil phis
dlm(method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws~\%zeros(2,2)),$
   y=driver,a=a,c=cf,exact,scale,sv=1.0,sw=sw,freeparms=2,type=smooth) / xstates vstates
*
* The final two components of xstates are the estimated coefficients based upon
* the full sample. The bottom right 2 x 2 corner of vstates will have their
* estimated variances after being scaled by the concentrated variance.
*
compute beta=%xsubvec(xstates(1984:12),13,14)
compute xx=%xsubmat(vstates(1984:12),13,14,13,14)*%variance
disp beta(1) sqrt(xx(1,1))
disp beta(2) sqrt(xx(2,2))
*
set seasonal = %dot(cs,%xsubvec(xstates(t),2,12))
set reglevel = %dot(cf,xstates(t))-seasonal
set irreg    = driver-seasonal-reglevel
*
spgraph(footer="Estimated components with intervention and regression effects",vfields=3)
graph(header="Level") 2
# reglevel
# driver
graph(header="Seasonal")
# seasonal
graph(header="Irregular")
# irreg
spgraph(done)