您的目的是对时间序列执行固定间隔(FI) 平滑。观察时间的平滑值$t$被定义为条件期望

$$\widehat{Y}_{t} := \mathbb{E}(Y_t|\mathbf{Y}_{1:r},\,\mathbf{Y}_{s:n})$$ 符号在哪里$\mathbf{Y}_{u:v} := [Y_u,\,Y_{u+1}, \, \dots,\,Y_v]$是从时间观察的向量$u$到时间$v$. 上面，假设间隙是从时间$r+1$到$s-1$， 和$n$ 是整个系列的长度。时间 $t$在差距中，期望可以写成$\widehat{Y}_{t|1:r, s:n}$ 回忆它的条件性质。

平滑值没有您猜测的简单形式。对于已知协方差结构的高斯平稳时间序列，估计的 $\widehat{Y}_{t}$为了$t$在间隙中可以通过求解线性系统来找到。

当时间序列模型可以置于状态空间 (SS) 形式时，FI 平滑是基于卡尔曼滤波的标准操作，并且可以使用可用的 R 函数来完成。您只需指定间隙中的值缺失即可。平滑算法估计隐藏状态 ，其中包含间隙中的所有有关的相关信息。ARIMA 模型可以采用 SS 形式。$\boldsymbol{\alpha}_t$$Y_t$$t$

有趣的是，FI 平滑可以写成两个过滤器的组合：一个前向和一个后向，导致您期望的那种公式，但是对于隐藏状态估计 （预测和反推），但不适用于观察。这称为Rauch-Tung-Striebel 过滤。$\boldsymbol{\alpha}_t$$Y_t$

至少在乘法版本中，像 Holt-Winters 这样的“临时”预测程序依赖于没有简单 FI 算法的随机模型，因为它们不能以 SS 形式表示。平滑公式可能可以通过使用 SS 模型来近似，但是使用 带有对数变换的结构时间序列模型要简单得多。R stats包的 'KalmanSmooth'、'tsSmooth' 和 'StructTS' 函数可以完成这项工作。您应该查看 R 帮助页面中引用的 Harvey 或 Durbin 和 Koopman 的书籍。提供条件方差$Y_t$并且可用于构建平滑间隔，通常在间隙中间较大。但是请注意，结构模型的估计可能很困难。

AP <- log10(AirPassengers) 
## Fit a Basic Structural Model
fit <- StructTS(AP, type = "BSM")

## Fit with a gap
AP.gap <- AP
AP.gap[73:96] <- NA
fit.gap <- StructTS(AP.gap, type = "BSM", optim.control = list(trace = TRUE))

# plot in orginal (non-logged) scale
plot(AirPassengers, col = "black", ylab = "AirPass")
AP.missing <- ts(AirPassengers[73:96], start=1955, , freq=12)
lines(AP.missing, col = "grey", lwd = 1)

## smooth and sum 'level' and 'sea' to retrieve series
sm <- tsSmooth(fit.gap)
fill <- apply(as.matrix(sm[ , c(1,3)]), 1, sum)
AP.fill <- ts(fill[73:96], start=1955, , freq=12)
lines(10^AP.fill, col = "red", lwd = 1)

我发现你建议的方法，即采用前后投射的方法，很有趣。

可能值得指出的一件事是，在任何表现出混乱结构的系统中，预测可能在较短的时期内更准确。并非所有系统都是如此，例如阻尼摆可以由具有错误周期的函数建模，在这种情况下，所有中期预测都可能是错误的，而长期预测都将是非常准确，因为系统收敛到零。但在我看来，从问题中的图表来看，这可能是一个合理的假设。

这意味着我们最好更多地依赖缺失期早期的预测数据，而更多地依赖后期的回溯数据。最简单的方法是使用线性递减的权重进行预测，反之则相反：

> n <- [number of missing datapoints] 
> w <- seq(1, 0, by = -1/(n+1))[2:(n+1)]

这在第一个元素上给出了一点权重。如果您只想使用第一个插值点的预测值，也可以使用 n-1，末尾不带下标。

> w
 [1] 0.92307692 0.84615385 0.76923077 0.69230769 0.61538462 0.53846154
 [7] 0.46153846 0.38461538 0.30769231 0.23076923 0.15384615 0.07692308

我没有你的数据，所以让我们在 R 中的 AirPassenger 数据集上试试这个。我将删除中心附近的两年期：

> APearly <- ts(AirPassengers[1:72], start=1949, freq=12)
> APlate <- ts(AirPassengers[97:144], start=1957, freq=12)
> APmissing <- ts(AirPassengers[73:96], start=1955, freq=12)
> plot(AirPassengers)
# plot the "missing data" for comparison
> lines(APmissing, col="#eeeeee")
# use the HoltWinters algorithm to predict the mean:
> APforecast <- hw(APearly)[2]$mean
> lines(APforecast, col="red")
# HoltWinters doesn't appear to do backcasting, so reverse the ts, forecast, 
# and reverse again (feel free to edit if there's a better process)
> backwards <- ts(rev(APlate), freq=12)
> backcast <- hw(backwards)[2]$mean
> APbackcast <- ts(rev(backcast), start=1955, freq=12)
> lines(APbackcast, col='blue')
# now the magic: 
> n <- 24 
> w <- seq(1, 0, by=-1/(n+1))[2:(n+1)]
> interpolation = APforecast * w + (1 - w) * APbackcast
> lines(interpolation, col='purple', lwd=2)

还有你的插值。

当然，它并不完美。我猜这是因为数据前半部分的模式与后半部分的模式不同（早些年 7 月至 8 月的峰值并不那么强烈）。但正如您从图像中看到的那样，这显然比仅预测或仅靠背投要好。我想您的数据可能会得到稍微不太可靠的结果，因为没有如此强烈的季节性变化。

我的猜测是你也可以尝试这个，包括置信区间，但我不确定这样做的有效性。

假设您有单独预测和回溯的平方预测误差，我建议您这样做：设 w 为长度为 12 的向量，设 m 为您感兴趣的月份。

w=rep(NA,12);
for(w in 1:12){
w[m]=SPE_Backcast[m]/(SPE_Backcast[m]+SPE_Forecast[m]);
}

现在 w 是预测的权重，1-w 是回溯的权重。