ARMA 模型:推导、识别与估计

AR、MA、ARMA 模型怎么从白噪声推出来。ACF/PACF 的识别规则。为什么 MLE 能用,似然曲面长什么样。

第 9 篇 · 共 12 篇约 12 分钟

ARMA 模型是平稳时间序列分析的主力。把它搞懂,应用计量里带时间戳的数据你就懂了一大半。数学上主要就是多项式代数,看一遍你就能在餐巾纸上重新推。难的是识别（怎么知道是 AR(2) 还是 MA(3)）和诊断检验（怎么知道你拟合的模型没问题）。

这篇文章从零推导 AR、MA、ARMA。参考 Hamilton (1994, 第 3 章) 和 Tsay (2010, 第 3 章)。想看更偏金融应用的可以翻 Brooks (2019, 第 5 章)。

起点:白噪声

我们从白噪声 $\varepsilon_t \sim \text{WN}(0, \sigma^2)$ 开始:零均值、常方差 $\sigma^2$ 、无自相关。我们要用它构造平稳过程。两种自然的变换给出两个基本积木。

自回归 (AR) 模型

AR( $p$ ) 过程说:今天的值是过去 $p$ 个值的线性组合,再加一个冲击。

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

更紧凑地写,

\Phi(L) y_t = c + \varepsilon_t

其中 $\Phi(L) = 1 - \phi_1 L - \phi_2 L^2 - \cdots - \phi_p L^p$ 是滞后算子 $L$ 下的自回归多项式, $L y_t = y_{t-1}$ 。

平稳性条件: $\Phi(z) = 0$ 的所有根必须在复平面的单位圆外。等价地说,多项式 $1 - \phi_1 z - \cdots - \phi_p z^p \neq 0$ 对所有 $|z| \leq 1$ 成立。

如果 $\Phi(L)$ 可逆(根在单位圆外时就可逆),我们可以写:

y_t = \mu + \Psi(L) \varepsilon_t = \mu + \sum_{j=0}^\infty \psi_j \varepsilon_{t-j}

其中 $\mu = c / \Phi(1) = c / (1 - \phi_1 - \cdots - \phi_p)$ 是无条件均值,MA 系数 $\psi_j$ 几何衰减。

Yule-Walker 方程给出自协方差的闭式解。把 AR 方程两边乘 $y_{t-h}$ 再取期望:

\gamma_h = \phi_1 \gamma_{h-1} + \phi_2 \gamma_{h-2} + \cdots + \phi_p \gamma_{h-p} \quad \text{for } h \geq 1

除以 $\gamma_0$ 得到自相关递推:

\rho_h = \phi_1 \rho_{h-1} + \phi_2 \rho_{h-2} + \cdots + \phi_p \rho_{h-p} \quad \text{for } h \geq 1

其中 $\rho_0 = 1$ 。

特例. AR(1): $y_t = c + \phi_1 y_{t-1} + \varepsilon_t$ , $|\phi_1| < 1$ 。

均值: $\mu = c / (1 - \phi_1)$ 。方差: $\gamma_0 = \sigma^2 / (1 - \phi_1^2)$ 。自相关: $\rho_h = \phi_1^h$ 。几何衰减。

如果 $\phi_1 > 0$ ,自相关为正、平滑衰减。如果 $\phi_1 < 0$ ,正负交替。无论哪种情况,ACF 都是拖尾。

特例. AR(2): $y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \varepsilon_t$ 。

特征方程是 $1 - \phi_1 z - \phi_2 z^2 = 0$ 。平稳性要求两个根都在单位圆外。三种状态:

两个实根(都正、都负、或异号)→ ACF 按两个指数混合衰减。
共轭复根 → ACF 呈现阻尼正弦振荡。

这就是宏观数据里”伪周期”的来源:一个带复根的平稳 AR(2) 看起来就是周期性的,虽然它根本没有商业周期机制。

移动平均 (MA) 模型

MA( $q$ ) 过程说:今天的值是当前和过去 $q$ 个冲击的线性组合。

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

即,

y_t = \mu + \Theta(L) \varepsilon_t

其中 $\Theta(L) = 1 + \theta_1 L + \cdots + \theta_q L^q$ 。

MA 过程总是平稳的,对 $\theta_1, \ldots, \theta_q$ 任何取值都成立(因为 MA 表示的 $|\psi_j|$ 是有限和)。唯一约束是可逆性:要从观测到的 $y_t$ 还原冲击, $\Theta(z) = 0$ 的所有根必须在单位圆外。

