ARIMAX 模型与单位根（DF/ADF）检验

多元时间序列分析：ARIMAX 模型与单位根检验

小象实战讲义 · 时间序列分析

在之前的学习中，我们专注于分析单个时间序列的内部动态结构，例如使用 AR、MA 或 ARMA 模型。然而，现实世界中的变量往往相互关联。例如，居民消费支出会受到收入水平的显著影响。为了更精准地预测和分析，我们需要将多个相关的时间序列变量同时纳入模型框架。本节将开启多元时间序列分析的大门，首先学习如何为平稳序列建立包含外生变量的动态回归模型（ARIMAX），并深入探讨一个至关重要的前置步骤——如何科学地检验序列的平稳性，以避免“虚假回归”的陷阱。

💡 核心导读

本节你将掌握：

1. ARIMAX 模型：理解如何将外生变量引入 ARIMA 框架，构建动态回归模型以分析多变量间的动态关系。
2. 虚假回归：认识对非平稳序列直接进行回归分析可能导致的严重统计谬误。
3. 单位根检验：学习使用 DF 检验这一统计方法，严谨地判断一个序列是否平稳，其核心是检验序列的特征根是否落在单位圆上。
4. DF 检验的三种类型：了解如何根据序列可能包含的确定性成分（如常数项、趋势项）选择合适的检验模型。
5. Python 实战：使用 statsmodels 库构建 ARIMAX 模型，并进行单位根检验，将理论应用于实际数据分析。

ARIMAX 模型：平稳序列的动态回归

在一元时间序列分析中，ARMA 模型用序列自身的过去值（AR 部分）和过去冲击的线性组合（MA 部分）来解释当前值。当我们希望引入其他时间序列变量（外生变量）来解释目标序列时，就需要扩展到多元框架。

对于一个平稳的响应序列 $y_t$ 和 $k$ 个平稳的输入序列 $x_{1t}, x_{2t}, \dots, x_{kt}$，我们可以建立 ARIMAX 模型（Autoregressive Integrated Moving Average with eXogenous inputs），也称为动态回归模型。其一般形式为：

$$y_t = \mu + \sum_{i=1}^{k} \frac{\Theta_i(B)}{\Phi_i(B)} B^{l_i} x_{it} + \eta_t$$

其中：

• $\mu$ 为常数项。
• $x_{it}$ 是第 $i$ 个外生输入序列。
• $B$ 是滞后算子（$B x_t = x_{t-1}$）。
• $l_i$ 是第 $i$ 个输入变量的延迟阶数。
• $\frac{\Theta_i(B)}{\Phi_i(B)}$ 是对第 $i$ 个输入变量的传递函数，用于刻画 $x_{it}$ 对 $y_t$ 影响的动态结构（可能包含自回归和移动平均成分）。
• $\eta_t$ 是回归模型的残差序列。由于 $y_t$ 和所有 $x_{it}$ 都是平稳的，$\eta_t$ 理论上也是一个平稳序列。

这个残差序列 $\eta_t$ 可能仍然包含自相关信息，因此可以进一步用一个 ARMA 模型来刻画：

$$\Phi(B) \eta_t = \Theta(B) a_t$$

其中 $a_t$ 是零均值的白噪声序列，即 $a_t \sim WN(0, \sigma_a^2)$。将两部分结合，就得到了完整的 ARIMAX 模型结构。

模型直观理解：ARIMAX 模型可以看作两部分之和。第一部分是外生变量 $x_t$ 及其滞后项对 $y_t$ 的线性影响（动态回归部分）；第二部分是对回归后残差序列 $\eta_t$ 中剩余自相关结构的 ARMA 建模。这比单独对 $y_t$ 拟合 ARMA 模型或简单地将 $y_t$ 对 $x_t$ 做回归都更有效，因为它同时捕捉了外部影响和内部动态。

虚假回归：对非平稳序列建模的陷阱

ARIMAX 模型要求所有序列（$y_t$ 和 $x_{it}$）都是平稳的。如果忽略这一前提，对非平稳序列直接应用最小二乘回归，可能导致 虚假回归。

