残差自回归模型

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

在上一节中，我们系统学习了 ARIMA 模型及其在非平稳序列建模中的应用。然而，当我们使用确定性因素分解方法（如提取趋势项和季节项）后，得到的残差序列可能依然存在自相关性，这意味着模型未能充分提取序列中的信息。本节将聚焦于 残差自回归模型 (Auto-Regressive Model)，其核心思想是：先通过确定性模型提取主要趋势和季节效应，再对残差序列拟合自回归模型，以充分挖掘其内部的自相关结构。掌握此模型，你将能更精细地处理那些具有复杂确定性成分和随机波动的时间序列数据。

💡 核心导读

• 模型动机：理解为何在提取确定性因素后，残差序列的自相关性仍需建模。
• 模型结构：掌握残差自回归模型的数学表达式及其构成部分（确定性部分 + 残差自回归部分）。
• 检验方法：学习使用 DW 检验 和 Durbin h 检验 来判断残差序列是否存在一阶自相关。
• 建模流程：熟悉从趋势拟合、残差检验到自回归模型拟合的完整步骤。
• 模型比较：了解如何通过 AIC 和 BIC 准则在不同结构的模型（如 ARIMA 与残差自回归模型）间进行选择。

模型构造思想与结构

为什么需要残差自回归模型？

在经典的时间序列分解中，我们常将序列 $x_t$ 分解为确定性部分（如趋势 $T_t$、季节 $S_t$）和随机部分（$\varepsilon_t$）： $$x_t = T_t + S_t + \varepsilon_t$$ 理想情况下，残差 $\varepsilon_t$ 应是一个纯随机序列（白噪声），即 $\varepsilon_t \sim WN(0, \sigma^2)$，且不同期之间不相关：$Cov(\varepsilon_t, \varepsilon_{t-j}) = 0, \forall j \ge 1$。

然而，在实际数据分析中，尤其是经济、金融等领域，确定性因素的提取往往不够充分或精确。这导致残差序列 $\varepsilon_t$ 中可能仍包含未被提取的短期自相关信息。如果忽略这种自相关性，直接使用普通最小二乘法 (OLS) 进行趋势拟合，会导致参数估计的标准误有偏，进而影响假设检验（如 t 检验、F 检验）的可靠性，最终降低模型的预测精度。

残差自回归模型 正是为了解决这一问题而提出的。其基本思路是：对确定性模型拟合后的残差序列 $\varepsilon_t$，再建立一个自回归模型 (AR(p))，以提取其内部的自相关结构。最终，使得新的扰动项 $a_t$ 满足白噪声假设。

模型数学结构

残差自回归模型的一般形式如下：

$$\begin{cases} x_t = T_t + S_t + \varepsilon_t, & \text{(确定性部分)} \\ \varepsilon_t = \phi_1 \varepsilon_{t-1} + \phi_2 \varepsilon_{t-2} + \cdots + \phi_p \varepsilon_{t-p} + a_t, & \text{(残差自回归部分)} \\ a_t \sim WN(0, \sigma_a^2), \quad Cov(a_t, a_{t-j}) = 0, \forall j \ge 1. & \text{(白噪声假定)} \end{cases}$$

其中：

• $T_t$ 和 $S_t$ 分别代表趋势效应和季节效应，可以通过多种方法拟合（下文详述）。
• $\varepsilon_t$ 是提取确定性因素后的残差序列。
• $\phi_1, \phi_2, \dots, \phi_p$ 是自回归系数，$p$ 为自回归阶数。
• $a_t$ 是经过自回归模型提取后的最终残差，要求是零均值、同方差且无自相关的白噪声序列。

该模型表明，序列 $x_t$ 不仅受确定性因素影响，其随机波动部分 ($\varepsilon_t$) 自身也具有记忆性，可以通过其历史值进行预测。

确定性趋势与季节效应的拟合方法

在构建残差自回归模型时，首先需要选择合适的方法来拟合确定性部分 $T_t$ 和 $S_t$。

趋势效应 $T_t$ 的拟合

1. 时间 $t$ 的幂函数模型： 将趋势视为时间的函数，通常采用多项式形式。 $$T_t = \beta_0 + \beta_1 t + \beta_2 t^2 + \cdots + \beta_k t^k + \eta_t$$ 其中 $\eta_t$ 为随机误差。当序列呈现线性趋势时，取 $k=1$；呈现曲线趋势时，可取 $k=2$ 或更高。