对 $|\theta_1| < 1$ 的 MA(1):

y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1}

ACF 在滞后 1 之后截尾: $\rho_0 = 1$ , $\rho_1 = \theta_1 / (1 + \theta_1^2)$ , $\rho_h = 0$ 对 $h \geq 2$ 。

这个截尾是 MA 过程的诊断指纹。跟 AR 正好相反,AR 的 ACF 是拖尾。

ARMA( $p, q$ ) 模型

把两者组合起来:

y_t = c + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}

即,

\Phi(L) y_t = c + \Theta(L) \varepsilon_t

平稳性要求 $\Phi(z)$ 所有根在单位圆外。可逆性要求 $\Theta(z)$ 所有根在单位圆外。

ARMA( $p, q$ ) 的 ACF 在滞后 $q$ 之后拖尾(像 MA),PACF 在滞后 $p$ 之后拖尾(像 AR)。这种双重行为是 ARMA 灵活的原因。

识别: $(p, q)$ 阶数必须从数据识别出来。同一个过程没有唯一的 ARMA 表示。给 $\Phi(L)$ 加个因子、再从 $\Theta(L)$ 里去掉同一个因子,得到的过程还是一样的。约定是要求两个多项式的根都在单位圆外,这样因果性（平稳性 + 可逆性）成立。在这个约束下,表示是唯一的。

ARMA 的 Yule-Walker 方程

对 ARMA( $p, q$ ),自协方差对 $h > q$ 满足扩展 Yule-Walker 递推:

\gamma_h = \phi_1 \gamma_{h-1} + \phi_2 \gamma_{h-2} + \cdots + \phi_p \gamma_{h-p}

对 $1 \leq h \leq q$ ,关系里同时有 $\gamma$ 和 $y$ 与 $\varepsilon$ 的互协方差。闭式解需要求块 Toeplitz 矩阵的逆。实际上你不会手算这个,让软件去做。

识别:ACF 和 PACF

Box-Jenkins (1976) 识别规则:

模型	ACF	PACF
AR( $p$ )	拖尾	滞后 $p$ 之后截尾
MA( $q$ )	滞后 $q$ 之后截尾	拖尾
ARMA( $p, q$ )	滞后 $q$ 之后拖尾	滞后 $p$ 之后拖尾
白噪声	无显著自相关	无显著偏自相关

PACF(偏自相关函数)在滞后 $h$ 处,是去掉中间滞后项的线性影响之后, $y_t$ 和 $y_{t-h}$ 的相关。计算方法是逐个拟合 AR( $k$ ) 模型( $k = 1, 2, \ldots$ ),读 $y_{t-k}$ 的系数。等价地,解 Yule-Walker 方程组:

\begin{pmatrix} \rho_1 \\ \rho_2 \\ \vdots \\ \rho_h \end{pmatrix} = \begin{pmatrix} 1 & \rho_1 & \cdots & \rho_{h-1} \\ \rho_1 & 1 & \cdots & \rho_{h-2} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{h-1} & \rho_{h-2} & \cdots & 1 \end{pmatrix} \begin{pmatrix} \phi_{h,1} \\ \phi_{h,2} \\ \vdots \\ \phi_{h,h} \end{pmatrix}

最后一个系数 $\phi_{h,h}$ 就是滞后 $h$ 的 PACF。

直观理解:PACF 告诉你 $y_{t-h}$ 在 $y_{t-1}, \ldots, y_{t-h+1}$ 之外额外的预测能力。如果真实模型是 AR( $p$ ),超过 $p$ 的滞后就什么也不加了。

估计:OLS、Yule-Walker 和 MLE

三种主要估计策略:

1. OLS(最小二乘)。对 AR( $p$ ),把 $y_t$ 对 $y_{t-1}, \ldots, y_{t-p}$ 加常数做回归。标准 OLS 理论适用,因为回归子是预定的(滞后因变量),误差项 $\varepsilon_t$ 条件均值为零。估计量一致且渐近正态,但小样本偏差众所周知、可能很严重(Marriott & Pope, 1954; Shaman & Stine, 1988)。对 MA( $q$ ),OLS 不行,因为不可观测的冲击出现在回归子中。

2. Yule-Walker(矩估计)。用样本自协方差替换 Yule-Walker 方程里的理论自协方差然后求解。计算便宜但效率低,只用二阶矩信息,没用上整个分布。

