在具有季节性周期的时间序列中插入缺失值答案

【问题标题】：Interpolate missing values in a time series with a seasonal cycle在具有季节性周期的时间序列中插入缺失值
【发布时间】：2011-02-11 00:12:27
【问题描述】：

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

这是一个示例，其中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)

【问题讨论】：

此代码不运行。 index 和 obs 未定义。 zoo包中的na.approx、na.spline、na.locf等na.*函数可以填写NA的值。
谢谢，粘贴了正确的块。
请显示您用于创建情节的代码并解释“不工作”的含义。
@G。 Grothendieck：这三种插值方法不起作用，因为它们仅基于时间序列中的邻居，而不考虑每日模式。

标签： r interpolation time-series

【解决方案1】：

试试这个：

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

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

为了使用 StructTS，我必须将您的 zoo 对象转换为频率为 4 的 ts 对象。您可能想再次将拟合值更改回 zoo。

【讨论】：

谢谢，几乎完全一样。但这里有些奇怪：该图显示fit 的第一个点相距甚远（相差 0.85），而 (x-fit)^2 的总和约为 0.96。但是，如果您将 x 替换为 x <- ts(rev(myzoo), f = 4) ，则匹配变得完美。知道发生了什么吗？
zoo::na.StructTS 函数更容易执行第 2-3 行：fit2 <- na.StructTS(x) 创建一个与 x 相同的序列，其 NA 通过季节性卡尔曼滤波器填充（30.66，与 fit 相同的值在这个答案中）。

【解决方案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) 非常接近，但这只是运气。使用较长的时间序列，您应该能够校准季节性校正，以便为您提供良好的预测。对信号和季节性校正的基本机制有一个很好的先验知识会有助于更好地发挥样本性能。

【讨论】：

添加了一张图片。绘图没有问题：points(na.locf(myzoo)[10], col = "blue")
@jonw -- 哦！我误解了。我认为问题在于获得 a 点。问题在于通过这种“季节性”得到很好的估计，这实际上是一个日常的季节性。我应该尝试绘图（我刚刚尝试过?points.zoo）。

【解决方案3】：

forecast::na.interp 是一个好方法。来自documentation

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

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

This paper 针对实时序列评估几种插值方法，发现na.interp 既准确又高效：

在本文测试的 R 实现中，forecast 包中的 na.interp 和 zoo 包中的 na.StructTS 显示出最佳的整体结果。

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

另外值得注意的是，Rob Hyndman 编写了 forecast 包，并在提供了对这个问题的回答后包含了 na.interp。 na.interp 很可能是对这种方法的改进，即使它在这种情况下表现更差（可能是由于在 StructTS 中指定了时间段，而 na.interp 可以计算出来）。

【讨论】：

【解决方案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 onject 转换回 ts，因为 arima() 只接受 ts。

【讨论】：