我知道这个问题已经存在多年了，但是，以下想法可能有用：

1. 如果变量之间存在联系（并且理论公式效果不佳），则可以使用 PCA 以系统的方式寻找（线性）依赖关系。我将证明这对于这个问题中的给定数据很有效。

2. 鉴于没有太多数据（总共 112 个数字），只能估计几个模型参数（例如，无法拟合完整的季节性效应），尝试自定义模型可能是有意义的。

以下是我将如何根据以下原则进行预测：

第 1 步：我们可以使用 PCA 来揭示数据中的依赖关系。使用 R，数据存储在x：

> library(jvcoords)
> m <- PCA(x)
> m
PCA: mapping p = 4 coordinates to q = 4 coordinates

                              PC1         PC2          PC3          PC4
standard deviation     0.18609759 0.079351671 0.0305622047 0.0155353709
variance               0.03463231 0.006296688 0.0009340484 0.0002413477
cum. variance fraction 0.82253436 0.972083769 0.9942678731 1.0000000000

这表明前两个主成分解释了 97% 的方差，使用三个主成分覆盖了 99.4% 的方差。因此，为前两台或三台 PC 制作模型就足够了。（数据大约满足。）$W = 0.234\, wd - 1.152\, wc - 8.842 \,p$

进行 PCA 涉及找到一个正交矩阵。这种矩阵的空间是 6 维的，所以我们估计了 6 个参数。（由于我们下面只真正使用PC1，这可能是“有效”参数较少。）$4\times 4$

Step 2. PC1有明显趋势：

> t <- 1:28
> plot(m$y[,1], type = "b", ylab = "PC1")
> trend <- lm(m$y[,1] ~ t)
> abline(trend)

我创建了 PC 分数的副本，删除了这个趋势：

> y2 <- m$y
> y2[,1] <- y2[,1] - fitted(trend)

绘制其他 PC 的分数并没有显示出明显的趋势，所以我保持不变。

由于 PC 分数居中，因此趋势通过 PC1 样本的质心，拟合趋势仅对应于估计一个参数。

第 3 步。一对散点图显示没有清晰的结构，因此我将 PC 建模为独立的：

> pairs(y2, asp = 1, oma = c(1.7, 1.7, 1.7, 1.7))

第 4 步。 PC1 有明显的周期性，滞后 13（如问题所示）。这可以通过不同的方式看到。例如，滞后 13 自相关在相关图中显示为与 0 显着不同：

> acf(y2[,1])

（将数据与移动副本一起绘制时，周期性在视觉上更加引人注目。）

由于我们希望将估计参数的数量保持在较低水平，并且由于相关图显示滞后 13 是唯一具有显着贡献的滞后，因此我将 PC1 建模为，其中是独立且标准正态分布的（即这是一个大多数系数固定为 0 的 AR(13) 过程）。估计和的一种简单方法是使用以下函数：$y^{(1)}_{t+13} = \alpha_{13} y^{(1)}_t + \sigma \varepsilon_{t+13}$$\varepsilon_t$$\alpha_{13}$$\sigma$lm()

> lag13 <- lm(y2[14:28,1] ~ y2[1:15,1] + 0)
> lag13

Call:
lm(formula = y2[14:28, 1] ~ y2[1:15, 1] + 0)

Coefficients:
y2[1:15, 1]  
     0.6479  

> a13 <- coef(lag13)
> s13 <- summary(lag13)$sigma

作为合理性测试，我绘制了给定的数据（黑色），以及我的 PC1 模型的随机轨迹（蓝色），范围为未来一年：

t.f <- 29:41
pc1 <- m$y[,1]
pc1.f <- (predict(trend, newdata = data.frame(t = t.f))
          + a13 * y2[16:28, 1]
          + rnorm(13, sd = s13))
plot(t, pc1, xlim = range(t, t.f), ylim = range(pc1, pc1.f),
     type = "b", ylab = "PC1")
points(t.f, pc1.f, col = "blue", type = "b")

蓝色的模拟路径看起来像是数据的合理延续。PC2 和 PC3 的相关图没有显示出显着的相关性，因此我将这些分量建模为白噪声。PC4 确实显示了相关性，但对总方差的贡献很小，以至于它似乎不值得建模，我还将这个组件建模为白噪声。

在这里，我们又拟合了两个参数。这给我们带来了模型中总共 9 个参数（包括 PCA），当我们开始使用由 112 个数字组成的数据时，这似乎并不荒谬。

预报。 我们可以通过排除噪声（以获得平均值）并反转 PCA 来获得数字预测：

> pc1.f <- predict(trend, newdata = data.frame(t = t.f)) + a13 * y2[16:28, 1]
> y.f <- data.frame(PC1 = pc1.f, PC2 = 0, PC3 = 0, PC4 = 0)
> x.f <- fromCoords(m, y.f)
> rownames(x.f) <- t.f
> x.f
          W       wd       wc         p
29 4.456825 4.582231 3.919151 0.5616497
30 4.407551 4.563510 3.899012 0.5582053
31 4.427701 4.571166 3.907248 0.5596139
32 4.466062 4.585740 3.922927 0.5622955
33 4.327391 4.533055 3.866250 0.5526018
34 4.304330 4.524294 3.856824 0.5509898
35 4.342835 4.538923 3.872562 0.5536814
36 4.297404 4.521663 3.853993 0.5505056
37 4.281638 4.515673 3.847549 0.5494035
38 4.186515 4.479533 3.808671 0.5427540
39 4.377147 4.551959 3.886586 0.5560799
40 4.257569 4.506528 3.837712 0.5477210
41 4.289875 4.518802 3.850916 0.5499793

不确定带可以通过解析或简单地使用 Monte Carlo 获得：

N <- 1000 # number of Monte Carlo samples
W.f <- matrix(NA, N, 13)
for (i in 1:N) {
    y.f <- data.frame(PC1 = (predict(trend, newdata = 
                 data.frame(t = t.f))
              + a13 * y2[16:28, 1]
              + rnorm(13, sd = s13)),
              PC2 = rnorm(13, sd = sd(y2[,2])),
              PC3 = rnorm(13, sd = sd(y2[, 3])),
              PC4 = rnorm(13, sd = sd(y2[, 4])))
    x.f <- fromCoords(m, y.f)
    W.f[i,] <- x.f[, 1]
}
bands <- apply(W.f, 2,
               function(x) quantile(x, c(0.025, 0.15, 0.5, 
                        0.85, 0.975)))
plot(t, x$W, xlim = range(t, t.f), ylim = range(x$W, bands),
     type = "b", ylab = "W")
for (b in 1:5) {
    lines(c(28, t.f), c(x$W[28], bands[b,]), col = "grey")
}

该图显示了的实际数据，以及使用拟合模型进行预测的 60%（内三行）和 95%（外两行）不确定性带。$W$