*
* Example 14.2 from pp 230-233
*
cal(m) 1969
open data van.dat
data(format=free,skips=1,org=columns) 1969:1 1984:12 van level83_2
*
@SeasonalDLM(type=additive) as cs sws
@LocalDLM(type=level) al cl swl
*
* This is used to get guess values for the variances of the components for a
* Gaussian linear model. The sequence of calculations is described in
* durkp162.prg.
*
nonlin phil phis
filter(remove=seasonal) van / deseas
filter(type=centered,width=5) deseas / trend
set irreg = deseas-trend
stats(noprint) irreg
compute vi=%variance
set difft = trend-trend{1}
stats(noprint) difft
compute vl=%variance
compute phil=log(sqrt(vl/vi))
compute phis=-4.0
*
* This includes the regression coefficient as a freely estimated state variable.
*
dec frml[rect] cf
frml cf = cl~~cs~~level83_2
compute a=al~\as~\1.0
dlm(method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws~\0.0),$
   y=van,a=a,c=cf,exact,scale,sv=1.0,sw=sw,free=1) / xstates vstates
*
* Poisson model
* Start out with an expansion around theta-tilde = log(mean of van). Iterate
* until the change in the log likelihood is small
*
nonlin phil phis
set theta = log(9.0)
compute log0=%na,lastone=0
do iters=1,10
   dlm(print=lastone,method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws~\0.0),$
     y=van-exp(theta)+exp(theta)*theta,a=a,c=cf*exp(theta),free=1,exact,sv=exp(theta),sw=sw,type=smooth) / xstates vstates
   if lastone
      break
   set theta = %dot(xstates(t),cf)
   compute lastone=fix((%logl-log0)<1.e-4),log0=%logl
end do iters
*
set seasonal = %dot(cs,%xsubvec(xstates(t),2,12))
set levelreg = theta-seasonal
set mean     = exp(levelreg)
*
spgraph(footer="Fig 14.1",vfields=2)
graph(header="Number of van drivers killed and estimated level including intervention") 2
# mean
# van
set seasadj = van/exp(seasonal)
graph(header="Seasonally adjusted number of van drivers killed and estimated level") 2
# mean
# seasadj
spgraph(done)