如何使用“KFAS”R 包和 AR(1) 转移方程估计卡尔曼滤波器？

Question

我正在使用 R 中的“KFAS”包来估计带有卡尔曼滤波器的状态空间模型。我的测量和转移方程是：

y_t = Z_t * x_t + ps_t（测量）

x_t = T_t * x_{t-1} + R_t * ta_t（转变），

ps_t ~ N(0,H_t) 和 ta_t ~ N(0,Q_t)。

所以，我想估计方差 H_t 和 Q_t，还有 T_t，即 AR(1) 系数。我的代码如下：

library(KFAS)

set.seed(100)

eps <- rt(200, 4, 1)
meas <- as.matrix((arima.sim(n=200, list(ar=0.6), innov = rnorm(200)*sqrt(0.5)) + eps), 
ncol=1)

Zt <- 1
Ht <- matrix(NA)
Tt <- matrix(NA)
Rt <- 1
Qt <- matrix(NA)

ss_model <- SSModel(meas ~ -1 + SSMcustom(Z = Zt, T = Tt, R = Rt, 
                                             Q = Qt), H = Ht)
fit <- fitSSM(ss_model, inits = c(0,0.6,0), method = 'L-BFGS-B')

但它返回：“is.SSModel(do.call(updatefn, args = c(list(inits, model), update_args)))中的错误：系统矩阵（不包括Z）包含NA或无限值，协方差矩阵包含值大于 1e+07"

如软件包论文中所述，方差的 NA 定义效果很好。然而，对于 AR 系数似乎不能这样做。有谁知道我该怎么做？

请注意，我知道 SSMarima 函数，它简化了 ARIMA 模型中转移方程的定义。虽然我能够估计 AR(1) 系数。和Q_t这样，我仍然无法估计ps_t方差（H_t）。此外，我正在将我的卡尔曼滤波器代码从EViews迁移到R，因此我需要学习SSMcustom以用于其他更复杂的模型。

谢谢！

Answer 1

您的示例中似乎遗漏了一些内容，因为您的错误消息来自函数

fitSSM

。如果你想使用

fitSSM

来估计一般状态空间模型，你需要提供自己的模型更新函数。默认行为只能处理协方差矩阵 H 和 Q 中的 NA。

fitSSM

的主要目标只是从简单的东西开始。对于复杂模型和/或大数据，我建议手动使用您自己编写的目标函数（在

logLik

方法的帮助下）和您最喜欢的数值优化例程，以获得最大性能。像这样的东西：

library(KFAS)

set.seed(100)

eps <- rt(200, 4, 1)
meas <- as.matrix((arima.sim(n=200, list(ar=0.6), innov = rnorm(200)*sqrt(0.5)) + eps), 
ncol=1)

Zt <- 1
Ht <- matrix(NA)
Tt <- matrix(NA)
Rt <- 1
Qt <- matrix(NA)

ss_model <- SSModel(meas ~ -1 + SSMcustom(Z = Zt, T = Tt, R = Rt, 
                                             Q = Qt), H = Ht)

objf <- function(pars, model, estimate = TRUE) {
  model$H[1] <- pars[1]
  model$T[1] <- pars[2]
  model$Q[1] <- pars[3]
  if (estimate) {
  -logLik(model)
  } else {
    model
  }
}

opt <- optim(c(1, 0.5, 1), objf, method = "L-BFGS-B", 
  lower = c(0, -0.99, 0), upper = c(100, 0.99, 100), model = ss_model)

ss_model_opt <- objf(opt$par, ss_model, estimate = FALSE)

与

fitSSM

相同：

updatefn <- function(pars, model) {
  model$H[1] <- pars[1]
  model$T[1] <- pars[2]
  model$Q[1] <- pars[3]
  model
}


fit <- fitSSM(ss_model, c(1, 0.5, 1), updatefn, method = "L-BFGS-B", 
  lower = c(0, -0.99, 0), upper = c(100, 0.99, 100))

identical(ss_model_opt, fit$model)

Answer 2

我正在使用最新的 KFAS 包，它是 SSMarima 方法来构建 AR(1)。我按照此https://www.jstatsoft.org/article/view/v078i10中给出的说明进行操作。下面是我的代码片段。

model_arima <- SSModel(ssm_total_vol$total_vol ~ SSMarima(ar = 1, d = 0, Q = 1, 
                                                          stationary = FALSE))

update_model <- function(pars, model) {
  tmp <- SSMarima(ar = pars[1], d = 0, Q = pars[2])
  model["R", states = "arima"] <- tmp$R
  model["Q", states = "arima"] <- tmp$Q
  model["P1", states = "arima"] <- tmp$P1
  model
  }

fit_arima <- fitSSM(model_arima, inits = c(0, 1), updatefn = update_model,
                    method = "L-BFGS-B", lower = c(-1, 0), upper = c(1, 100))

但是，运行 fitSSM 方法后，我收到错误“is.SSModel(do.call(updatefn, args = c(list(inits, model), update_args)), 中的错误： 矩阵 P1 包含用于漫反射状态的非零对角线元素。”我的数据集是非平稳的每周数据集。请提供帮助吗？