时间序列模型（一）：传统统计模型

下一篇：LSTM 深度解析 –>

本章要点

平稳性为什么是整个 ARIMA 家族的入场券，差分如何换来它。
像 Box-Jenkins 学派那样阅读 ACF / PACF：用 “截尾 vs 拖尾” 这条规则识别 $p$ 与 $q$。
ARIMA / SARIMA 的完整机器，以及季节性如何通过滞后 $s$ 算子被纳入模型。
VAR、GARCH、指数平滑、Prophet 与卡尔曼滤波如何被装进同一张地图：均值动态 vs. 方差动态 vs. 状态空间递推。
一条决策规则：什么时候传统模型就够了，什么时候必须升级到本系列后面的深度模型。

前置知识

基本的概率与统计（均值、方差、协方差、相关系数）。
熟悉 NumPy 和 pandas 的时间索引。
VAR / 卡尔曼小节会用到一点线性代数（矩阵乘法、特征值）。

1. 为什么传统模型仍然重要

在深度学习时代之前，时间序列工具箱已经相当完备。ARIMA 抓线性自相关，SARIMA 把日历效应补上，VAR 推广到多元，GARCH 描述方差动态，卡尔曼滤波则在状态空间框架下统一了上面所有人。它们共享三条深度模型并不免费提供的优点：

可解释性。 每个参数都有意义 – “昨天的水平以权重 $\phi_1$ 影响今天”，“两个月前的冲击以权重 $\theta_2$ 衰减”。
校准过的不确定性。 置信区间直接来自极大似然，而不是 Dropout 之类的近似手法。
样本效率。 几百个观测值就够，不需要 GPU。

如果你的序列短、平滑、或具备清晰的日历结构，传统模型经常打败 LSTM，而且更容易调试。请把本文当作 必须先打败的基线，再谈别的复杂方法。

2. 分解视角

任何建模之前，先建立一张图：大多数序列都可以写成

$$ y_t = T_t + S_t + R_t, $$

一条缓慢运动的趋势 $T_t$、一段周期为 $s$ 的 季节性 $S_t$，以及结构被剥离后看起来像噪声的残差 $R_t$。经典加法分解就是把它具体化：

图 1 – 经典加法分解。一旦趋势和季节性被减掉，残差就应当近似平稳白噪声 – 而这正是 ARMA 类模型可以良定义的区域。

整个 ARIMA 纲领可以用一句话概括：先变换数据让它看起来平稳，然后拟合一个带自相关误差的线性模型。

平稳性的形式定义

序列若 （弱）平稳，则其均值、方差与自协方差不依赖 $t$：

$$ \mathbb{E}[y_t] = \mu, \qquad \mathrm{Var}(y_t) = \sigma^2, \qquad \mathrm{Cov}(y_t, y_{t-k}) = \gamma_k. $$

绝大多数真实序列都不是平稳的 – 它们要么有趋势、要么漂移、要么方差随水平变大。两种标准疗法是：

差分：$\nabla y_t = y_t - y_{t-1}$ 去除线性趋势；$\nabla^2 y_t$ 去除二次趋势。
方差稳定化变换：$\log y_t$ 或 Box-Cox 变换可以驯服乘性增长。

ADF（增广 Dickey-Fuller）检验给出一个假设检验：$H_0$ 为 “存在单位根”，$p$ 值小就可以认为变换后的序列平稳。

3. AR、MA、ARMA – 三种记忆方式

ARIMA 是由两种基本组件搭起来的。它们看起来相像，但对 “序列记住了什么” 的看法完全不同。

图 2 – AR 记住过去的取值，MA 记住过去的冲击，ARMA 同时做这两件事。即便不算任何统计量，三者的样本路径在视觉上就有清晰差别。

自回归：AR($p$)

$$ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \varepsilon_t. $$

今天的取值是过去 $p$ 步取值的线性组合，再加一个白噪声冲击 $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$。“持久性"直接编码在系数 $\phi_k$ 里。$\phi_1$ 接近 1 的 AR(1) 会产生非常平滑、缓慢均值回复的轨迹。

移动平均：MA($q$)

$$ y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \cdots + \theta_q \varepsilon_{t-q}. $$

今天的取值是 最近 $q$ 个冲击 的线性组合。记忆窗口是有限的：超过 $q$ 步以前的事件对今天没有直接影响。

合体：ARMA($p$, $q$)

$$ \phi(B)\, y_t = \theta(B)\, \varepsilon_t, $$

其中 $B$ 是 滞后算子（$B y_t = y_{t-1}$），并且

$$ \phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p, \qquad \theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q. $$

