*
*  Illustration from section 2.2.2
*
open data nile.dat
calendar 1871
data(format=free,org=columns,skips=1) 1871:1 1970:1 nile
*
* Note that there is an important difference between DLM and the description in
* the book which will affect some of the calculations. The RATS definition of the
* state-space model is
*
*   x(t) = A x(t-1) + w(t)
*   y(t) = c x(t) + v(t)
*
* (The options on the DLM instruction take their names from the component names in
* this description)
*
* The book uses
*
*   a(t+1) = T a(t) + eta(t)
*   y(t)   = Z a(t) + eps(t)
*
* This is not just an inconsequential relabeling of matrices. Note the difference
* in the subscripts in the state equation (the first in each pair). The two forms
* end up having identical behavior, but the book's model has a(t+1) being the
* forecast of x(t+1) given information through time t. While the forecast of y
* and its predictive variance are the same with either model, the actual state
* variables and variances aren't.
*
* This is a very simple set up for DLM since there is only one state variable,
* and one observable, and all the variances are treated as known.
*
dlm(y=nile,a=1.0,c=1.0,sx0=1.e+7,sv=15099.0,sw=1469.1,vhat=vhat,svhat=fhat) / xstates vstates
*
* Because DLM can handle multiple inputs (vector Y) and because most models have
* more than one state, the outputs generally are series of VECTORs for states,
* errors and observables and series of SYMMETRICs for variances. To use these in
* further calculations, you usually will have to extract the information into
* simpler series of real numbers. That's what the SET instructions are doing. Note
* that the two that are pulling out variances (from vstates which gives the
* posterior variances of the states and from fhat which are the prediction error
* variances) request the 1,1 element (since those will, in general be NxN
* matrices), while the others are N-vectors and only need the (1) subscript.

set a = xstates(t)(1)
set p = vstates(t)(1,1)
set v = vhat(t)(1)
set f = fhat(t)(1,1)
set lower = a-sqrt(p)*%invnormal(.95)
set upper = a+sqrt(p)*%invnormal(.95)
*
* The graphs all leave out the first data point (1871) which will be poorly
* estimated and will show up as a huge outlier on the state and variance graphs.
*
spgraph(footer="Figure 2.1. Nile data and output of Kalman filter",vfields=2,hfields=2)
graph(header="Filtered State and 90% Confidence Intervals") 4
# a 1872:1 *
# lower 1872:1 * 2
# upper 1872:1 * 2
# nile 1871:1 * 3
graph(header="Prediction Errors")
# v 1872:1 *
graph(header="Filtered State Variance")
# p 1872:1 *
graph(header="Prediction Variance")
# f 1872:1 *
spgraph(done)