Search code examples
rkalman-filterestimationstate-space

How to simulate the posterior filtered estimates of a Kalman Filter using the DSE package in R


How do I call for the posterior (refined) state estimates from a Kalman Filter simulation in R using the DSE package?

I have added an example below. Assume that I have created a simple random walk state space with the error being a standard normal distribution. The model is created using the SS function, with initialised state and covariance estimates of zero. The theoretical model form is thus: X(t) = X(t-1) + e(t)~N(0,1) for state evolution Y(t) = X(t) + w(t)~N(0,1)

We now implement this in R by following the instructions on page 6 and 7 of the "Kalman Filtering in R" article in the Journal of Statistical Software. First we create the state space model using the SS() function and store it in the variable called kalman.filter:

kalman.filter=dse::SS(F = matrix(1,1,1), 
                  Q = matrix(1,1,1),
                  H = matrix(1,1,1),
                  R = matrix(1,1,1),
                  z0 = matrix(0,1,1),
                  P0 = matrix(0,1,1)
                  )

Then we simulate a 100 observations from the model form using simulate() and put them in a variable called simulate.kalman.filter:

simulate.kalman.filter=simulate(kalman.filter, start = 1, freq = 1, sampleT = 100)

Then we run the kalman filter against the measurements using l() and store it under the variable called test:

test=l(kalman.filter, simulate.kalman.filter)

From the outputs, which ones are my filtered estimates?


Solution

  • I have found the answer to this question.

    Firstly, the filtered estimates of the model are not given in the l() function. This function only gives the one step ahead predictions. The above framing of my problem was coded as:

    kalman.filter=dse::SS(F = matrix(1,1,1), 
                  Q = matrix(1,1,1),
                  H = matrix(1,1,1),
                  R = matrix(1,1,1),
                  z0 = matrix(0,1,1),
                  P0 = matrix(0,1,1)
                  )
    simulate.kalman.filter=simulate(kalman.filter, start = 1, freq = 1, sampleT = 100)
    test=l(kalman.filter, simulate.kalman.filter)
    

    The one step ahead predictions are given by:

    predictions = test$estimates$pred
    

    A quick way to visualize this is given by:

    tfplot(test)
    

    This allows you to quickly plot your one step ahead predictions against the actual data. To get your filtered estimates, you need to use the smoother() function, in the same dse package. It inputs the state model as well as the data, in this case it is kalman.filter and simulate.kalman.filter respectively. The output is smoothed estimates for all the time points. But note that it does this after considering the full data set, so it does not do this as each observation comes in. See code below. The first line of the code gives you your smoothed estimates, the following lines plot the example:

    smooth = smoother(test, simulate.kalman.filter)
    plot(test$estimates$pred, ylim=c(max(test$estimates$pred,smooth$filter$track,simulate.kalman.filter$outpu), min(test$estimates$pred,smooth$filter$track,simulate.kalman.filter$output)))
    points(smooth$smooth$state, col = 3)
    points(simulate.kalman.filter$output, col = 4)
    

    The above plot plots all your actual data, model estimates and smoothed estimates against one another.