*
*  Gas consumption in UK example from section 14.3
*
open data gas.dat
cal(q) 1960
data(format=free,skip=1) 1 107 gas
*
@SeasonalDLM(type=additive) as cs sws
@LocalDLM(type=trend) al cl swl
*
nonlin phil phis
*
compute phil=-2.0,phis=-4.0
compute c=cl~~cs
compute a=al~\as
dlm(method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
   y=gas,a=a,c=c,exact,var=concentrate,sv=1.0,sw=sw,vhat=vhat,type=smooth) / xstates vstates
*
*  Pull out seasonal and irregular
*
set seasonal = xstates(t)(3)
set irreg    = vhat(t)(1)
spgraph(footer="Fig 14.4 Analyses of gas data based upon Gaussian model",vfields=2)
graph(header="Gaussian conditional seasonal")
# seasonal
graph(header="Gaussian conditional irregular")
# irreg
spgraph(done)
*
compute nu=50.0
compute veps=%variance
set eps = 0.0
compute lastone=0
do iters=1,20
   compute lastnu=nu
   set fiddle = ((nu-2)+eps**2/veps)/(nu+1)
   dlm(print=lastone,method=bfgs,start=(sw=exp(2*phil)*swl~\exp(2*phis)*sws),$
      y=gas,a=a,c=c,exact,var=concentrate,sv=fiddle,sw=sw,vhat=vhat,type=smooth) / xstates vstates
   if lastone
      break
   set eps = vhat(t)(1)
   stats(noprint) eps
   *
   * Use method of moments to estimate nu and sigma epsilon
   *
   compute nu=(zf=(%kurtosis+3)/3.0),(4*zf-2)/(zf-1)
   compute veps=%variance
   compute lastone=fix((abs(nu-lastnu)<.01))
end do iters
*
set seasonal = xstates(t)(3)
set irreg    = vhat(t)(1)
spgraph(footer="Fig 14.5 Analyses of gas data based upon t-model ",vfields=2)
graph(header="t conditional seasonal")
# seasonal
graph(header="t conditional irregular")
# irreg
spgraph(done)