*
* Example 11.3 from pp 535-540
*
open data jnj.dat
calendar(q) 1960
data(format=free,org=columns) 1960:1 1980:4 jnj
set jnj = log(jnj)
*
* This generates the system matrices for the local level model and the additive
* seasonal.
*
@localdlm(type=level) al cl swl
@seasonaldlm(type=additive) as cs sws
*
* This concatenates the system "A" and "C" matrices. "A" is a block diagonal
* concatenation, while "C" is constructed by placing one over the other.
*
compute a=al~\as
compute c=cl~~cs
*
* Again, we'll concentrate out a scale factor. It turns out that in this case,
* we'd be better off pegging one of the state variances because we're normalizing
* on a number which is coming in close to zero, but the results end up being
* almost identical regardless.
*
nonlin leta lomega
compute leta=lomega=0.0
*
dlm(method=bfgs,y=jnj,a=a,c=c,start=(sw=exp(leta)*swl~\exp(lomega)*sws),sw=sw,$
   sv=1.0,var=concentrate,exact,type=smooth) / xstates vstates
*
* Note that the estimates shown in the text are incorrect. The GetSsfStsm function
* takes standard deviations, not variances as inputs, so the jnjest definition
* shouldn't include the sqrt(...). Because those are then fed back into the model,
* the upper to lower spreads on the subsequent graphs are wider than they should
* be.
*
disp "Component Std Dev"
disp "Measurement" @20 sqrt(%variance)
disp "Level" @20 sqrt(%variance*exp(leta))
disp "Seasonal" @20 sqrt(%variance*exp(lomega))
*
spgraph(footer="Figure 11.6 Smoothed components of fitting 11.85 to log of\\quarterly earnings per share of Johnson & Johnson",vfields=2)
set trend = xstates(t)(1)
set lower = trend+%invnormal(.025)*sqrt(vstates(t)(1,1)*%variance)
set upper = trend+%invnormal(.975)*sqrt(vstates(t)(1,1)*%variance)
graph(header="Trend Component") 3
# trend
# lower 5 * 2
# upper 5 * 2
*
set seasonal = xstates(t)(2)
set lower = seasonal+%invnormal(.025)*sqrt(vstates(t)(2,2)*%variance)
set upper = seasonal+%invnormal(.975)*sqrt(vstates(t)(2,2)*%variance)
graph(header="Seasonal Component") 3
# seasonal
# lower 5 * 2
# upper 5 * 2
spgraph(done)