使用季节性周期对时间序列中的缺失值进行插值

Question

我有一个时间序列，我想智能地插入缺失值。特定时间的价值受到多日趋势及其在每日周期中的位置的影响。

这是一个示例，其中 myzoo 中缺少第十个观察值。

start <- as.POSIXct("2010-01-01") 
freq <- as.difftime(6, units = "hours") 
dayvals <- (1:4)*10 
timevals <- c(3, 1, 2, 4) 
index <- seq(from = start, by = freq, length.out = 16)
obs <- (rep(dayvals, each = 4) + rep(timevals, times = 4))
myzoo <- zoo(obs, index)
myzoo[10] <- NA

如果我必须实现此功能，我会使用附近日期的关闭时间的某种加权平均值，或者添加该天的值到适合更大趋势的函数线，但我希望已经存在一些适用于这种情况的包或函数？

编辑：稍微修改代码以澄清我的问题。有一些 na.* 方法可以从最近的邻居进行插值，但在这种情况下，它们无法识别缺失值是在当天最低值的时间。也许解决方案是将数据重塑为宽格式，然后进行插值，但我不想完全忽略同一天的连续值。值得注意的是，diff(myzoo, lag = 4) 返回一个 10 的向量。解决方案可能在于 reshape、na.spline 和 diff.inv 的某种组合，但我就是无法弄清楚。

以下是三种不起作用的方法：

编辑2。使用以下代码生成的图像。

myzoo <- zoo(obs, index)
myzoo[10] <- NA # knock out the missing point
plot(myzoo, type="o", pch=16) # plot solid line
points(na.approx(myzoo)[10], col = "red")
points(na.locf(myzoo)[10], col = "blue")
points(na.spline(myzoo)[10], col = "green")
myzoo[10] <- 31 # replace the missing point
lines(myzoo, type = "o", lty=3, pch=16) # dashed line over the gap
legend(x = "topleft", 
       legend = c("na.spline", "na.locf", "na.approx"), 
       col=c("green","blue","red"), pch = 1)

Answer 1

试试这个：

x <- ts(myzoo,f=4)
fit <- ts(rowSums(tsSmooth(StructTS(x))[,-2]))
tsp(fit) <- tsp(x)
plot(x)
lines(fit,col=2)

这个想法是使用时间序列的基本结构模型，该模型使用卡尔曼滤波器很好地处理缺失值。然后使用卡尔曼平滑来估计时间序列中的每个点，包括任何遗漏的点。

我必须将你的 Zoo 对象转换为频率为 4 的 ts 对象才能使用 StructTS。您可能想再次将拟合值更改回动物园。

Answer 2

在这种情况下，我认为您需要在 ARIMA 模型中进行季节性修正。这里没有足够的日期来适应季节性模型，但这应该可以帮助您开始。

library(zoo)
start <- as.POSIXct("2010-01-01") 
freq <- as.difftime(6, units = "hours") 
dayvals <- (1:4)*10 
timevals <- c(3, 1, 2, 4) 
index <- seq(from = start, by = freq, length.out = 16)
obs <- (rep(dayvals, each = 4) + rep(timevals, times = 4))
myzoo <- myzoo.orig <- zoo(obs, index)
myzoo[10] <- NA

myzoo.fixed <- na.locf(myzoo)

myarima.resid <- arima(myzoo.fixed, order = c(3, 0, 3), seasonal = list(order = c(0, 0, 0), period = 4))$residuals
myzoo.reallyfixed <- myzoo.fixed
myzoo.reallyfixed[10] <- myzoo.fixed[10] + myarima.resid[10]

plot(myzoo.reallyfixed)
points(myzoo.orig)

在我的测试中，ARMA(3, 3) 非常接近，但这只是运气。对于较长的时间序列，您应该能够校准季节性校正，以便为您提供良好的预测。事先了解信号和季节性校正的潜在机制将有助于获得更好的样本性能。

Answer 3

forecast::na.interp 是一个很好的方法。来自文档

对非季节性序列使用线性插值，并对季节性序列使用周期性 stl 分解来替换缺失值。

library(forecast)
fit <- na.interp(myzoo)
fit[10]  # 32.5, vs. 31.0 actual and 32.0 from Rob Hyndman's answer

本文针对实时序列评估了几种插值方法，并发现na.interp< /code> 既准确又高效：

从本文测试的 R 实现来看，预测包中的 na.interp 和 Zoo 包中的 na.StructTS 显示了最佳的总体结果。

na.interp 函数也不比
na.approx[最快的方法]，所以黄土分解对计算时间的要求似乎不是很高。

另外值得注意的是，Rob Hyndman 编写了 forecast 包，并在提供此问题的答案后包含了 na.interp。 na.interp 很可能是对此方法的改进，尽管它在本例中表现较差（可能是由于在 StructTS 中指定了周期，其中 na .interp 计算出来）。

Answer 4

imputeTS 包有一种对 ARIMA 模型的状态空间表示进行卡尔曼平滑的方法 - 这可能是解决此问题的一个好方法。

library(imputeTS)
na_kalman(myzoo, model = "auto.arima")

也可以直接与动物园时间序列对象一起使用。您还可以在此函数中使用您自己的 ARIMA 模型。如果您认为您可以做得更好，那么“auto.arima”。可以这样完成：

library(imputeTS)
usermodel <- arima(myts, order = c(1, 0, 1))$model
na_kalman(myts, model = usermodel)

但在这种情况下，您必须将 Zoo 对象转换回 ts，因为 arima() 只接受 ts。