2. 延迟因变量回归模型： 将当期趋势与历史观测值建立回归关系。这是一种动态的拟合方式。 $$T_t = \beta_0 + \beta_1 x_{t-1} + \beta_2 x_{t-2} + \cdots + \beta_k x_{t-k} + \eta_t$$ 例如，$T_t = \beta_0 + \beta_1 x_{t-1}$ 表示趋势由上一期值线性决定。

季节效应 $S_t$ 的拟合

1. 季节指数法： 这是确定性季节分析中的经典方法。为每个季节（如月度、季度）计算一个固定的季节指数 $S_i'$，然后赋值给对应的时间点。 $$S_t = S_i', \quad \text{当 } t \text{ 属于第 } i \text{ 个季节时}.$$

2. 季节自回归模型： 适用于季节效应与历史同期值相关的情况。假设周期为 $s$（如月度数据 $s=12$）。 $$S_t = \alpha_0 + \alpha_1 x_{t-s} + \alpha_2 x_{t-2s} + \cdots + \alpha_l x_{t-l \cdot s} + \eta_t$$ 该模型捕捉的是以周期 $s$ 为步长的长期季节相关性。

选择建议：对于趋势明显且简单的序列，方法1（时间幂函数）更直观易解释。对于趋势复杂或与自身历史值紧密相关的序列，方法2（延迟因变量）可能拟合效果更好。季节部分同理，需根据序列的周期性特征选择。

残差自相关性的检验：DW 检验与 Durbin h 检验

在决定是否需要对残差序列 $\varepsilon_t$ 进行自回归建模前，必须先检验其是否存在自相关性。这里我们重点介绍两种针对一阶自相关的检验方法。

Durbin-Watson (DW) 检验

原假设 $H_0$：残差序列不存在一阶自相关性，即 $\rho_1 = Corr(\varepsilon_t, \varepsilon_{t-1}) = 0$。 备择假设 $H_1$：残差序列存在一阶自相关性，即 $\rho_1 \ne 0$。

DW 统计量 构造如下： $$DW = \frac{\sum_{t=2}^n (\varepsilon_t - \varepsilon_{t-1})^2}{\sum_{t=1}^n \varepsilon_t^2}$$ 可以证明，在大样本下，$DW \approx 2(1 - \hat{\rho}_1)$，其中 $\hat{\rho}_1$ 是残差一阶自相关系数的估计值。

判定准则：

• $DW \approx 2$：无法拒绝 $H_0$，认为无自相关。
• $DW$ 显著小于 2（接近 0）：存在正自相关。
• $DW$ 显著大于 2（接近 4）：存在负自相关。 具体的临界值需要查 Durbin-Watson 检验表，通过比较 $DW$ 值与 $d_L$, $d_U$ 等临界值来判断。

DW 检验的局限：当回归模型的自变量中包含延迟因变量（如 $x_{t-1}$）时，DW 统计量是有偏的，倾向于低估自相关的程度，可能导致犯第二类错误（即实际存在自相关却未能检出）。此时应使用 Durbin h 检验。

Durbin h 检验

为了解决 DW 检验在模型包含滞后因变量时的偏误问题，Durbin 提出了 h 检验。

Dh 统计量： $$Dh = \hat{\rho}_1 \sqrt{\frac{n}{1 - n \cdot \widehat{Var}(\hat{\beta}_1)}}$$ 其中，$\hat{\rho}_1$ 是残差一阶自相关系数，$\hat{\beta}_1$ 是模型中滞后因变量（如 $x_{t-1}$）的系数估计值，$\widehat{Var}(\hat{\beta}_1)$ 是其方差的估计。在大样本下，$Dh$ 统计量近似服从标准正态分布 $N(0,1)$。

检验步骤：计算 $Dh$ 值后，与标准正态分布的临界值（如 $\pm 1.96$，对应 $\alpha=0.05$）比较，或直接计算 p 值进行判断。

注意：DW 和 Dh 检验主要针对一阶自相关。若检验不显著（接受 $H_0$），仅能说明不存在一阶自相关，但仍可能存在高阶自相关，此时需要借助残差序列的 ACF/PACF 图 或 Ljung-Box 检验 进行更全面的诊断。