ARMA 的好处是简洁：一组小的 $(p, q)$ 就能近似一个 MA($\infty$) 或 AR($\infty$) 的自相关结构。平稳性要求 $\phi(B) = 0$ 的所有根落在单位圆外；可逆性对 $\theta(B)$ 要求相同的事。

从 ARMA 到 ARIMA

如果原始序列不平稳，先取 $d$ 阶差分，再在差分序列上拟合 ARMA。完整记号是 ARIMA($p, d, q$)：

$$ \phi(B)\, (1-B)^d\, y_t = \theta(B)\, \varepsilon_t. $$

实践中 $d = 1$ 处理线性漂移，$d = 2$ 处理曲率，几乎不需要 $d > 2$。

4. ACF 与 PACF – 模型识别的显微镜

怎么挑 $p$ 和 $q$？两张诊断图基本就够了。

自相关函数 ACF：$\rho_k = \mathrm{Corr}(y_t, y_{t-k})$。
偏自相关 PACF：$y_t$ 与 $y_{t-k}$ 在 剔除了 $y_{t-1}, \ldots, y_{t-k+1}$ 的线性效应之后的相关。

Box-Jenkins 识别规则：

过程	ACF	PACF
AR($p$)	拖尾（光滑衰减）	在 $p$ 阶后截尾
MA($q$)	在 $q$ 阶后截尾	拖尾
ARMA($p$, $q$)	拖尾	拖尾

图 3 – AR(2) 的 PACF 在 2 阶后掉到 0（这种棒图是它的标志）。对偶地，MA(2) 的 ACF 在 2 阶后掉到 0。拿不准时优先选更简单的模型 – 第 5 节的诊断检验会抓住欠拟合。

如果两张图都没有干净的截尾，那基本就是 ARMA 过程；最干净的做法是对 $(p, q)$ 做网格搜索，按信息准则选最好的一对。

AIC 与 BIC

$$ \mathrm{AIC} = -2\,\ell(\hat{\theta}) + 2k, \qquad \mathrm{BIC} = -2\,\ell(\hat{\theta}) + k\log n, $$

其中 $\ell$ 是极大化的对数似然，$k$ 是自由参数个数，$n$ 是样本量。BIC 对复杂度惩罚更重，样本量小时优先用它。

5. Box-Jenkins 工作流

ARIMA 不是一次拟合就完事的，它是 Box 和 Jenkins 在 1970 年正式提出的 迭代循环。后来所有的统计预测工具（包括 auto.arima）都不过是在自动化同样这四个方框。

Box-Jenkins 方法论：识别、估计、诊断检验与预测，残差诊断不通过时回到识别。

图 4 – 真正的关键是那条反馈箭头。如果残差通不过 Ljung-Box 检验，或者残差 ACF 仍有结构，那就是模型设定错了 – 应回头改 $(p, d, q)$，而不是相信预测结果。

识别。 画图、做 ADF 看平稳性、研究 ACF/PACF、提出候选 $(p, d, q)$。
估计。 极大似然或条件最小二乘，主流库都帮你做好。
诊断检验。 残差是不是白噪声？
- Ljung-Box 统计量 $Q = n(n+2)\sum_{k=1}^h \hat{\rho}_k^2 / (n-k)$ 应不能拒绝 $H_0$。
- 残差 ACF 应落在 $\pm 1.96/\sqrt{n}$ 的带内。
- QQ 图与 Jarque-Bera 检验正态性（这关系到预测区间的可信度）。
预测。 残差看起来干净之后，再给出点预测和区间预测。

6. ARIMA 的代码实现

用 statsmodels 的最小完整流程：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.diagnostic import acorr_ljungbox

# 1. 模拟一个非平稳序列（带漂移的随机游走）
rng = np.random.default_rng(0)
n = 200
eps = rng.normal(0, 1.0, n)
y = np.zeros(n)
for t in range(1, n):
    y[t] = y[t - 1] + 0.05 + 0.6 * eps[t]

series = pd.Series(y, index=pd.date_range("2018-01-01", periods=n, freq="D"))

# 2. 识别：分别对原序列和一阶差分做 ADF
print("ADF on level :", adfuller(series)[1])                    # p 值大 -> 非平稳
print("ADF on diff  :", adfuller(series.diff().dropna())[1])    # p 值小 -> 平稳

# 3. 估计：用前 170 个点训练
train, test = series.iloc[:-30], series.iloc[-30:]
model = ARIMA(train, order=(2, 1, 1)).fit()

# 4. 诊断检验
lb = acorr_ljungbox(model.resid, lags=[10], return_df=True)
print(lb)  # p 值高 -> 残差是白噪声

# 5. 给出 95% 区间预测
fc = model.get_forecast(steps=30)
mean_fc = fc.predicted_mean
ci = fc.conf_int(alpha=0.05)

