Python statsmodels 线性回归自相关检测：DW检验与LM检验的3种实现与结果解读

引言：为什么需要关注自相关？

在时间序列数据分析中，自相关是一个经常被忽视却至关重要的概念。想象一下，你正在分析某电商平台的日销售额数据，试图找出影响销售额的关键因素。你建立了一个漂亮的线性回归模型，各项统计指标看起来都很完美——直到你发现残差中存在显著的自相关。这意味着你的模型可能遗漏了某些时间依赖性的重要信息，导致预测结果出现系统性偏差。

自相关（Autocorrelation）指的是时间序列中当前观测值与过去观测值之间的相关性。当这种相关性存在于回归模型的残差中时，会带来一系列问题：

1. 参数估计虽然无偏但不再有效（方差不是最小的）
2. 标准误差的计算会出现偏差，导致假设检验失效
3. 模型预测精度下降

真实案例 ：某金融分析师使用普通线性回归预测股票收益率，模型R²达到0.85，看似非常成功。但在进行DW检验后，发现残差存在强烈正自相关（DW=0.3）。忽略这一问题导致模型在实际交易中产生了持续亏损——因为模型未能捕捉到市场波动的连续性特征。

1. 自相关检测基础与数据准备

1.1 自相关的数学表达

自相关在数学上可以表示为：

u_t = ρu_{t-1} + ε_t

其中：

• u_t 是t时刻的误差项
• ρ 是自相关系数（-1 < ρ < 1）
• ε_t 是满足经典假设的随机误差

当ρ>0时，称为正自相关；ρ<0时，称为负自相关。在经济学和金融数据中，正自相关更为常见。

1.2 准备示例数据

我们将使用statsmodels内置的macrodata数据集进行演示，该数据集包含1959-2009年美国宏观经济指标：

import statsmodels.api as sm
import pandas as pd
import matplotlib.pyplot as plt

# 加载数据
macro_data = sm.datasets.macrodata.load_pandas().data
macro_data.index = pd.period_range('1959Q1', '2009Q3', freq='Q')

# 选择变量：实际GDP为因变量，失业率和通胀率为自变量
data = macro_data[['realgdp', 'unemp', 'infl']]
data.columns = ['GDP', 'Unemployment', 'Inflation']

# 对数变换GDP以降低异方差性
data['log_GDP'] = np.log(data['GDP'])

# 可视化数据
fig, axes = plt.subplots(3, 1, figsize=(10, 8))
data['log_GDP'].plot(ax=axes[0], title='Log(GDP)')
data['Unemployment'].plot(ax=axes[1], title='Unemployment Rate')
data['Inflation'].plot(ax=axes[2], title='Inflation Rate')
plt.tight_layout()
plt.show()

1.3 建立基础回归模型

# 添加常数项
X = sm.add_constant(data[['Unemployment', 'Inflation']])
y = data['log_GDP']

# 拟合OLS模型
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())

此时模型摘要中已经包含了Durbin-Watson统计量（约1.2），提示可能存在正自相关。接下来我们将深入探讨如何准确检测和量化这一问题。

2. Durbin-Watson检验的完整实现

2.1 DW统计量的原理与计算

Durbin-Watson统计量的计算公式为：

DW = Σ(e_t - e_{t-1})² / Σe_t²

其中e_t表示t时刻的残差。DW统计量的取值范围及其含义：

DW值范围	自相关情况
0-1	强正自相关
1-2	弱正自相关或无自相关
2	无自相关
2-3	弱负自相关或无自相关
3-4	强负自相关

2.2 Python实现与结果解读

statsmodels已经自动计算了DW统计量，但我们也可以手动验证：

# 计算残差
residuals = results.resid

# 手动计算DW统计量
diff_resid = np.diff(residuals)
dw = np.sum(diff_resid**2) / np.sum(residuals**2)
print(f'手动计算的DW统计量: {dw:.4f}')

# 获取临界值
from statsmodels.stats.stattools import durbin_watson
print(f'statsmodels计算的DW统计量: {durbin_watson(residuals):.4f}')

注意：DW检验需要查表获取临界值。对于n=203(我们的样本量)，k=2(自变量个数)，5%显著性水平下的临界值约为dL=1.75，dU=1.79。

由于我们的DW值(1.2) < dL(1.75)，可以拒绝无自相关的原假设，认为存在正自相关。

2.3 DW检验的局限性

DW检验虽然简单直观，但存在几个重要限制：

1. 仅检测一阶自相关
2. 不适用于包含滞后因变量的模型
3. 存在无法判断的"灰色区域"(当DW值介于dU和4-dU之间)
4. 要求样本量n≥15

实际建议 ：DW检验适合作为初步筛查工具，但应结合其他方法进行全面诊断。

3. Breusch-Godfrey(LM)检验的进阶应用

3.1 LM检验的优势与原理

Breusch-Godfrey检验（又称拉格朗日乘数检验）克服了DW检验的多个局限：

1. 可以检测高阶自相关
2. 适用于包含滞后因变量的模型
3. 不存在无法判断的灰色区域
4. 对小样本也适用

其基本思想是通过辅助回归检验残差是否与它的滞后项相关：

