在具有季节性循环的时间序列中插值缺失值

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)

- J. Win.

---

这段代码无法运行。index和obs未定义。zoo包中的na.approx、na.spline、na.locf和其他na.*函数可以填充NA值。 - G. Grothendieck

谢谢，已经粘贴了正确的代码块。 - J. Win.

请展示您用于创建图表的代码，并解释“不起作用”是什么意思。 - G. Grothendieck

@G. Grothendieck：这三种插值方法不可行，因为它们仅基于时间序列中的邻居数据，而没有考虑每日模式。 - J. Win.

Answer 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格式。

- Rob Hyndman

---

谢谢，这几乎完美解决了问题。但是有一些奇怪的事情：图表显示fit的第一个点相差很远（0.85），而且（x-fit）^2的总和约为0.96。但是，如果您将x替换为x <- ts(rev(myzoo), f = 4)，拟合就变得完美了。您知道发生了什么吗？ - J. Win.

zoo::na.StructTS函数可以更轻松地执行第2-3行：fit2 <- na.StructTS(x)通过季节性卡尔曼滤波器（30.66，与此答案中的fit相同）填充NA，创建与x相同的系列。 - Max Ghenis

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)非常接近，但那只是运气。对于更长的时间序列，您应该能够校准季节性修正以获得良好的预测结果。有一个关于信号和季节性修正的基本机制的良好先验将有助于提高样本外表现。

- Richard Herron

---

添加了一张图片。绘图很容易：points(na.locf(myzoo)[10], col = "blue")。 - J. Win.

@jonw -- 哦！我误解了。我以为问题是获取一个点。问题是如何获得这个“季节性”的良好估计值，这实际上是每日季节性。我应该尝试绘图（我刚刚尝试了?points.zoo）。 - Richard Herron

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既准确又高效：

在本文测试的R实现中，来自预测包的na.interp和来自zoo包的na.StructTS表现最佳。

na.interp函数的速度也不比最快的na.approx慢多少，因此loess分解在计算时间方面似乎并不是很苛刻。

值得注意的是，Rob Hyndman编写了forecast包，并在回答这个问题后包含了na.interp。尽管它在这个实例中表现更差（可能是由于在StructTS中指定周期，而na.interp则自动计算），但na.interp很可能是这种方法的改进。

- Max Ghenis

Answer 4

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

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

这个函数也可以直接处理zoo时间序列对象。您还可以在此功能中使用自己的ARIMA模型。如果您认为您能做得更好，超过“auto.arima”，可以采用以下方式完成：

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

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

- Steffen Moritz