3. 极大似然。对高斯 ARMA,用创新递推写出精确似然:

\hat{\varepsilon}_t = y_t - \hat{c} - \hat{\phi}_1 y_{t-1} - \cdots - \hat{\phi}_p y_{t-p} - \hat{\theta}_1 \hat{\varepsilon}_{t-1} - \cdots - \hat{\theta}_q \hat{\varepsilon}_{t-q}

把回溯值 $\hat{\varepsilon}_0, \hat{\varepsilon}_{-1}, \ldots$ 初始化为零,或者用无条件矩初始化。对数似然:

\ell(\theta) = -\frac{T}{2} \log(2\pi) - \frac{T}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \sum_{t=1}^T \hat{\varepsilon}_t^2

对 $\theta = (c, \phi_1, \ldots, \phi_p, \theta_1, \ldots, \theta_q, \sigma^2)$ 求极大。把 $\sigma^2$ 集中出去得到:

\hat{\sigma}^2 = \frac{1}{T} \sum_{t=1}^T \hat{\varepsilon}_t^2

把优化降到 $(p + q + 1)$ 维参数。条件 MLE 让 $\hat{\varepsilon}_t = 0$ 对 $t \leq 0$ 。精确 MLE 用卡尔曼滤波或 diffuse 初始化算精确似然,对近单位根或短序列效果更好。

MA 和 ARMA 的似然曲面出了名的粗糙（多个局部极大、山脊、单位圆附近的近平坦区域）。Box-Jenkins 建议先用 Yule-Walker 给初值,再用 MLE 精修。现代软件大多用 CSS(条件平方和)起手再接 MLE。

诊断检验

拟合完之后,检查残差是不是看起来像白噪声。标准套件:

1. 残差 ACF / PACF。算 $\hat{\varepsilon}_t$ 的样本自相关。在模型正确的零假设下, $\hat{\rho}_h \sim N(0, 1/T)$ 对 $h \geq 1$ 。画出 $\pm 1.96/\sqrt{T}$ 的置信带。哪个滞后的自相关穿出去,就是模型误设的信号。

2. Ljung-Box 检验。Q 统计量:

Q(m) = T(T+2) \sum_{h=1}^m \frac{\hat{\rho}_h^2}{T - h} \xrightarrow{d} \chi^2_{(m - p - q)}

检验前 $m$ 个残差自相关为零的联合零假设。 $Q$ 超过临界值就拒绝。用 $m \approx \log(T)$ 、 $m = 10$ 或 $m = 20$ 。

3. 正态性。对残差跑 Jarque-Bera、Shapiro-Wilk。ARMA 的推断假设高斯性;如果残差是肥尾的,你需要肥尾创新分布或者换模型类。

4. 异方差。对残差跑 ARCH-LM 检验(Engle, 1982)。如果被拒绝,说明你的 ARMA 没吸收完条件异方差,要在它上面再叠一层 GARCH。

5. 参数稳定性。递推最小二乘、CUSUM 检验(Brown, Durbin, Evans, 1975)。参数在样本里漂移的话,说明有结构突变,该考虑状态切换模型。

代码:Python 拟合 ARMA

import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA


rng = np.random.default_rng(0)
n = 1000
phi = [0.5, -0.3]
theta = [0.4]
eps = rng.standard_normal(n)
y = np.zeros(n)
for t in range(2, n):
    y[t] = 0.5*y[t-1] - 0.3*y[t-2] + eps[t] + 0.4*eps[t-1]

# 用 MLE 拟合 ARMA(2,1)
model = ARIMA(y, order=(2, 0, 1), trend='c').fit()
print(model.summary())

# 诊断:残差 ACF 和 Ljung-Box
from statsmodels.stats.diagnostic import acorr_ljungbox
lb = acorr_ljungbox(model.resid, lags=[10, 20], return_df=True)
print(lb)

拟合的 $\hat{\phi}_1, \hat{\phi}_2, \hat{\theta}_1$ 应该接近真实的 $(0.5, -0.3, 0.4)$ 。Ljung-Box 在常规水平下不应该拒绝。

ARMA 做不到的事

ARMA 用于平稳、线性、短记忆过程。它处理不了:

