(Fast) Stationary Kalman Filter and Smoother
SKFS(X, A, C, Q, R, F_0, P_0, loglik = FALSE)
numeric data matrix (\(T \times n\)).
transition matrix (\(rp \times rp\)).
observation matrix (\(n \times rp\)).
state covariance (\(rp \times rp\)).
observation covariance (\(n \times n\)).
initial state vector (\(rp \times 1\)).
initial state covariance (\(rp \times rp\)).
logical. Compute log-likelihood?
All results from SKF
and FIS
, and additionally
a \(rp \times rp \times T\) matrix PPm_smooth
, which is equal to the estimate of \(Cov(F^smooth_t, F^smooth_{t-1} | T)\) and needed for EM iterations.
See 'Property 6.3: The Lag-One Covariance Smoother' in Shumway & Stoffer (2017).
Shumway, R. H., & Stoffer, D. S. (2017). Time Series Analysis and Its Applications: With R Examples. Springer.
library(collapse)
#> collapse 1.9.5, see ?`collapse-package` or ?`collapse-documentation`
#>
#> Attaching package: ‘collapse’
#> The following object is masked from ‘package:vars’:
#>
#> B
#> The following object is masked from ‘package:stats’:
#>
#> D
## Two-Step factor estimates from monthly BM (2014) data
X <- fscale(diff(qM(BM14_M))) # Standardizing as KF has no intercept
r <- 5L # 5 Factors
p <- 3L # 3 Lags
n <- ncol(X)
## Initializing the Kalman Filter with PCA results
X_imp <- tsnarmimp(X) # Imputing Data
v <- eigen(cov(X_imp))$vectors[, 1:r] # PCA
F_pc <- X_imp %*% v # Principal component factor estimates
C <- cbind(v, matrix(0, n, r*p-r)) # Observation matrix
res <- X - tcrossprod(F_pc, v) # Residuals from static predictions
R <- diag(fvar(res)) # Observation residual covariance
var <- .VAR(F_pc, p) # VAR(p)
A <- rbind(t(var$A), diag(1, r*p-r, r*p))
Q <- matrix(0, r*p, r*p) # VAR residual matrix
Q[1:r, 1:r] <- cov(var$res)
F_0 <- var$X[1L, ] # Initial factor estimate and covariance
P_0 <- ainv(diag((r*p)^2) - kronecker(A,A)) %*% unattrib(Q)
dim(P_0) <- c(r*p, r*p)
## Run standartized data through Kalman Filter and Smoother once
kfs_res <- SKFS(X, A, C, Q, R, F_0, P_0, FALSE)
## Two-step solution is state mean from the Kalman Smoother
F_kal <- kfs_res$F_smooth[, 1:r, drop = FALSE]
colnames(F_kal) <- paste0("f", 1:r)
## See that this is equal to the Two-Step estimate by DFM()
all.equal(F_kal, DFM(X, r, p, em.method = "none", pos.corr = FALSE)$F_2s)
#> [1] TRUE
## Same in two steps using SKF() and FIS()
kfs_res2 <- with(SKF(X, A, C, Q, R, F_0, P_0, FALSE), FIS(A, F, F_pred, P, P_pred))
F_kal2 <- kfs_res2$F_smooth[, 1:r, drop = FALSE]
colnames(F_kal2) <- paste0("f", 1:r)
all.equal(F_kal, F_kal2)
#> [1] TRUE
rm(X, r, p, n, X_imp, v, F_pc, C, res, R, var, A, Q, F_0, P_0, kfs_res, F_kal, kfs_res2, F_kal2)