如何使用 auto.arima 估算缺失值

机器算法验证 有马
2022-01-27 07:53:13

我有一个包含许多缺失值的动物园系列。我读到auto.arima可以估算这些缺失值?谁能教我怎么做?多谢!

这是我尝试过的,但没有成功:

fit <- auto.arima(tsx)
plot(forecast(fit))
2个回答

首先,请注意forecast计算样本外预测,但您对样本内观察感兴趣。

卡尔曼滤波器处理缺失值。因此,您可以从forecast::auto.arimaor返回的输出中获取 ARIMA 模型的状态空间形式,stats::arima并将其传递给KalmanRun.

编辑 (根据 stats0007 的回答修复代码)

在以前的版本中,我采用了与观察到的序列相关的过滤状态列,但是我应该使用整个矩阵并对观察方程进行相应的矩阵运算(感谢@stats0007 的评论。)下面我更新代码并相应地绘制。yt=Zαt

我使用一个ts对象作为示例系列而不是zoo,但它应该是相同的:

require(forecast)
# sample series
x0 <- x <- log(AirPassengers)
y <- x
# set some missing values
x[c(10,60:71,100,130)] <- NA
# fit model
fit <- auto.arima(x)
# Kalman filter
kr <- KalmanRun(x, fit$model)
# impute missing values Z %*% alpha at each missing observation
id.na <- which(is.na(x))
for (i in id.na)
  y[i] <- fit$model$Z %*% kr$states[i,]
# alternative to the explicit loop above
sapply(id.na, FUN = function(x, Z, alpha) Z %*% alpha[x,], 
  Z = fit$model$Z, alpha = kr$states)
y[id.na]
# [1] 4.767653 5.348100 5.364654 5.397167 5.523751 5.478211 5.482107 5.593442
# [9] 5.666549 5.701984 5.569021 5.463723 5.339286 5.855145 6.005067

您可以绘制结果(对于整个系列和整个年份,样本中间缺少观测值):

par(mfrow = c(2, 1), mar = c(2.2,2.2,2,2))
plot(x0, col = "gray")
lines(x)
points(time(x0)[id.na], x0[id.na], col = "blue", pch = 19)
points(time(y)[id.na], y[id.na], col = "red", pch = 17)
legend("topleft", legend = c("true values", "imputed values"), 
  col = c("blue", "red"), pch = c(19, 17))
plot(time(x0)[60:71], x0[60:71], type = "b", col = "blue", 
  pch = 19, ylim = range(x0[60:71]))
points(time(y)[60:71], y[60:71], col = "red", pch = 17)
lines(time(y)[60:71], y[60:71], col = "red")
legend("topleft", legend = c("true values", "imputed values"), 
  col = c("blue", "red"), pch = c(19, 17), lty = c(1, 1))

原始系列的图和归因于缺失观测值的值

您可以使用卡尔曼平滑器而不是卡尔曼滤波器重复相同的示例。所有你需要改变的是这些行:

kr <- KalmanSmooth(x, fit$model)
y[i] <- kr$smooth[i,]

通过卡尔曼滤波器处理缺失的观测值有时被解释为序列的外推;当使用卡尔曼平滑器时,缺失的观察被认为是通过观察序列中的插值来填充的。

这将是我的解决方案:

# Take AirPassengers as example
data <- AirPassengers

# Set missing values
data[c(44,45,88,90,111,122,129,130,135,136)] <- NA


missindx <- is.na(data)

arimaModel <- auto.arima(data)
model <- arimaModel$model

#Kalman smoothing
kal <- KalmanSmooth(data, model, nit )
erg <- kal$smooth  

for ( i in 1:length(model$Z)) {
       erg[,i] = erg[,i] * model$Z[i]
}
karima <-rowSums(erg)

for (i in 1:length(data)) {
  if (is.na(data[i])) {
    data[i] <- karima[i]
  }
}
#Original TimeSeries with imputed values
print(data)

@贾夫拉卡勒:

谢谢你的帖子,很有趣!

我对您的解决方案有两个问题,希望您能帮助我:

  1. 你为什么使用 KalmanRun 而不是 KalmanSmooth ?我读到 KalmanRun 被认为是外推,而平滑是估计。

  2. 我也没有得到你的身份证部分。为什么不使用 .Z 中的所有组件?我的意思是,例如 .Z 给出 1, 0,0,0,0,1,-1 -> 7 个值。这意味着 .smooth (在你的情况下是 KalmanRun 状态)给了我 7 列。据我了解,所有 1 或 -1 的列都进入模型。

    假设 AirPass 中缺少第 5 行。然后我会像这样取第 5 行的总和:我会从第 1 列添加值(因为 Z 给了 1),我不会添加第 2-4 列(因为 Z 说 0),我会添加第 5 列,我会添加第 7 列的负值(因为 Z 表示 -1)

    我的解决方案错了吗?或者他们俩都好吗?您能否进一步向我解释一下?