留出尾段做预测的效果是这样：

图 5 – 点预测会快速回到长期漂移线，区间则按 $\sigma\sqrt{h}$ 的速率随预测步长 $h$ 张开。这条逐渐张开的喇叭，是线性模型对未来诚实的认知边界。

7. 季节性：SARIMA

很多序列带着 ARIMA 自己处理不了的日历 – 12 月达峰的零售月度数据、有周周期的日度流量、有 24 小时周期的小时负荷。SARIMA($p, d, q$)($P, D, Q$)$_s$ 把周期为 $s$ 的季节滞后纳入：

$$ \Phi(B^s)\, \phi(B)\, (1-B)^d (1-B^s)^D y_t \;=\; \Theta(B^s)\, \theta(B)\, \varepsilon_t. $$

你做两种差分 – 常规差分 $(1-B)$ 去趋势，季节差分 $(1-B^s)$ 去周期 – 然后在常规滞后和 $s$ 的整数倍滞后处分别加 AR / MA 项。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
from statsmodels.tsa.statespace.sarimax import SARIMAX

model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),  # P, D, Q, s
    enforce_stationarity=False,
    enforce_invertibility=False,
).fit(disp=False)

forecast = model.get_forecast(steps=24)

图 6 – 季节周期在样本外被干净复现：季节差分算子 $(1-B^{12})$ 把年内模式从残差中剥掉，季节 AR/MA 项又把它在正确的相位上加回到预测里。

季节部分识别小抄：把 季节差分后 的序列在 $s, 2s, 3s, \ldots$ 处看 ACF/PACF，截尾 / 拖尾的判别规则与非季节情形一样。

8. 超出 ARIMA：家族其它成员

上面的思想沿着四个有用的方向被推广。它们并不是平行宇宙 – 每一个都只在 ARIMA 的某一维上做了专门化。

8.1 VAR – 多元动态

当你有多条相互影响的序列（GDP 与失业率、用电量与气温），就把标量 AR 升级成 向量自回归：

$$ \mathbf{y}_t = \mathbf{c} + A_1 \mathbf{y}_{t-1} + A_2 \mathbf{y}_{t-2} + \cdots + A_p \mathbf{y}_{t-p} + \boldsymbol{\varepsilon}_t. $$

矩阵 $A_k$ 的每一项都有意义：$(A_k)_{ij}$ 是变量 $j$ 在滞后 $k$ 处对今天变量 $i$ 的边际效应。这让 VAR 在宏观经济学里非常受欢迎，因为 Granger 因果（"$x$ 在 $y$ 自己的过去之外，是否还有助于预测 $y$？"）就是核心问题。

实务陷阱：$K$ 条序列、滞后 $p$ 时模型有 $K + pK^2$ 个自由参数。$K = 10, p = 4$ 已经是 410 个数，要从几百个观测里估出来。所以高维场景才有正则化的近亲 – 贝叶斯 VAR、因子模型。

8.2 GARCH – 方差动态

ARIMA 建模均值，GARCH 建模 条件方差。基础的 GARCH(1,1) 是：

$$ \sigma_t^2 = \omega + \alpha\, \varepsilon_{t-1}^2 + \beta\, \sigma_{t-1}^2, \qquad \varepsilon_t = \sigma_t\, z_t,\quad z_t \sim \mathcal{N}(0, 1). $$

$\alpha$ 项让昨天的大冲击能把今天的方差顶起来（“ARCH 效应”），$\beta$ 项让方差具有持续性。平稳性要求 $\alpha + \beta < 1$；金融数据里 $\alpha + \beta$ 通常落在 0.95-0.99，说明波动率是个慢变量。

GARCH 是风险管理（VaR、期权定价）的标配，并且天然可以与 ARIMA 均值模型搭配：先用 ARIMA 拟合收益率，再用 GARCH 拟合残差平方。

8.3 指数平滑与 Holt-Winters

如果说 ARIMA 是 “显式随机” 的视角，指数平滑就是 “加权平均” 的视角。简单指数平滑假设有一个水平 $\ell_t$ 随新信息漂移：

$$ \ell_t = \alpha\, y_t + (1-\alpha)\, \ell_{t-1}. $$

Holt 在此基础上加了趋势 $b_t$；Holt-Winters 又加了季节 $s_t$。加法形式为：

$$ \ell_t = \alpha (y_t - s_{t-s}) + (1-\alpha)(\ell_{t-1} + b_{t-1}),\\ b_t = \beta(\ell_t - \ell_{t-1}) + (1-\beta) b_{t-1},\\ s_t = \gamma(y_t - \ell_t) + (1-\gamma) s_{t-s}. $$

这是 R / statsmodels 中 ETS 家族背后的主力，赢得 M-竞赛的许多方法都是它的变种。

