*
* Example 8.5.2 from pp 281-283
*
cal(m) 1949
open data airpass.dat
data(format=free) 1949:1 1960:12 airpass
*
graph(footer="Figure 8-3 International airline passengers; monthly totals")
# airpass
*
* The model is y(t)=c(t)+s(t)+e(t) where the
* trend term c(t) is given by:
*  c(t)=c(t-1)+tau(t-1)+sig1(t),
*  tau(t)=tau(t-1)+sig2(t)
* and the seasonal s(t) by s(t)+s(t-1)+...+s(t-11)=sig3
*
* The A and C matrices are fixed. We use the LocalDLM and SeasonalDLM
* procedures to create the trend and seasonal parts and combine them.
*
@LocalDLM(type=trend) at ct
@SeasonalDLM(type=additive) as cs
compute a=at~\as
compute c=ct~~cs
*
* The state equation variance matrix is created as
* a FRML so the variances can be varied easily.
*
dec real sig1 sig2 sig3 sigv
dec frml[symm] sw
frml sw = %diag(||sig1,sig2,sig3,0,0,0,0,0,0,0,0,0,0||)
*
* Rather than estimate initial conditions, we use the "exact" option to allow
* for a diffuse prior. As is typical of models like this, freely estimating the
* component variances tends to produce an overly optimistic "fit" by having the
* trend term adjust to fit the data almost perfectly, while forcing the
* observation equation variance to zero. sigv is pegged to a small number and the
* model is estimated, concentrating out a scale factor on the variances.
*
compute sigv=.000001,sig2=0.0,sig1=sig3=10.0
nonlin sig1 sig2 sig3 sig1>=0.0 sig2=0.0 sig3>=0.0
dlm(pmethod=simplex,piters=10,method=bfgs,a=a,c=c,sw=sw,sv=sigv,$
  y=airpass,exact,var=concentrate,yhat=yhat) / xstates
*
* These are the components off the KF with estimated parameters. xstates(t) is the
* state vector at time t, so we want to pull out the 1st and 3rd components of
* it. The second component is the slope component, which will be (almost)
* constant across the data set, since its variance is estimated at zero.
*
set seasonal = xstates(t)(3)
set slope    = xstates(t)(2)
set trend    = xstates(t)(1)
*
* Graph the components
*
graph(footer="Figure 8-4 State components (contemporaneous)",key=upleft) 3
# seasonal
# trend
# slope
*
set onestep = yhat(t)(1)
graph(footer="Figure 8-5 One-step predictors and actual data",overlay=dots,ovsame) 2
# onestep
# airpass
*
* This redoes the model with estimated initial conditions as is done in the book.
* Note that RATS uses a different labelling scheme for the state equation with
* X(t)=FX(t-1)+V(t) instead of X(t+1)=FX(t)+V(t). Ordinarily, this has little
* consequence, but here, the difference is that the initial state for the RATS
* formulation would be X(0) rather than X(1). The means are related by X(1)=FX(0).
* However, a zero variance for X(0) isn't the same as a zero variance for X(1).
* In order to come close to matching up the results in the text, it's necessary
* to back out the solution for Omega(0) in AOmega(0)A'+Q=0.
*
dec vect x0(13)
dec frml[vect] x0f
dec frml[symm] sx0f
frml x0f = x0
frml sx0f = -1.0*inv(a)*sw*tr(inv(a))
nonlin sig1 sig2 sig3 sig2=0.0 x0
compute sigv=.000001,sig2=0.0,sig1=sig3=10.0
dec vect x1(13)
input x1
 146.9 2.171 -34.92 -34.12 -47.00 -16.98 22.99 53.99 58.34 33.65 2.204 -4.053 -6.894
compute x0=inv(a)*x1
*
dlm(pmethod=simplex,piters=10,method=bfgs,iters=100,$
   a=a,c=c,sw=sw,sv=sigv,y=airpass,x0=x0f,sx0=sx0f,var=conditional) / xstates