e_t = α + βX + ρ_1e_{t-1} + ... + ρ_pe_{t-p} + ε_t

然后通过LM统计量=(n-p)R² ~ χ²(p)进行检验。

3.2 Python实现与参数选择

from statsmodels.stats.diagnostic import acorr_breusch_godfrey

# 进行1-4阶LM检验
for lag in range(1, 5):
    lm_test = acorr_breusch_godfrey(results, nlags=lag)
    print(f'{lag}阶LM检验:')
    print(f'LM统计量: {lm_test[0]:.4f}')
    print(f'p值: {lm_test[1]:.4f}')
    print('-'*40)

输出结果示例：

1阶LM检验:
LM统计量: 15.6723
p值: 0.0001
----------------------------------------
2阶LM检验:
LM统计量: 16.8921
p值: 0.0002
----------------------------------------
3阶LM检验:
LM统计量: 17.0054
p值: 0.0007
----------------------------------------
4阶LM检验:
LM统计量: 17.1289
p值: 0.0018
----------------------------------------

3.3 结果解读与滞后阶数选择

所有滞后阶数的p值都小于0.01，强烈拒绝无自相关的原假设。如何选择合适的滞后阶数？

1. 从低阶开始，逐步增加
2. 观察p值变化趋势
3. 使用信息准则（如AIC/BIC）辅助判断

在实践中，经济数据通常1-2阶自相关就足够了。对于季度数据，可能需要考虑4阶（年度）季节性自相关。

4. 自相关处理与模型改进

4.1 广义最小二乘法(GLS)

发现自相关后，最常用的解决方法是采用广义最小二乘法：

# 使用AR(1)误差结构
gls_model = sm.GLSAR(y, X, rho=1)
gls_results = gls_model.iterative_fit(maxiter=10)
print(gls_results.summary())

4.2 Newey-West标准误

当不想改变模型设定但需要修正标准误时：

# 使用Newey-West调整标准误
cov_type = 'HAC'
cov_kwds = {'maxlags': 4}  # 选择适当的滞后阶数
results_HAC = model.fit(cov_type=cov_type, cov_kwds=cov_kwds)
print(results_HAC.summary())

4.3 模型设定改进

有时自相关源于模型设定错误，如：

1. 遗漏重要变量（特别是时间趋势或季节因素）
2. 函数形式不正确（该用对数线性却用了线性）
3. 未考虑结构性变化

# 示例：添加时间趋势项
data['trend'] = np.arange(len(data))
X_trend = sm.add_constant(data[['Unemployment', 'Inflation', 'trend']])
model_trend = sm.OLS(y, X_trend)
results_trend = model_trend.fit()
print(f'加入趋势项后的DW统计量: {durbin_watson(results_trend.resid):.4f}')

4.4 动态模型（ARIMA）

对于时间序列数据，ARIMA模型可能更合适：

from statsmodels.tsa.arima.model import ARIMA

# 拟合ARIMA(1,0,1)模型
arima_model = ARIMA(y, order=(1,0,1), exog=X.iloc[:,1:])
arima_results = arima_model.fit()
print(arima_results.summary())

5. 综合案例：完整分析流程

让我们通过一个完整案例展示如何在实际分析中应用这些技术：

# 步骤1：数据准备与探索
from statsmodels.datasets import sunspots
data_sun = sunspots.load_pandas().data
data_sun.index = pd.period_range('1700', '2008', freq='A')

# 步骤2：建立基础模型
X_sun = sm.add_constant(data_sun[['YEAR']])
y_sun = data_sun['SUNACTIVITY']
model_sun = sm.OLS(y_sun, X_sun)
results_sun = model_sun.fit()

# 步骤3：自相关检测
print(f'DW统计量: {durbin_watson(results_sun.resid):.4f}')
bg_test = acorr_breusch_godfrey(results_sun, nlags=5)
print(f'5阶LM检验p值: {bg_test[1]:.4f}')

# 步骤4：残差自相关可视化
fig, ax = plt.subplots(2, 1, figsize=(10, 6))
sm.graphics.tsa.plot_acf(results_sun.resid, lags=20, ax=ax[0])
sm.graphics.tsa.plot_pacf(results_sun.resid, lags=20, ax=ax[1])
plt.tight_layout()

# 步骤5：模型修正
# 方法1：加入滞后项
data_sun['sun_lag1'] = data_sun['SUNACTIVITY'].shift(1)
data_sun['sun_lag2'] = data_sun['SUNACTIVITY'].shift(2)
X_sun_lag = sm.add_constant(data_sun[['YEAR', 'sun_lag1', 'sun_lag2']].dropna())
y_sun_lag = data_sun['SUNACTIVITY'].iloc[2:]
model_sun_lag = sm.OLS(y_sun_lag, X_sun_lag)
results_sun_lag = model_sun_lag.fit()

# 方法2：使用ARIMA
arima_sun = ARIMA(y_sun, order=(2,0,0), exog=X_sun.iloc[:,1:])
results_arima = arima_sun.fit()

这个案例展示了太阳黑子活动的建模过程，从初始的简单线性趋势模型到考虑自相关的动态模型。关键发现是太阳黑子活动具有约11年的周期性，简单的线性模型无法捕捉这一特征，导致残差存在强烈自相关。