实战演练：中国农业实际国民收入指数分析

我们以 1952-1988 年中国农业实际国民收入指数序列为例，演示残差自回归模型的完整建模流程，并与 ARIMA 模型进行比较。该序列具有明显的线性递增趋势，无显著季节效应。

步骤 1：趋势拟合与残差计算

首先，我们尝试两种趋势拟合方法。

方法一：时间 $t$ 的线性函数 拟合模型：$x_t = \beta_0 + \beta_1 t + \varepsilon_t$

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.stats.stattools import durbin_watson
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 假设已有数据 `df`，包含 'year' 和 'income_index' 列
# 为演示，我们生成一个具有类似特征的模拟序列
np.random.seed(42)
n = 37  # 1952-1988共37年
t = np.arange(1, n+1)
# 生成具有线性趋势和自相关残差的序列
true_trend = 66.1491 + 4.5158 * t
# 生成具有 AR(2) 结构的残差
resid_ar = np.zeros(n)
phi1, phi2 = 1.4859, -0.5848
for i in range(2, n):
    resid_ar[i] = phi1 * resid_ar[i-1] + phi2 * resid_ar[i-2] + np.random.normal(0, 3)
income_index = true_trend + resid_ar
df = pd.DataFrame({'year': np.arange(1952, 1989), 't': t, 'income_index': income_index})

# 方法一：线性趋势拟合
X1 = sm.add_constant(df['t'])  # 添加常数项
model1 = sm.OLS(df['income_index'], X1).fit()
df['epsilon1'] = model1.resid  # 残差序列 ε_t
print("=== 方法一：线性趋势模型 ===")
print(model1.summary())
print(f"DW 统计量: {durbin_watson(model1.resid):.4f}")

运行结果摘要（关键部分）：

                            OLS Regression Results
==============================================================================
Dep. Variable:          income_index   R-squared:                       0.981
Model:                            OLS   Adj. R-squared:                  0.981
Method:                 Least Squares   F-statistic:                     1816.
Date:                ...               Prob (F-statistic):           1.24e-33
Time:                        ...        Log-Likelihood:                -126.42
No. Observations:                  37   AIC:                             256.8
Durbin-Watson:                   0.1378
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         66.1491      0.778     84.997      0.000      64.569      67.729
t              4.5158      0.106     42.618      0.000       4.300       4.732
==============================================================================

DW 统计量为 0.1378，远小于 2，强烈提示残差存在正自相关。

方法二：一阶延迟因变量回归 拟合模型：$x_t = \beta_0 + \beta_1 x_{t-1} + \varepsilon_t$

# 方法二：延迟因变量模型
df['income_lag1'] = df['income_index'].shift(1)
df2 = df.dropna().copy()  # 丢弃第一个NA值
X2 = sm.add_constant(df2['income_lag1'])
model2 = sm.OLS(df2['income_index'], X2).fit()
df2['epsilon2'] = model2.resid
print("\n=== 方法二：延迟因变量模型 ===")
print(model2.summary())
# 由于包含滞后因变量，使用 Durbin h 检验
from statsmodels.stats.diagnostic import acorr_lm
# acorr_lm 可用于检验自相关，这里我们手动计算 Dh 统计量近似值
# 更严谨的做法是使用专门的函数，此处为演示逻辑
rho1 = acf(model2.resid, nlags=1, fft=False)[1]
n_obs = len(model2.resid)
var_beta1 = model2.bse['income_lag1']**2
dh_stat = rho1 * np.sqrt(n_obs / (1 - n_obs * var_beta1))
print(f"残差一阶自相关系数 ρ1: {rho1:.4f}")
print(f"Durbin h 统计量 (近似): {dh_stat:.4f}")
# 计算 p-value (双侧检验)
from scipy import stats
p_val = 2 * (1 - stats.norm.cdf(abs(dh_stat)))
print(f"Durbin h 检验 p-value: {p_val:.4f}")

运行结果摘要（关键部分）：

                            OLS Regression Results
