The following examples all consider the local linear "smooth" trend or integrated random walk model. This model is used to show how easy Kalman filtering, moment smoothing and simulation smoothing can be implemented using SsfPack. The time series consists of nine observations: 1,9,2,5,8,4,6,7,3. More details can be found in the SsfPack documentation.

kalman filter | smoothing | simulation

ssfkf.ox (back)

#include < oxstd.h >
#include < packages/ssfpack/ssfpack.h >

main()
{
   decl mphi = < 1,1;0,1;1,0 > ;
   decl momega = diag(< 0,0.1,1 >);
   decl msigma = < 0,0;0,0;1,.5 > ;
  
   decl mr = sqrt(momega) * rann(3, 21);
   decl md = SsfRecursion(mr, mphi, momega, msigma);
   decl myt = md[2][1:];    // 20 observations
  
   decl mkf = KalmanFil(myt, mphi, momega);
   print("mKF\' (t=10)", "%c", {"v", "K(1,1)", "K(2,1)", "F^-1"},
         mkf[][9]');

}

mKF' (t=10)

          v     K(1,1)     K(2,1)       F^-1

    -3.0347    0.76491    0.21161    0.44669

ssfsmo.ox (back)

#include < oxstd.h >
#include < packages/ssfpack/ssfpack.h >

main()
{
   decl mphi = < 1,1;0,1;1,0 > ;
   decl momega = diag(< 0,0.1,1 >);
   decl msigma = < 0,0;0,0;1,.5 > ;
   
   decl mr = sqrt(momega) * rann(3, 21);
   decl md = SsfRecursion(mr, mphi, momega, msigma);
   decl myt = md[2][1:];    // 20 observations

   decl mkf = KalmanFil(myt, mphi, momega);
   print("mKF\' (t=10)", "%c", {"v", "K(1,1)", "K(2,1)", "F^-1"},
         mkf[][9]');

   decl mks = KalmanSmo(mkf, mphi, momega);
   print("Basic smoother output: mKS\' (t=10)", "%c",
         {"r","e(1,1)","e(2,1)","N","D(1,1)","D(2,2)"}, mks[][10]');

   decl msmodist = mks[0:2][0] ~ momega * mks[0:2][1:];
   print("Smoothed disturbances (t=10)", "%c",
         {"E[H.eps](1,1)","E[H.eps](2,1)","E[G.eps]"}, msmodist[][10]');

   decl msmostat = SsfRecursion(msmodist, mphi, momega);
   print("Smoothed states (t=10)",	"%c",
         {"alphahat(1,1)[t+1]","alphahat(2,1)[t+1]","y[t]"}, msmostat[][10]');
}

mKF' (t=10)

          v     K(1,1)     K(2,1)       F^-1

    -3.0347    0.76491    0.21161    0.44669

Basic smoother output: mKS' (t=10)

          r     e(1,1)     e(2,1)          N     D(1,1)     D(2,2)

   -0.64966     1.3099    -1.1358    0.60208     2.0578    0.79365

Smoothed disturbances (t=10)

E[H.eps](1,1)  E[H.eps](2,1)     E[G.eps]

      0.00000        0.13099      -1.1358

Smoothed states (t=10)

alphahat(1,1)[t+1]  alphahat(2,1)[t+1]         y[t]

           0.31542            -0.27487     -0.41452

ssfsim.ox (back)

#include < oxstd.h >
#include < packages/ssfpack/ssfpack.h >
							
main()
{
   decl mphi = < 1,1;0,1;1,0 > ;
   decl momega = diag(< 0,0.1,1 >);
   decl msigma = < 0,0;0,0;1,.5 > ;
   
   decl md = SsfRecursion(sqrt(momega) * rann(3, 21), mphi, momega, msigma);
   decl myt = md[2][1:];
   decl mkf = KalmanFil(myt, mphi, momega);

   decl ct = columns(myt);    // 20 observations
   decl mgamma = diag(< 0,1,0 >);
   decl mwgt = SimSmoWgt(mgamma, mkf, mphi, momega);

   print("Simulation smoother weights (t=10)", "%c",
         {"A(1,1)","A(1,2)","A(1,3)","B(1,1)"}, mwgt[][9]');

   // draw 1
   md = SimSmoDraw(mgamma, rann(1, ct), mwgt, mkf, mphi, momega);
   print("Draw 1 for slope disturbances (t=10)", "%c",
         {"H.eps(1,1)","H.eps(2,1)","G.eps"}, md[][10]');
   md = SsfRecursion(md, mphi, momega);
   print("Draw 1 for state and signal (t=10)", "%c",
         {"alpha(1,1)[t+1]","alpha(2,1)[t+1]","y[t]"}, md[][10]');

   // draw 2
   decl md2 = SimSmoDraw(mgamma, rann(1, ct), mwgt, mkf, mphi, momega);
   md2 = SsfRecursion(md2, mphi, momega);

   // draw 3
   decl md3 = SimSmoDraw(mgamma, rann(1, ct), mwgt, mkf, mphi, momega);
   md3 = SsfRecursion(md3, mphi, momega);
}

Simulation smoother weights (t=10)

     A(1,1)     A(1,2)     A(1,3)     B(1,1)

   -0.40350     1.5248  -0.014001    0.24905

Draw 1 for slope disturbances (t=10)

   H.eps(1,1)   H.eps(2,1)        G.eps

      0.00000     -0.25514      0.00000

Draw 1 for state and signal (t=10)

alpha(1,1)[t+1]   alpha(2,1)[t+1]         y[t]

      1.1744           -0.78113         1.7004