考虑一个简单的线性回归模型： $$y_t = \beta_0 + \beta_1 x_t + u_t$$ 我们通常用 $t$ 检验来考察 $\beta_1$ 是否显著不为零（$H_0: \beta_1=0$）。

当 $x_t$ 和 $y_t$ 都是平稳序列且不相关时，$t$ 统计量服从 $t$ 分布，我们能够以预设的显著性水平（如 $\alpha=0.05$）控制犯第一类错误（拒真）的概率。

然而，如果 $x_t$ 和 $y_t$ 是非平稳序列（例如，都是随机游走过程），并且它们之间本来没有真正的相关关系，情况会截然不同。理论研究和模拟实验（如 Granger 和 Newbold 在 1974 年的著名研究）表明，在这种情况下：

• 回归系数 $\hat{\beta}_1$ 的标准误会被严重低估。
• 导致 $t$ 统计量的值异常地大。
• 使得我们拒绝 $\beta_1=0$ 这个真原假设的概率远高于设定的 $\alpha$（在模拟中可能高达 75%，而非 5%）。

这意味着，我们很容易错误地得出“$x_t$ 对 $y_t$ 有显著影响”的结论，而实际上这种关系是虚假的、由数据本身的非平稳性“伪造”出来的。因此，在建立任何时间序列回归模型之前，对序列进行平稳性检验是至关重要的第一步。

单位根检验：DF 检验

之前我们通过观察时序图来直观判断平稳性。现在引入严格的统计检验方法——单位根检验。其核心思想是：检验序列的特征根是否在单位圆上。若存在单位根，则序列非平稳。

最基础的检验是 DF 检验。它从一个一阶自回归模型 AR(1) 出发： $$x_t = \phi_1 x_{t-1} + \varepsilon_t, \quad \varepsilon_t \sim WN(0, \sigma^2)$$ 我们知道，该序列平稳的充要条件是 $|\phi_1| < 1$；非平稳（存在单位根）则对应 $\phi_1 = 1$。

因此，DF 检验的假设设置为：

• 原假设 $H_0: \phi_1 = 1$ （序列非平稳，存在单位根）
• 备择假设 $H_1: |\phi_1| < 1$ （序列平稳）

为了构造检验统计量，将模型改写为： $$\Delta x_t = \rho x_{t-1} + \varepsilon_t$$ 其中 $\Delta x_t = x_t - x_{t-1}$， $\rho = \phi_1 - 1$。此时的假设等价于：

• $H_0: \rho = 0$ （存在单位根）
• $H_1: \rho < 0$ （序列平稳）

检验统计量构造为： $$\tau = \frac{\hat{\rho}}{S(\hat{\rho})}$$ 其中 $\hat{\rho}$ 是 $\rho$ 的最小二乘估计，$S(\hat{\rho})$ 是其标准误。

在 $H_0$ 成立（序列有单位根）时，统计量 $\tau$ 的极限分布不是标准的 $t$ 分布或正态分布，而是一个由维纳过程函数构成的复杂分布。该分布无法求出解析的分位数，因此 Dickey 和 Fuller 通过蒙特卡洛模拟计算了其临界值表。检验时，将计算得到的 $\tau$ 值与 DF 临界值表比较，若 $\tau$ 值小于临界值（更负），则拒绝 $H_0$，认为序列平稳。

DF 检验的三种类型

实际序列可能包含确定性成分，因此 DF 检验扩展为三种类型，对应不同的数据生成过程：

1. 类型一（无常数项无趋势项）： $$\Delta x_t = \rho x_{t-1} + \varepsilon_t$$ 适用于不带漂移项的随机游走过程。序列均值非平稳，方差时变，但一阶差分后平稳。

2. 类型二（包含常数项）： $$\Delta x_t = \mu + \rho x_{t-1} + \varepsilon_t$$ 适用于带漂移项（常数项 $\mu$）的随机游走过程。序列既有趋势，波动性也随时间增强。

