---

作者

---

可下载资源

从这张图中可以清晰地看出，乘客数量存在强烈的季节性波动，但这种模式也似乎在随时间变化。一种常见的检查可能季节性模式的方法是计算和绘制STL分解（使用Loess的季节性和趋势分解），该方法利用Loess平滑来迭代地平滑和移除数据的某些成分（即可能是非线性的长期趋势和可能随时间变化的季节性）。

	  autoplot()

可以通过修改t.window和s.window参数来调整趋势和季节性成分被平滑的周期，但这些对于本文来说并不是很重要。关键信息是，似乎存在一个长期的非线性趋势和一个随时间变化的循环季节性模式，这种变化不仅体现在幅度上，可能还在形状上略有不同。

自相关函数（ACF）图清晰地显示了强烈的自相关性和季节性证据，这正是我们希望在我们的模型中捕捉到的。

使用张量积交互项拟合mvgam模型

现在我们可以拟合一个模型来尝试捕捉这两个组成部分，即趋势和季节性模式。

尽管此模型的样本外预测表现尚可，但不确定性区间可能稍显宽泛。如果你花些时间理解样条如何超出训练数据进行外推，这并不令人意外。在实际应用中，可能需要进一步调整模型或考虑其他方法来缩小预测的不确定性区间。

---

---

dplyr::glimpse(fou_terms)

 

接下来，我们可以将这些新的预测变量添加到训练和测试数据集中。重要的是，我们也需要为样本外测试集添加这些项，因为这将用于帮助我们从模型中计算预测值。

然后，我们可以拟合一个mvgam()模型，该模型包括平滑的非线性时间趋势和时变季节性，但使用傅里叶项代替上面使用的循环样条。这首先通过使用时间的单调样条来实现，确保时间趋势永远不会下降，这在当前情况下似乎是合适的。时变季节性通过允许傅里叶项的系数随时间变化来实现，这通过使用s()包装器中的'by'参数并利用单独的样条来实现。这实际上设置了时间的独立样条，然后与傅里叶项进行交互，为每个傅里叶项制定了一个时变系数模型。

mod2 <- mvgam(y ~) 
                # 时间单调平滑，以确保趋势
                # 要么增加要么变平，但不会 
                # 下降
                s(time, bs = 'moi', k = 10) +
                
                # 时变傅立叶系数，以捕捉
                # 可能变化的季节性
                s(time, by = s1_12, k = 5) + # 时间变化的傅立叶系数，以捕捉 # 可能变化的季节性。

当然，该模型仍然使用时间样条来模拟趋势，但单调性限制确保了预测结果至少在一定程度上是可控的。基于样本内拟合的效果，我们更倾向于选择第二个模型，并可以通过近似留一交叉验证（Approximate Leave-One-Out Cross-Validation）轻松地比较两个模型的性能。在这里，我们更倾向于具有更高期望对数预测密度（ELPD）的模型。

查看平滑局部效应图，了解某些傅立叶系数的估计值随时间的变化情况

绘制谐波回归的趋势和季节成分

绘制该模型的图需要更多的工作，因为我们需要使用该模型的目标预测来直观地显示两个主要效应（即非线性趋势和随时间变化的季节性）。

首先是趋势图，在这个模型中使用参数可以很容易地计算出模型的期望值（即结果标度上的预测值，但忽略了负二项采样分布的超分散参数）。

现在是时变季节模式，需要从总预测值中减去条件趋势预测值

season_preds <- with_season - agg_r_season

现在，我们将这些效果绘制在一起

这幅图清楚地显示了模型是如何估算出这两个组成部分的，以及季节模式在形状和幅度上随时间的变化情况

时变周期性

上例说明了如何利用傅立叶项来捕捉时变季节性。但如果周期性本身随时间变化呢？换句话说，如果一个序列开始时每 8 周出现一次有规律的周期性变化，但随后周期性窗口不断缩小或扩大，我们该如何建立模型呢？下面是一个非常简单的示例，说明我们如何扩展傅立叶策略，以对这种情况进行建模

首先，我定义了构建单调递增或递减系数的函数。这将有助于模拟时变周期性

  x <- sort(as.vector(scale(times)))
  exp(k * x) / (1 + exp(k * x))

现在进行模拟，包括构建两个不同周期的傅里叶级数（本例中为 “k = 12 “和 “k = 6″），每个周期有三个傅里叶项。然后，我模拟这些项的时变系数，迫使周期为 12 的项随时间下降为零，而周期为 6 的项随时间上升。这些系数加上模拟的傅立叶项本身，就可以模拟出一个周期性随时间从 ~12 缩小到 ~6 的时间序列。

绘制随时间变化的傅立叶系数图

现在计算线性预测因子并绘制曲线

最后，添加一些高斯随机漫步噪声，绘制模拟时间序列图

现在我们可以构建用于建模的数据；这需要我们仔细考虑问题，并将两个周期性的傅立叶级数作为预测因子纳入其中

现在我们可以拟合动态 GAM，它使用状态空间表示法来捕捉随时间变化的季节性、周期性和随机漫步自相关过程

和之前一样，绘制时变傅立叶系数图

draw(trend_mgcv_model)

和以前一样，我们的模型成功地捕捉到了这些项的时变系数。请注意，许多 12 期项（即 s1_12 和 c3_12）估计会随时间推移而下降，而许多 6 期项（即 s2_6 和 c3_6）估计会随时间推移而上升。令人鼓舞的是，该模型在处理随机漫步自相关过程和两个方差项（过程误差和观测误差）的同时，还能发现这些问题。

我们还可以使用后验期望绘制季节性时变图，就像上面第一个例子那样。但要注意的是，该图的预测不包括随机漫步部分。

现在，我们可以提取估算出的 RW 趋势，并减去季节性预测，这样就可以看到模型估算出的潜在非线性趋势了

treneds <- hinst(mod3, type = 'expected')$hiasts[[1]] - 
  sepreds

将二者绘制在一起。

参考文献

以下论文和资源提供了有关动态谐波回归和使用 GAM 进行时间序列建模/预测的有用资料

Autin, M. A., & Edwards, D. 河口水质数据的非参数谐波回归 Environmetrics 21.6 (2010)： 588-605.

Karunarathna, K. A. N. K., Wells, K., & Clark, N. J. (2024). 用分层动态广义加法模型模拟沙漠啮齿动物物种对环境变化的非线性响应. Ecological Modelling, 490, 110648.

Young, P. C., Pedregal, D. J., & Tych, W. (1999). Dynamic harmonic regression. Journal of Forecasting, 18, 369-394.