单位根: $y_t = y_{t-1} + \varepsilon_t$ 有 $\phi_1 = 1$ ,非平稳。对一阶差分跑 ARMA 可以(ARIMA,下一篇文章)。
长记忆:自相关双曲线衰减而不是几何衰减(ARFIMA,分数阶积分)。
非线性:用 GARCH 处理方差,马尔可夫切换处理状态,门槛模型、神经网络也行。
结构突变:参数在已知或未知日期发生变化。
多元动态:序列之间的互相关(VAR、VECM)。

每种情况都有专门的模型。接下来 11 篇大多跟其中之一相关。

小样本偏差问题

AR 系数的 OLS 估计在有限样本下有偏,正 AR 系数是向下偏。偏差对 AR(1) 系数 $\phi$ 大约是 $-2\phi/T$ (Kendall, 1954; Marriott & Pope, 1954)。 $T = 100$ 、 $\phi = 0.9$ 时偏差是 $-0.018$ ,很小。 $T = 30$ 、 $\phi = 0.9$ 时偏差是 $-0.06$ ,不能忽视。

对近单位根过程( $\phi \to 1$ ),偏差很严重。Shaman & Stine (1988) 显示,高持续性序列在小样本里偏差可能超过真实系数的一半。意思是:用 30 个观测、真实 $\phi = 0.95$ 的序列拟合 AR(1),OLS 估计可能给出 0.4 或 0.5,离真相远得很。

偏差修正是可能的(参见 Marriott & Pope, 1954,或 Shaman & Stine, 1988 的刀切法偏差修正),但标准做法就是承认偏差、尽量用长样本。

ARMA 的 MLE 渐近也有同样的偏差。对 MA 和 ARMA,偏差可能比 AR 还糟,因为 MA 多项式可以接近单位圆,而且似然曲面在边界附近有个尖峰。

用 ARMA 做预测

样本末端的 $h$ 步预测是:

\hat{y}_{T+h} = E[y_{T+h} \mid y_T, y_{T-1}, \ldots]

递归地算。对 AR(2):

\hat{y}_{T+1} = \hat{c} + \hat{\phi}_1 y_T + \hat{\phi}_2 y_{T-1}

\hat{y}_{T+2} = \hat{c} + \hat{\phi}_1 \hat{y}_{T+1} + \hat{\phi}_2 y_T

以此类推。 $h > p$ 时,所有滞后值本身就是预测,所以递推里用的是 $\hat{y}_{T+h-j}$ 对 $j > p$ 。

预测误差方差:

\text{Var}(y_{T+h} - \hat{y}_{T+h}) = \sigma^2 \sum_{j=0}^{h-1} \hat{\psi}_j^2

其中 $\hat{\psi}_j$ 是无限 MA 表示里的 MA 系数。对平稳 ARMA, $\sum \hat{\psi}_j^2$ 收敛到 $\sigma^{-2} \gamma_0$ (无条件方差),所以预测方差有渐近上限。

对近单位根 AR,预测方差按 $h \sigma^2$ 增长,比平稳情况快得多。

95% 预测区间:

\hat{y}_{T+h} \pm 1.96 \sqrt{\text{Var}(y_{T+h} - \hat{y}_{T+h})}

这假设创新项高斯。创新项肥尾的话,区间该更宽。Bootstrap 预测区间的处理方法:

拟合模型,拿残差 $\hat{\varepsilon}_t$ 。
有放回重抽样残差得 $\varepsilon_t^*$ 。
模拟 $\hat{y}_t^* = \hat{c} + \hat{\phi}_1 \hat{y}_{t-1}^* + \cdots + \hat{\varepsilon}_t^*$ 。
重复 $B = 1000$ 次。
$\hat{y}_{T+h}^*$ 的 2.5% 和 97.5% 分位数就是 bootstrap 区间。

ARMA 与有效市场假说

一个微妙但重要的点:ARMA 模型对资产收益是否可预测这件事没有发言权。如果收益是 iid(白噪声),最佳 ARMA 拟合是 ARMA(0, 0),也就是一个常数均值加白噪声残差。系数估计为零,残差方差等于样本方差,预测是平的。

如果 ARMA 拟合给出显著的 AR 或 MA 系数,说明存在某种可预测性:过去收益有助于预测未来收益。这种可预测性能不能用来赚钱,取决于交易成本、可预测性的幅度、以及它是否在样本外持续。

有效市场假说说收益不能用公开信息预测。ARMA 拟合是检验方法之一。但 EMH 是关于收益与条件信息联合分布的论断,不只是关于线性可预测性。市场可以是有效的(没有可利用的模式),即便 ARMA 系数非零,只要可预测性小到跨不过交易成本。

