SKF function

(Fast) Stationary Kalman Filter

(Fast) Stationary Kalman Filter

A simple and fast C++ implementation of the Kalman Filter for stationary data with time-invariant system matrices and missing data.

SKF(X, A, C, Q, R, F_0, P_0, loglik = FALSE)

Arguments

  • X: numeric data matrix (TxnT x n).
  • A: transition matrix (rpxrprp x rp).
  • C: observation matrix (nxrpn x rp).
  • Q: state covariance (rpxrprp x rp).
  • R: observation covariance (nxnn x n).
  • F_0: initial state vector (rpx1rp x 1).
  • P_0: initial state covariance (rpxrprp x rp).
  • loglik: logical. Compute log-likelihood?

Returns

Predicted and filtered state vectors and covariances. - F: TxrpT x rp filtered state vectors.

  • P: rpxrpxTrp x rp x T filtered state covariances.

  • F_pred: TxrpT x rp predicted state vectors.

  • P_pred: rpxrpxTrp x rp x T predicted state covariances.

  • loglik: value of the log likelihood.

Details

The underlying state space model is:

xt=CFt+etN(0,R)x(t)=CF(t)+e(t) N(0,R) \textbf{x}_t = \textbf{C} \textbf{F}_t + \textbf{e}_t \sim N(\textbf{0}, \textbf{R})x(t) = C F(t) + e(t) ~ N(0, R) Ft=A Ft1+utN(0,Q)F(t)=AF(t1)+u(t) N(0,Q) \textbf{F}_t = \textbf{A F}_{t-1} + \textbf{u}_t \sim N(\textbf{0}, \textbf{Q})F(t) = A F(t-1) + u(t) ~ N(0, Q)

where x(t)x(t) is X[t, ]. The filter then first performs a time update (prediction)

Ft=A Ft1F(t)=AF(t1) \textbf{F}_t = \textbf{A F}_{t-1}F(t) = A F(t-1) Pt=A Pt1A+QP(t)=AP(t1)A+Q \textbf{P}_t = \textbf{A P}_{t-1}\textbf{A}' + \textbf{Q}P(t) = A P(t-1) A' + Q

where P(t)=Cov(F(t))P(t) = Cov(F(t)). This is followed by the measurement update (filtering)

Kt=PtC(C PtC+R)1K(t)=P(t)Cinv(CP(t)C+R) \textbf{K}_t = \textbf{P}_t \textbf{C}' (\textbf{C P}_t \textbf{C}' + \textbf{R})^{-1}K(t) = P(t) C' inv(C P(t) C' + R) Ft=Ft+Kt(xtC Ft)F(t)=F(t)+K(t)(x(t)CF(t)) \textbf{F}_t = \textbf{F}_t + \textbf{K}_t (\textbf{x}_t - \textbf{C F}_t)F(t) = F(t) + K(t) (x(t) - C F(t)) Pt=PtKtC PtP(t)=P(t)K(t)CP(t) \textbf{P}_t = \textbf{P}_t - \textbf{K}_t\textbf{C P}_tP(t) = P(t) - K(t) C P(t)

If a row of the data is all missing the measurement update is skipped i.e. the prediction becomes the filtered value. The log-likelihood is computed as

1/2tlog(St)etStetnlog(2π)1/2sumt[log(det(S(t)))e(t)S(t)e(t)nlog(2pi)] 1/2 \sum_t \log(|St|)-e_t'S_te_t-n\log(2\pi)1/2 sum_t[log(det(S(t))) - e(t)' S(t) e(t) - n log(2 pi)]

where S(t)=inv(CP(t)C+R)S(t) = inv(C P(t) C' + R) and e(t)=x(t)CF(t)e(t) = x(t) - C F(t) is the prediction error.

For further details see any textbook on time series such as Shumway & Stoffer (2017), which provide an analogous R implementation in astsa::Kfilter0. For another fast (C-based) implementation that also allows time-varying system matrices and non-stationary data see FKF::fkf.

Examples

# See ?SKFS

References

Shumway, R. H., & Stoffer, D. S. (2017). Time Series Analysis and Its Applications: With R Examples. Springer.

Harvey, A. C. (1990). Forecasting, structural time series models and the Kalman filter.

Hamilton, J. D. (1994). Time Series Analysis. Princeton university press.

See Also

FIS SKFS

  • Maintainer: Sebastian Krantz
  • License: GPL-3
  • Last published: 2024-06-09