图 7 – 分解结构是 Holt-Winters 的迷人之处：每个分量本身就可读，平滑系数 ($\alpha, \beta, \gamma$) 告诉你每个分量适应得多快，递推形式让拟合开销只是 $O(n)$。

8.4 Prophet

Prophet（Facebook，2017）是一个面向业务分析师、刻意保持简单的 加法分解 模型：

$$ y(t) = g(t) + s(t) + h(t) + \varepsilon_t, $$

$g(t)$ 是带变点的分段线性或逻辑趋势，$s(t)$ 是用傅里叶展开表达多重季节性，$h(t)$ 编码用户提供的节假日。它通过 Stan 的 MAP 或全贝叶斯采样拟合，对外只暴露少量旋钮，对缺失值与离群点鲁棒。它在标准基准上并非最强，但 API 质量和节假日处理让它在产品分析里仍然是默认选择。

8.5 卡尔曼滤波与状态空间��角

上面所有模型，都是 线性高斯状态空间模型 的特例：

$$ \mathbf{x}_t = F_t \mathbf{x}_{t-1} + \mathbf{w}_t, \qquad \mathbf{w}_t \sim \mathcal{N}(0, Q_t),\\ \mathbf{y}_t = H_t \mathbf{x}_t + \mathbf{v}_t, \qquad \mathbf{v}_t \sim \mathcal{N}(0, R_t). $$

卡尔曼滤波 是精确的贝叶斯递推，在数据到达时维护 $p(\mathbf{x}_t \mid \mathbf{y}_{1:t}) = \mathcal{N}(\hat{\mathbf{x}}_t, P_t)$：

$$ \hat{\mathbf{x}}_{t|t-1} = F_t\, \hat{\mathbf{x}}_{t-1},\\ P_{t|t-1} = F_t\, P_{t-1}\, F_t^\top + Q_t,\\ K_t = P_{t|t-1} H_t^\top (H_t P_{t|t-1} H_t^\top + R_t)^{-1},\\ \hat{\mathbf{x}}_t = \hat{\mathbf{x}}_{t|t-1} + K_t (\mathbf{y}_t - H_t \hat{\mathbf{x}}_{t|t-1}),\\ P_t = (I - K_t H_t)\, P_{t|t-1}. $$

为什么这件事重要：ARIMA、指数平滑、动态线性回归、带季节哑元的结构模型 – 全部可以转写成上面的形式，免费继承卡尔曼递推。这正是 statsmodels 的 SARIMAX 在底层的实现路径，也是它能够无痛处理缺失观测和外生回归元的原因。

9. 怎么选模型

模型	适用	不适用
ARIMA	单序列、线性自相关、无日历	强季节性、多输入
SARIMA	单一明显季节周期、序列长度适中	多个重叠季节性
VAR	几条平稳序列存在反馈关系	高维（$K \gtrsim 10$）
GARCH	收益率 / 波动率数据，存在聚集	平滑、方差恒定的序列
Holt-Winters / ETS	鲁棒的季节基线、需要快速重训	高度非线性的状态切换
Prophet	带节假日、有缺失的业务序列	亚日级高频数据
卡尔曼 / 状态空间	在线更新、缺失数据、需要自定义结构	数据充足且可接受黑箱时

务实的处方：先用 ETS 或 SARIMA 当基线，关心方差时再加 GARCH；只有当残差呈现明显非线性结构、或你有大量并行序列需要共享信息时，再上深度模型（LSTM、TCN、Transformer、N-BEATS、Informer – 也就是本系列后面几篇）。

10. 局限与下一步

传统模型在以下情形会触顶：

动态本身就是 非线性 的（机制切换、阈值、硬饱和）。
你有 几百到几千条并行序列，应当共享参数。
相关上下文窗口非常长，ARMA 的根没法简洁地表示它。
信号 非高斯 到连 Box-Cox 都救不回来。

那正是后面七篇文章的起点。我们从 LSTM 这个非线性序列模型的代表开始，然后一路走到 GRU、注意力、Transformer、TCN、N-BEATS，最后到面向超长序列的 Informer。它们每一个都在推广上面某个线性模型；请把 ARIMA / SARIMA / Holt-Winters 当成 必须跨过去的基线，而不是可以跳过去的过去。

参考资料

Box, G., Jenkins, G., Reinsel, G., & Ljung, G. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press.
Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. https://otexts.com/fpp3/
Durbin, J., & Koopman, S. J. (2012). Time Series Analysis by State Space Methods (2nd ed.). Oxford University Press.
Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007.
Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3), 307-327.
Taylor, S. J., & Letham, B. (2018). Forecasting at scale. The American Statistician, 72(1), 37-45.

系列导航