*
* Example 8.6 from pp 257-263
*
cal(m) 1979:5
open data table8_1.xls
data(format=xls,org=columns) 1979:5 1984:9 gas
*
graph(footer="Figure 8.2 Gas consumption data")
# gas
*
@SeasonalDLM(type=fourier) as cs sws
@LocalDLM al cl swl
*
compute c=cl~~cs
compute a=al~\as
*
* The reference analysis can be done using the options EXACT (for exact initial
* Kalman filter) and VAR=CONCENTRATE which concentrates out the observation
* variance. This is, in fact, the method of handling unobservable components
* models presented in most other textbooks. This also uses the default setting of
* SW=zero matrix, which means that the trend and seasonal components are fixed
* (unknown) constants.
*
dlm(y=gas,a=a,c=c,exact,var=concentrate,sv=1.0) / xstates vstates
*
*  Create the table of F's
*
report(action=define)
report(atrow=1,atcol=1,fillby=cols) "Harmonic" "F-value"
*
* jlower to jupper is the range of components of the state vector that covers a
* harmonic. This is a consecutive pair for all harmonics but the last one. Since
* this was estimated with a concentrated out "V", we need to divide by %VARIANCE
* to get the F.
*
do j=1,6
   compute jlower=2*j,jupper=%imin(2*j+1,12)
   compute r=%xsubvec(xstates(1984:9),jlower,jupper)
   compute v=%xsubmat(vstates(1984:9),jlower,jupper,jlower,jupper)
   compute f=%qform(inv(v),r)/((jupper-jlower+1)*%variance)
   report(atrow=1,atcol=j+1,fillby=cols) j f
end do j
report(action=format,picture="*.##")
report(action=show)
*
* Redo the model with a discount factor. Save the one-step predictions.
*
dlm(y=gas,a=a,c=c,yhat=yhat,svhat=vhat,dfhist=dfhist,$
 exact,discount=.95,$
 var=concentrate,sv=1.0) / xstates vstates
set onestep = yhat(t)(1)
set stderr  = sqrt(vhat(t)(1,1)*%variance)
set lower   = onestep-%invttest(.10,53)*stderr
set upper   = onestep+%invttest(.10,53)*stderr
graph(footer="Figure 8.3 Gas consumption and the one-step point forecasts",overlay=dots,ovsame) 4
# onestep
# lower / 2
# upper / 2
# gas
*
* Redo the model to get smoothed values. Compute and graph the components
*
dlm(y=gas,a=a,c=c,exact,discount=.95,var=concentrate,sv=1.0,type=smoothed) / xstates vstates
dec vector cx(12)
ewise cx(i)=%if(i==1,0.0,cs(i-1,1))
set mean   = %dot(cx,xstates(t))
set stderr = sqrt(%qform(vstates(t),cx)*%variance)
set lower  = mean-%invttest(.10,53)*stderr
set upper  = mean+%invttest(.10,53)*stderr
graph(style=vertical,footer="Figure 8.4 Estimated seasonal pattern in consumption series") 2
# lower
# upper
*
do j=1,6
   compute jlower=2*j
   set mean   = xstates(t)(jlower)
   set stderr = sqrt(vstates(t)(jlower,jlower)*%variance)
   set lower  = mean-%invttest(.10,53)*stderr
   set upper  = mean+%invttest(.10,53)*stderr
   graph(style=vertical,footer="Figure 8."+(4+j)+" Harmonic "+j) 2
   # lower
   # upper
end do j