季节 ARMA

对带强季节性的数据(月度、季度),ARMA 模型自然延伸。最简单的形式是乘积季节 ARMA:

\Phi(L) \tilde{\Phi}(L^s) y_t = c + \Theta(L) \tilde{\Theta}(L^s) \varepsilon_t

其中 $s$ 是季节周期(月度数据为 12、季度数据为 4), $\tilde{\Phi}, \tilde{\Theta}$ 是季节 AR 和 MA 多项式。

滞后 $s$ 的季节 AR(1) 长这样:

y_t = c + \tilde{\phi}_s y_{t-s} + \varepsilon_t

意思是:今天的值依赖于 $s$ 期前的值。对月度零售数据(年度周期),这把今年 1 月和去年 1 月联系起来。

滞后 $s$ 的季节 MA(1):

y_t = c + \varepsilon_t + \tilde{\theta}_s \varepsilon_{t-s}

意思是:今天的值依赖于 $s$ 期前发生的冲击,且持续 $s$ 期。

识别:季节 ARMA 的 ACF 在季节滞后(月度年度周期是 12、24、36)上有尖峰。如果只有滞后 12 显著,就是季节 MA(1)。如果滞后 12 和 24 都显著,可能就是季节 AR(2)。

带差分的季节 ARIMA,在模型里加上季节差分 $(1 - L^s)$ :

\Phi(L) \tilde{\Phi}(L^s) (1 - L)^d (1 - L^s)^D y_t = c + \Theta(L) \tilde{\Theta}(L^s) \varepsilon_t

这就是 SARIMA $(p, d, q) \times (P, D, Q)_s$ 模型。它处理趋势(通过 $d$ )、季节趋势(通过 $D$ )、常规滞后和季节滞后上的 ARMA 动态。

对大多数月度经济数据,SARIMA $(0, 1, 1)(0, 1, 1)_{12}$ 是一个很强的默认起点。对强季节性的季度数据,SARIMA $(0, 1, 1)(0, 1, 1)_4$ 通常拟合得不错。

分数阶 ARIMA,简单提一下

ARIMA 类假设 $d$ 是整数。ARFIMA(Autoregressive Fractionally Integrated Moving Average,分数阶整合移动平均)允许 $d$ 不是整数:

\Phi(L) (1 - L)^d y_t = c + \Theta(L) \varepsilon_t

其中 $(1 - L)^d$ 是分数阶差分算子,用二项式级数定义:

(1 - L)^d = \sum_{j=0}^\infty \binom{d}{j} (-1)^j L^j = 1 - dL + \frac{d(d-1)}{2!} L^2 - \cdots

对 $0 < d < 0.5$ ,过程平稳但是长记忆:自相关按双曲线衰减( $\rho_h \sim h^{2d-1}$ )而不是几何衰减。这跟波动率、河流流量、互联网流量的经验证据相符。

代价是:ARFIMA 比 ARIMA 难估计, $d$ 估计量的渐近分布是非标准的(向正态收敛得很慢)。多数应用里,整数 $d$ 的 ARIMA 就够了。

参考文献

Box, G. E. P., & Jenkins, G. M. (1976). Time Series Analysis: Forecasting and Control. Holden-Day.
Brooks, C. (2019). Introductory Econometrics for Finance (4th ed.). Cambridge University Press. Ch. 5-6.
Brown, R. L., Durbin, J., & Evans, J. M. (1975). Techniques for testing the constancy of regression relationships over time. Journal of the Royal Statistical Society B, 37(2), 149-192.
Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007.
Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Ch. 3.
Ljung, G. M., & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297-303.
Marriott, F. H. C., & Pope, J. A. (1954). Bias in the estimation of autocorrelations. Biometrika, 41(3/4), 390-402.
Shaman, P., & Stine, R. A. (1988). The bias of autoregressive coefficient estimators. Journal of the American Statistical Association, 83(403), 842-848.
Tsay, R. S. (2010). Analysis of Financial Time Series (3rd ed.). Wiley. Ch. 3.
Kendall, M. G. (1954). Note on bias in the estimation of autocorrelation. Biometrika, 41(3/4), 403-404.
Granger, C. W. J., & Joyeux, R. (1980). An introduction to long-memory time series models and fractional differencing. Journal of Time Series Analysis, 1(1), 15-29.