==============================================================================
Dep. Variable:          income_index   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                 2.444e+04
Date:                ...               Prob (F-statistic):          1.21e-49
Time:                        ...        Log-Likelihood:                -122.32
No. Observations:                  36   AIC:                             248.6
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.0330      0.346      2.983      0.005       0.329       1.737
income_lag1    1.0365      0.007    156.321      0.000       1.023       1.050
==============================================================================

Durbin h 检验 p-value 远小于 0.05，同样拒绝“无自相关”的原假设。

步骤 2：残差序列的自回归模型拟合

以方法一的残差序列 epsilon1 为例，我们对其建立 AR 模型。

1. 定阶：观察残差序列的 ACF 和 PACF 图。

# 绘制方法一残差序列的 ACF 和 PACF
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(df['epsilon1'].dropna(), lags=15, ax=axes[0], title='残差序列 ACF')
plot_pacf(df['epsilon1'].dropna(), lags=15, ax=axes[1], title='残差序列 PACF')
plt.tight_layout()
plt.show()

从 PACF 图可以看到，滞后 1 阶、2 阶显著，从滞后 3 阶起截尾，提示适合 AR(2) 模型。 2. 拟合 AR(2) 模型：

from statsmodels.tsa.arima.model import ARIMA
# 对残差序列 epsilon1 拟合 AR(2) 模型
ar_model = ARIMA(df['epsilon1'].dropna(), order=(2, 0, 0), trend='n').fit()
print("\n=== 残差序列 AR(2) 模型拟合结果 ===")
print(ar_model.summary())

运行结果摘要（关键部分）：

                               SARIMAX Results
==============================================================================
Dep. Variable:                epsilon1   No. Observations:                   37
Model:                 ARIMA(2, 0, 0)   Log Likelihood                 -54.423
Date:                ...               AIC                            114.845
Time:                        ...        BIC                            119.289
Sample:                        0        HQIC                           116.393
                                - 37
Covariance Type:                  opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          1.4859      0.153      9.728      0.000       1.187       1.785
ar.L2         -0.5848      0.151     -3.864      0.000      -0.881      -0.288
sigma2         6.8893      1.594      4.321      0.000       3.764      10.014
==============================================================================

参数 $\phi_1=1.4859$, $\phi_2=-0.5848$ 均显著。模型残差 $a_t$ 应近似为白噪声，可通过 Ljung-Box 检验验证。

# 对 AR(2) 模型的残差进行白噪声检验
resid_ar = ar_model.resid
from statsmodels.stats.diagnostic import acorr_ljungbox
lb_test = acorr_ljungbox(resid_ar, lags=[6, 12], return_df=True)
print("\nAR(2) 模型残差 Ljung-Box 检验:")
print(lb_test)

如果 p 值大于 0.05，则接受残差为白噪声的原假设，说明 AR(2) 模型充分提取了信息。

步骤 3：模型比较与选择

我们将拟合的三个模型进行比较：

1. ARIMA(0,1,1) 模型（上节课结果）。
2. 残差自回归模型一：$x_t = 66.1491 + 4.5158 t + \varepsilon_t$, $\varepsilon_t = 1.4859 \varepsilon_{t-1} - 0.5848 \varepsilon_{t-2} + a_t$。
3. 残差自回归模型二：$x_t = 1.0330 + 1.0365 x_{t-1} + \varepsilon_t$, $\varepsilon_t = 0.4615 \varepsilon_{t-1} + a_t$。

使用 AIC (Akaike Information Criterion) 和 BIC (Bayesian Information Criterion) 进行模型选择，准则为：值越小，模型相对越好。

# 假设我们已经拟合了三个模型，并获取了其 AIC 和 BIC
# 这里用上节课的 ARIMA(0,1,1) 结果和本节两个残差自回归模型的结果进行对比
model_comparison = pd.DataFrame({
    'Model': ['ARIMA(0,1,1)', 'Auto-Regressive 模型一', 'Auto-Regressive 模型二'],
    'AIC': [249.3305, 260.8454, 250.6317],  # 来自课件数据
    'BIC': [252.4976, 267.2891, 253.7987]   # 来自课件数据
})
print("\n=== 模型比较 (AIC & BIC) ===")
print(model_comparison)

运行结果：

=== 模型比较 (AIC & BIC) ===
                     Model      AIC      BIC
0              ARIMA(0,1,1)  249.3305  252.4976
1  Auto-Regressive 模型一    260.8454  267.2891
2  Auto-Regressive 模型二    250.6317  253.7987