3. 类型三（包含常数项和趋势项）： $$\Delta x_t = \mu + \beta t + \rho x_{t-1} + \varepsilon_t$$ 适用于趋势平稳过程。序列的非平稳性主要来自确定性的时间趋势 $t$，剔除趋势后即为平稳序列。

选择建议：通常从包含趋势项和常数项的模型（类型三）开始检验。如果趋势项不显著，则简化为类型二；如果常数项也不显著，再简化为类型一。软件输出会直接给出 $p$ 值，若 $p$ 值大于显著性水平（如 0.05），则不能拒绝“存在单位根”的原假设，认为序列非平稳。

案例：中国农村居民收入与消费的 DF 检验

考虑 1978-2002 年中国农村居民家庭人均纯收入对数序列 $\ln(x_t)$ 和生活消费支出对数序列 $\ln(y_t)$。时序图显示两者均有明显的上升趋势，初步判断为非平稳序列。 分别对两个序列进行三种类型的 DF 检验，结果发现所有检验的 $p$ 值均远大于 0.05。因此，不能拒绝“序列存在单位根”的原假设，证实了 $\ln(x_t)$ 和 $\ln(y_t)$ 均为非平稳序列。这警示我们，若直接对这两个序列建立回归模型，很可能陷入虚假回归的陷阱。

Python 实战：构建 ARIMAX 模型与平稳性检验

下面我们通过一个模拟案例，演示如何用 Python 构建 ARIMAX 模型，并对序列进行单位根检验（ADF 检验，DF 检验的增强版，原理相通）。

"""
时间序列分析 6.1 课时配套测试代码
内容：多元时间序列与 ARIMAX 模型 (外生变量引入、动态回归与外推预测)
"""

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX

np.random.seed(42)
n = 150

# 1. 构造外生驱动变量 X_t (如市场推广支出) 与具有自相关残差的目标变量 Y_t (销售量)
# X_t \sim AR(1): X_t = 0.5 * X_{t-1} + e_{x,t}
x = np.zeros(n)
e_x = np.random.normal(0, 1.5, n)
for t in range(1, n):
    x[t] = 0.5 * x[t-1] + e_x[t]

# 真实外生效应: beta = 2.5
# 扰动项为一阶平稳 AR(1) 过程: \eta_t = 0.6 * \eta_{t-1} + e_{y,t}
eta = np.zeros(n)
e_y = np.random.normal(0, 1.0, n)
for t in range(1, n):
    eta[t] = 0.6 * eta[t-1] + e_y[t]

y = 2.5 * x + eta + 10.0  # 基础销售量 10.0

dates = pd.date_range(start='2022-01-01', periods=n, freq='D')
df = pd.DataFrame({'sales': y, 'ad_spend': x}, index=dates)

train_df = df.iloc[:-5]
test_df = df.iloc[-5:]

# 2. 拟合 ARIMAX 模型: 引入外生变量并联合估计 AR(1) 扰动
model = SARIMAX(train_df['sales'], exog=train_df['ad_spend'], order=(1, 0, 0))
res = model.fit(disp=False)

# 3. 提取参数
beta_est = res.params['ad_spend']
ar_est = res.params['ar.L1']
beta_se = res.bse['ad_spend']
beta_pval = res.pvalues['ad_spend']

# 4. 基于测试集的外生变量进行 5 步外推预测
pred_res = res.get_forecast(steps=5, exog=test_df[['ad_spend']])
preds = pred_res.predicted_mean
actuals = test_df['sales']
rmse = np.sqrt(np.mean((preds.values - actuals.values) ** 2))

print(f"训练样本数: {len(train_df)}, 预测步数: 5")
print(f"真实外生系数 beta = 2.5000, 真实 AR(1) 系数 phi = 0.6000")
print(f"ARIMAX 模型拟合参数:")
print(f"  外生变量系数 (ad_spend): {beta_est:.4f} (标准误: {beta_se:.4f}, p值: {beta_pval:.4e})")
print(f"  自回归系数 (ar.L1):      {ar_est:.4f} (p值: {res.pvalues['ar.L1']:.4e})")
print(f"  外生变量具有显著正向影响: {beta_pval < 0.05 and beta_est > 0}")
print(f"\n未来 5 步外推预测值 vs 真实值:")
for i in range(5):
    print(f"  Day {i+1}: 预测值 = {preds.iloc[i]:.4f}, 真实值 = {actuals.iloc[i]:.4f}, 绝对误差 = {abs(preds.iloc[i] - actuals.iloc[i]):.4f}")
