R/kalman.states.R

In JannesRed/rStar: Natural Rate of Interest

#------------------------------------------------------------------------------#
# File:        kalman.states.R
#
# Description: The function kalman.states() calls the functions
#              kalman.states.filtered() and kalman.states.smoothed() to
#              apply the Kalman filter and smoother.
#              It takes as input the coefficient matrices for the given
#              state-space model, with notation matching Hamilton (1994),
#              as well as conditional expectation and covariance matrix
#              of the initial state, xi.tm1tm1 and P.tm1tm1 respectively.
#------------------------------------------------------------------------------#
kalman.states <- function(xi.tm1tm1, P.tm1tm1, F_matrix, Q, A, H, R, y, x) {
  filtered <- kalman.states.filtered(xi.tm1tm1, P.tm1tm1, F_matrix, Q, A, H, R, y, x)
  smoothed <- kalman.states.smoothed(filtered$xi.ttm1, filtered$P.ttm1, filtered$xi.tt, filtered$P.tt,
                                     F_matrix, Q, A, H, R, y, x)
  return(list("filtered"=filtered, "smoothed"=smoothed))
}
kalman.states.filtered <- function(xi.tm1tm1, P.tm1tm1, F_matrix, Q, A, H, R, y, x, t=1) {
  xi.ttm1 <- as.vector(F_matrix %*% xi.tm1tm1)
  P.ttm1 <- F_matrix %*% P.tm1tm1 %*% t(F_matrix) + Q
  prediction.error <- (as.vector(y[t,]) - as.vector(t(A) %*% as.vector(x[t,])) - as.vector(t(H) %*% xi.ttm1))
  HPHR <- t(H) %*% P.ttm1 %*% H + R
  xi.tt <- xi.ttm1 + as.vector(P.ttm1 %*% H %*% solve(HPHR, prediction.error))
  P.tt <- P.ttm1 - P.ttm1 %*% H %*% solve(HPHR, t(H) %*% P.ttm1)
  if (t == dim(y)[1]) {
      return(list("xi.ttm1"=xi.ttm1, "P.ttm1"=P.ttm1, "xi.tt"=xi.tt, "P.tt"=P.tt))
  } else {
      tmp <- kalman.states.filtered(xi.tt, P.tt, F_matrix, Q, A, H, R, y, x, t+1)
      return(list("xi.ttm1"=rbind(xi.ttm1, tmp$xi.ttm1),
                  "P.ttm1"=rbind(P.ttm1, tmp$P.ttm1),
                  "xi.tt"=rbind(xi.tt, tmp$xi.tt),
                  "P.tt"=rbind(P.tt, tmp$P.tt)))
  }
}
kalman.states.smoothed <- function(xi.ttm1.array, P.ttm1.array, xi.tt.array, P.tt.array,
                                   F_matrix, Q, A, H, R, y, x, t=dim(y)[1], xi.tp1T=NA, P.tp1T=NA) {
  n <- dim(xi.ttm1.array)[2]
  if (t == dim(y)[1]) {
    xi.tT <- xi.tt.array[t,]
    P.tT <- P.tt.array[((t-1)*n+1):(t*n),]
    tmp <- kalman.states.smoothed(xi.ttm1.array, P.ttm1.array, xi.tt.array, P.tt.array,
                                  F_matrix, Q, A, H, R, y, x, t-1, xi.tT, P.tT)
    return(list("xi.tT"=rbind(tmp$xi.tT, xi.tT),
                "P.tT" =rbind(tmp$P.tT, P.tT)))
  } else {
    P.tt <- P.tt.array[((t-1)*n+1):(t*n),]
    P.tp1t <- P.ttm1.array[(t*n+1):((t+1)*n),]

    #added for debug
    # P.tp1t matrix is singular. And the source of my pain
    # #-------------------------------------------------------------------------
    # if (!matrixcalc::is.singular(P.tp1t)) {
    #   P.tp1t <- MASS
    #   print("P.tp1t matrix singular. Found closest PD")
    # }

    J.t <- P.tt %*% t(F_matrix) %*% solve(P.tp1t)
    #-------------------------------------------------------------------------
    xi.tt <- xi.tt.array[t,]
    xi.tp1t <- xi.ttm1.array[t+1,]
    xi.tT <- xi.tt + as.vector(J.t %*% (xi.tp1T - xi.tp1t))
    P.tT <- P.tt + J.t %*% (P.tp1T - P.tp1t) %*% t(J.t)
    if (t > 1) {
      tmp <- kalman.states.smoothed(xi.ttm1.array, P.ttm1.array, xi.tt.array, P.tt.array,
                                    F_matrix, Q, A, H, R, y, x, t-1, xi.tT, P.tT)
      return(list("xi.tT"=rbind(tmp$xi.tT, xi.tT),
                  "P.tT" =rbind(tmp$P.tT, P.tT)))
    } else {
      return(list("xi.tT"=xi.tT, "P.tT"=P.tT))
    }
  }
}