结论：根据 AIC 和 BIC 准则，ARIMA(0,1,1) 模型的拟合效果最好，其次是 Auto-Regressive 模型二，最差的是 Auto-Regressive 模型一。这表明，对于本例数据，直接使用差分处理非平稳性并拟合 ARMA 模型（即 ARIMA），比先拟合确定性趋势再对残差建模的方法更有效。这很可能是因为确定性趋势（尤其是简单的时间线性趋势）未能充分、精确地提取序列中的长期变化信息。

尽管如此，残差自回归模型（尤其是模型二）在模型可解释性上具有优势。例如，模型一可以解释为：“中国农业实际国民收入指数存在一个每年增长约 4.52 个单位的长期线性趋势，其随机波动部分具有短期记忆性（AR(2)）”。

📝 动手练一练

1. 趋势拟合方法比较：对一个有明显指数增长趋势的模拟序列，分别用“时间 $t$ 的线性模型”和“延迟一阶自回归模型”进行趋势拟合。计算两个模型的残差，并绘制其序列图。直观上，哪个模型的残差幅度更小、更随机？
2. DW 检验实践：使用 statsmodels 的 durbin_watson 函数计算上述两个趋势模型的 DW 统计量。根据结果（可查表或根据经验：DW<1.5 提示强正相关），判断哪个模型的残差自相关问题更严重？这与你从残差序列图上观察到的结果一致吗？

参考答案：

1. 对于指数增长序列，“延迟一阶自回归模型”通常能更好地捕捉其动态增长路径，因此其残差通常更小、更接近白噪声。而“时间线性模型”可能产生系统性的偏差，导致残差呈现明显的自相关模式。
2. 拟合“时间线性模型”的 DW 统计量通常会非常小（如小于 1），表明强烈的正自相关。而“延迟自回归模型”的 DW 统计量可能更接近 2（但可能因模型设定仍存在自相关）。这与残差序列图中观察到的“线性模型残差呈现连续的正负 clusters”现象是一致的。

本章小结

本节深入探讨了 残差自回归模型 (Auto-Regressive Model)，这是一种处理非平稳序列的精细建模策略。其核心在于“分而治之”：先提取可解释的确定性成分（趋势、季节），再对剩余的随机波动部分进行自回归建模，以充分挖掘数据中的自相关信息。

要点回顾：

• 动机：确定性模型拟合后，残差若存在自相关，表明信息提取不充分，需进一步建模。
• 结构：模型由“确定性部分”和“残差 AR(p) 部分”叠加而成，最终扰动项要求是白噪声。
• 趋势拟合：可采用时间函数法（直观）或延迟因变量法（动态，可能更精确）。
• 自相关检验：使用 DW 检验（通用）或 Durbin h 检验（当模型包含滞后因变量时），用于诊断一阶自相关。
• 建模流程：趋势/季节拟合 → 残差自相关检验 → 若存在自相关，则对残差序列定阶 (ACF/PACF) 并拟合 AR 模型 → 模型检验与比较。
• 模型选择：通过 AIC、BIC 等信息准则，在残差自回归模型与 ARIMA 等模型间进行客观比较。

行动清单：

1. 检验习惯：在完成任何回归或确定性模型拟合后，养成第一时间检验残差自相关性的习惯（计算 DW 统计量或绘制 ACF 图）。
2. 方法对比：面对一个具有趋势的序列，尝试分别构建 ARIMA 模型和残差自回归模型，并比较它们的 AIC/BIC 及预测效果，理解各自的适用场景。
3. 代码实现：在 Python 中熟练使用 statsmodels 库的 OLS 进行趋势拟合，用 durbin_watson 和 acorr_ljungbox 进行检验，用 ARIMA 对残差序列建模，并提取 aic、bic 属性进行模型比较。

残差自回归模型为我们提供了一种结构清晰、解释性强的建模框架。然而，它假设残差的方差是恒定的（同方差）。在金融等许多领域，时间序列的波动性（方差）本身也随时间变化，即存在“异方差”。如何对这种“波动率聚集”现象进行建模？这将是下一节 条件异方差模型 (ARCH/GARCH) 要解决的核心问题。

— 小象教研组