print(f"5 步预测 RMSE: {rmse:.4f}")

运行结果输出：

训练样本数: 145, 预测步数: 5
真实外生系数 beta = 2.5000, 真实 AR(1) 系数 phi = 0.6000
ARIMAX 模型拟合参数:
  外生变量系数 (ad_spend): 2.5190 (标准误: 0.0631, p值: 0.0000e+00)
  自回归系数 (ar.L1):      0.9934 (p值: 0.0000e+00)
  外生变量具有显著正向影响: True

未来 5 步外推预测值 vs 真实值:
  Day 1: 预测值 = 13.5610, 真实值 = 12.7194, 绝对误差 = 0.8416
  Day 2: 预测值 = 7.2483, 真实值 = 7.5810, 绝对误差 = 0.3327
  Day 3: 预测值 = 3.7425, 真实值 = 4.2122, 绝对误差 = 0.4698
  Day 4: 预测值 = 8.9172, 真实值 = 9.9467, 绝对误差 = 1.0294
  Day 5: 预测值 = 10.6211, 真实值 = 11.8402, 绝对误差 = 1.2190
5 步预测 RMSE: 0.8468

📝 动手练一练

1. 概念辨析：请简述“虚假回归”产生的原因及其在实践中的危害。为什么对非平稳序列进行差分预处理是避免虚假回归的常见方法？
2. 检验选择：现有一个宏观经济时间序列（如 GDP），其时序图显示其随时间有大致线性的增长趋势，但围绕趋势线的波动幅度似乎相对稳定。在进行 DF 检验时，你认为应该优先选择三种类型中的哪一种？为什么？

参考答案：

1. 虚假回归产生的原因是，当非平稳序列（如随机游走）被用于经典回归模型时，参数估计量的标准误会严重低估，导致 $t$ 统计量膨胀，使得本不存在的相关性在统计上显得“显著”。其危害是导致研究者得出错误的变量关系结论，误导决策。差分可以消除序列中的单位根，使其变为平稳序列，从而满足经典回归模型的基本假设，避免虚假回归。
2. 应优先选择类型三（包含常数项和趋势项） 的 DF 检验。因为时序图显示了明确的线性趋势，类型三的模型设定 $\Delta x_t = \mu + \beta t + \rho x_{t-1} + \varepsilon_t$ 能够同时捕捉序列可能存在的趋势成分和单位根成分。如果趋势项 $\beta$ 不显著，可以再考虑简化的模型。

本章小结

本节我们迈出了多元时间序列分析的第一步。

要点回顾：

• ARIMAX 模型 是分析平稳多变量时间序列动态关系的强大工具，它结合了外生变量的回归影响和对残差序列的 ARMA 建模。
• 虚假回归 是对非平稳序列进行回归分析时潜伏的巨大风险，会导致统计推断完全失效。
• 单位根检验（DF检验） 是判断序列平稳性的统计基石。通过检验 AR 模型的特征根是否等于 1，来判断序列是否具有单位根（非平稳）。
• DF检验的三种类型 对应了数据可能存在的不同确定性成分（无常数/有常数/有趋势），正确选择类型是获得可靠检验结论的关键。

行动清单：

1. 检验先行：在尝试建立任何时间序列回归模型前，务必使用 ADF 检验等工具对所有涉及的序列进行平稳性检验。
2. 模型构建：对于平稳的多变量序列，尝试使用 statsmodels 的 SARIMAX 类构建 ARIMAX 模型，并解读外生变量的影响是否显著。
3. 规避陷阱：牢记虚假回归的教训，对于检验出的非平稳序列，考虑进行差分处理或学习下一节将介绍的“协整”分析方法，而非直接建模。

— 小象教研组