*
* Example 2.1, pp 40-43
*
data(format=free,unit=input) 1 10 kurit
 150 136 143 154 135 148 128 149 146 .
*
* Local level model
*
* The state equation is mu(t)=1.0 x mu(t-1) + omega(t) with variance of
* omega(t)=5. In translating this to DLM, the <A> is 1.0, and <SW> is 5.0.
*
* The measurement equation is Y(t)=1.0 x mu(t) + nu(t) with variance of
* nu(t)=100. This gives <C>=1.0 and <SV>=100, with <Y>=the series kurit
*
* The presample state is Normal with mean 130 and variance 400. This translates to
* <X0>=130 and <SX0>=400.
*
* 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 vstate which gives the
* posterior variances of the states and from svhat 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.
*
dlm(a=1.0,c=1.0,sv=100.0,sw=5.0,x0=130.0,sx0=400.0,y=kurit,$
  yhat=yhat,svhat=svhat) 1 9 xstate vstate
set postmean 1 9 = xstate(t)(1)
set postvar  1 9 = vstate(t)(1,1)
set predict  1 9 = yhat(t)(1)
set predvar  1 9 = svhat(t)(1,1)
print(picture="*.#") / predvar predict kurit postmean postvar
*
* Intervention
* This is done with a time-varying "sw" formula and a state
* mean shifter, both of which take a different value for t==10
* than for other data points.
*
dec frml[symm] swf
dec frml[vect] zf
frml swf = ||%if(t<>10,5.0,900.0)||
frml zf  = ||%if(t<>10,0.0,143.0)||
*
* Other than the 326, these numbers are guessed from the graph
*
data(unit=input) 10 15 kurit
 326 350 315 330 305 340
dlm(a=1.0,c=1.0,sv=100.0,sw=swf,x0=130.0,sx0=400.0,y=kurit,z=zf,yhat=yhat) 1 15
set predict 1 15 = yhat(t)(1)
graph(footer="Figure 2.3 KURIT examples",grid=(t==9)) 2
# kurit
# predict