如何构建带ARIMA类误差的回归模型？求statsmodels示例代码

1. 环境准备与模拟数据 #

先导入依赖库，同时生成带自回归误差的可复现时间序列，方便你验证效果：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox

# 设置随机种子保证结果可复现
np.random.seed(42)

# 生成带AR(1)误差的时间序列：y_t = 2 + 0.3*t + e_t，其中e_t = 0.8*e_{t-1} + ε_t（ε为正态白噪声）
n = 100
t = np.arange(n)
epsilon = np.random.normal(0, 1, n)
e = np.zeros(n)
for i in range(1, n):
    e[i] = 0.8 * e[i-1] + epsilon[i]
y = 2 + 0.3*t + e

# 转为pandas Series，方便后续处理
ts = pd.Series(y, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))

# 可视化原始序列
plt.figure(figsize=(10,4))
plt.plot(ts)
plt.title('带自回归误差的时间序列')
plt.show()

2. 模型拟合 #

带自回归误差的序列，对应ARIMA模型中q=0（无移动平均项），d根据序列平稳性选择（本例序列含趋势，选择d=1做差分平稳化）。这里我们拟合ARIMA(1,1,0)（对应误差为AR(1)）：

# 拟合ARIMA模型
model = ARIMA(ts, order=(1,1,0))
results = model.fit()

# 查看模型参数与统计信息
print(results.summary())

3. 自回归阶数p的选择思路 #

如果不确定误差的自回归阶数p，可以通过两种方式判断：

• 查看差分后序列的PACF图：PACF截尾的阶数即为p的候选值

# 生成差分后序列
diff_ts = ts.diff().dropna()

# 绘制PACF图
plt.figure(figsize=(10,4))
plot_pacf(diff_ts, lags=20)
plt.title('差分后序列的PACF图')
plt.show()

• 遍历不同p值，选择AIC最小的模型：

aic_values = []
for p in range(0,4):
    try:
        model = ARIMA(ts, order=(p,1,0))
        res = model.fit()
        aic_values.append((p, res.aic))
    except:
        continue

# 筛选最优p值
best_p = min(aic_values, key=lambda x: x[1])[0]
print(f'最优自回归阶数p: {best_p}')

4. 模型诊断 #

验证残差是否为白噪声（即模型已捕捉到误差的自相关性）：

# 提取模型残差
residuals = results.resid.dropna()

# 绘制残差ACF图
plt.figure(figsize=(10,4))
plot_acf(residuals, lags=20)
plt.title('残差ACF图')
plt.show()

# Ljung-Box检验：p值>0.05则认为残差是白噪声
lb_test = acorr_ljungbox(residuals, lags=[10], return_df=True)
print('Ljung-Box检验结果：')
print(lb_test)

5. 预测未来值 #

用拟合好的模型预测后续时间点的数值：

# 预测未来10个时间点
forecast = results.get_forecast(steps=10)
forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int()

# 可视化预测结果
plt.figure(figsize=(12,5))
plt.plot(ts, label='原始序列')
plt.plot(forecast_mean, label='预测值', color='red')
plt.fill_between(forecast_ci.index, forecast_ci.iloc[:,0], forecast_ci.iloc[:,1], color='pink', alpha=0.3)
plt.title('时间序列预测结果')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Deepak Yadav