Python时间序列预测全攻略：用statsmodels搞定ARMA/ARIMA建模

Python时间序列预测全攻略：用statsmodels搞定ARMA/ARIMA建模

时间序列分析，这个听起来有些学术的词汇，实际上早已渗透到我们日常决策的方方面面。无论是金融分析师预测下一季度的股价波动，还是电商运营预估未来的商品销量，甚至是能源公司预测电网负荷，背后都离不开对历史数据规律的挖掘和未来趋势的推断。对于Python数据分析师和金融、经济领域的从业者而言，掌握一套可靠、系统的时间序列预测方法，无异于手握一把打开数据宝藏的钥匙。

在众多预测工具中，ARMA（自回归滑动平均）及其扩展模型ARIMA（差分整合移动平均自回归）家族，凭借其坚实的数学理论基础和出色的实战表现，成为了经典且核心的选择。它们不像一些“黑箱”机器学习模型那样难以解释，其模型参数具有明确的统计意义，能帮助我们深刻理解数据的内在结构——是惯性使然（自回归），还是受随机冲击影响（移动平均），亦或存在长期趋势需要处理（差分）。

然而，从理论公式到一行行可运行的代码，中间往往隔着“理解”与“实践”的鸿沟。如何判断你的数据是否适合建模？怎样从一堆ACF、PACF图中读出隐藏的模型阶数？参数(p,d,q)到底该怎么选？本指南将彻底抛开晦涩的教科书式讲解，带你一步步走完从数据准备、检验、建模到评估优化的完整闭环。我们将深度依赖Python中强大的statsmodels库，结合Jupyter Notebook的交互式环境，把每一个环节都拆解成可执行、可复现的操作。你会发现，那些看似复杂的模型，其实用几行清晰的代码就能构建起来，而其中的关键，在于掌握一套系统的方法论和解读技巧。

1. 基石：理解时间序列的平稳性与ARMA家族核心思想

在动手写代码之前，我们必须先打好地基。时间序列建模有一个几乎不可妥协的前提：平稳性。你可以把一个平稳的时间序列想象成一个在固定轨道上运行的钟摆，虽然每一刻的位置都在变化，但它围绕的中心点（均值）和摆动的幅度（方差）是稳定的，并且任意两时刻之间的关联程度只取决于它们的时间间隔，而与具体的时间起点无关。

为什么平稳性如此重要？因为ARMA/ARIMA模型的理论推导和统计性质（如参数估计、预测区间）都建立在平稳序列的假设之上。用一个预测“随机游走”的模型去预测一个存在明显上升趋势的销量数据，结果注定是失准的。常见的非平稳现象包括：

• 趋势：数据随时间呈现长期的上升或下降。
• 季节性：数据以固定周期（如每年、每季度、每月）重复出现的波动。
• 突变：序列的均值或方差在某个时间点发生结构性变化。

注意：这里说的“平稳”通常指弱平稳或宽平稳，即均值、方差恒定，且自协方差仅与时间间隔有关。完全意义上的严平稳条件更为苛刻，在实际应用中我们主要关注弱平稳。

ARMA模型正是为描述这类平稳序列而生的。它巧妙地将两种影响序列的机制结合在一起：

• AR(p) - 自回归部分：当前值y_t可以看作是过去p个历史值(y_{t-1}, y_{t-2}, ..., y_{t-p})的线性组合，再加上一个随机扰动。这体现了序列的“记忆”或“惯性”。例如，昨天的高温很可能导致今天气温也偏高。
• MA(q) - 移动平均部分：当前值y_t受到当前以及过去q个随机冲击或误差(ε_t, ε_{t-1}, ..., ε_{t-q})的影响。这体现了外部偶然事件对序列的持续影响。例如，一个突发的市场消息不仅影响当天的股价，其效应可能还会残留几天。

将两者结合，ARMA(p, q)模型便同时捕捉了序列的内在记忆和外部随机冲击。而ARIMA(p, d, q)模型则是ARMA模型的升级版，它多了一个d参数，代表差分阶数。如果原始序列不平稳，我们可以通过d次差分运算将其转换为一个平稳序列，然后再对这个差分后的序列拟合ARMA模型。简单来说，ARIMA(p,d,q)≈对原始序列进行d阶差分+对差分后序列拟合ARMA(p,q)。

理解了这些核心概念，我们就知道建模的第一步永远是：让数据变得平稳。

2. 实战第一步：数据准备与平稳性检验

让我们从一个实际场景开始。假设你手头有一份某产品过去几年的月度销售额数据。你的目标是建立模型，预测未来半年的销售情况。

2.1 数据加载与初步观察

首先，我们导入必要的库，并加载、观察数据。

import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.arima.model import ARIMA import warnings warnings.filterwarnings('ignore') # 忽略一些不影响运行的警告 # 设置绘图风格 plt.style.use('seaborn-v0_8-darkgrid') sns.set_palette("husl") # 假设数据已加载为Pandas Series，索引为时间戳 # 这里我们模拟生成一组包含趋势和季节性的非平稳数据作为示例 np.random.seed(42) time_index = pd.date_range(start='2018-01-01', periods=60, freq='M') trend = np.linspace(0, 20, 60) seasonality = 10 * np.sin(2 * np.pi * np.arange(60) / 12) noise = np.random.normal(0, 2, 60) sales_data = 100 + trend + seasonality + noise sales_series = pd.Series(sales_data, index=time_index, name='Monthly_Sales') # 绘制原始序列 fig, ax = plt.subplots(figsize=(12, 5)) sales_series.plot(ax=ax, marker='o', linewidth=2) ax.set_title('原始月度销售额序列 (含趋势与季节性)', fontsize=14, fontweight='bold') ax.set_xlabel('日期') ax.set_ylabel('销售额') plt.tight_layout() plt.show()

运行这段代码，你会看到图表清晰地显示出上升趋势和以一年为周期的季节性波动。这直观地告诉我们：原始序列不平稳。

2.2 严苛的守门人：ADF检验

视觉判断有时会主观，我们需要一个严格的统计检验——Augmented Dickey-Fuller (ADF) 检验。它的原假设（H0）是：序列存在单位根，即序列非平稳。备择假设（H1）是：序列平稳。

我们使用statsmodels的adfuller函数来进行检验。

# 执行ADF检验 adf_result = adfuller(sales_series, autolag='AIC') # autolag='AIC'自动选择最佳滞后阶数 # 整理并打印结果 adf_output = pd.Series(adf_result[0:4], index=['ADF统计量', 'p值', '使用的滞后阶数', '观测数']) for key, value in adf_result[4].items(): adf_output[f'临界值 ({key})'] = value print("ADF检验结果:") print(adf_output.to_string()) print(f"\n结论：由于p值={adf_result[1]:.4f} > 0.05，无法拒绝原假设，该序列是非平稳的。")

如果p值显著大于0.05（例如0.8），如上所述，我们接受原假设，确认序列非平稳。接下来就需要进行差分处理。

2.3 让数据“平稳”下来：差分操作

差分是消除趋势和季节性的有力工具。一阶差分是计算当前值与前一时刻值的差值：y'_t = y_t - y_{t-1}。对于季节性数据，我们还可以进行季节性差分（如对于月度数据，进行12阶差分以消除年度季节性）。

# 进行一阶差分 sales_diff_1 = sales_series.diff().dropna() # 绘制差分后的序列 fig, axes = plt.subplots(1, 2, figsize=(14, 4)) sales_diff_1.plot(ax=axes[0], marker='o', linewidth=1.5, color='coral') axes[0].set_title('一阶差分后序列', fontsize=12) axes[0].axhline(y=0, color='gray', linestyle='--', alpha=0.7) # 对差分后序列再次进行ADF检验 adf_result_diff = adfuller(sales_diff_1, autolag='AIC') axes[0].text(0.02, 0.95, f'ADF p值: {adf_result_diff[1]:.4f}', transform=axes[0].transAxes, fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8)) # 尝试季节性差分（12阶，针对月度数据） sales_diff_seasonal = sales_series.diff(12).dropna() sales_diff_seasonal.plot(ax=axes[1], marker='s', linewidth=1.5, color='teal') axes[1].set_title('季节性差分(12阶)后序列', fontsize=12) axes[1].axhline(y=0, color='gray', linestyle='--', alpha=0.7) adf_result_seasonal = adfuller(sales_diff_seasonal, autolag='AIC') axes[1].text(0.02, 0.95, f'ADF p值: {adf_result_seasonal[1]:.4f}', transform=axes[1].transAxes, fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8)) plt.tight_layout() plt.show()

通过对比差分后序列的图表和ADF检验的p值，我们可以判断差分是否有效。通常，p值小于0.05即认为差分后序列已变得平稳。此时，我们便确定了ARIMA模型中的d（差分阶数）。例如，如果一阶差分就使序列平稳，则d=1；如果需要先做季节性差分再做一阶差分才平稳，则可能需要考虑季节性ARIMA模型（SARIMA），这更复杂一些。

3. 破解模型密码：ACF与PACF图解读技巧

确定了d之后，接下来要确定AR部分的阶数p和MA部分的阶数q。这里，自相关函数（ACF）图和偏自相关函数（PACF）图是我们的“罗塞塔石碑”。

• ACF图：展示时间序列y_t与其自身滞后k期值y_{t-k}之间的总体相关性。这种相关性包含了所有中间滞后项（y_{t-1}, ..., y_{t-k+1}）的间接影响。
• PACF图：展示在控制了中间滞后项（y_{t-1}, ..., y_{t-k+1}）的影响后，y_t与y_{t-k}之间的“纯粹”相关性。

对于平稳序列，识别p和q的经典法则如下表所示：

提示：“截断”并非指在某个阶数后ACF/PACF值绝对为0，而是在该阶数之后，大部分值落入置信区间（图中蓝色阴影区域）内，且没有明显的规律性。“拖尾”则是指以指数衰减或正弦波形式缓慢衰减，不会突然切断。

让我们绘制经过一阶差分后（假设已平稳）序列的ACF和PACF图。

# 假设 sales_diff_1 是平稳的 fig, axes = plt.subplots(1, 2, figsize=(14, 4)) # 绘制ACF图 plot_acf(sales_diff_1, lags=30, ax=axes[0], title='差分后序列的自相关函数 (ACF)') axes[0].set_xlabel('滞后阶数') # 绘制PACF图 plot_pacf(sales_diff_1, lags=30, ax=axes[1], title='差分后序列的偏自相关函数 (PACF)', method='ywm') axes[1].set_xlabel('滞后阶数') plt.tight_layout() plt.show()

现在，你需要像侦探一样审视这两张图：

1. 看ACF图：它是否在某个滞后阶数q之后迅速衰减并保持在置信区间内？如果是，可能暗示MA(q)。
2. 看PACF图：它是否在某个滞后阶数p之后迅速衰减并保持在置信区间内？如果是，可能暗示AR(p)。
3. 如果两者都拖尾，则可能是ARMA(p,q)，需要进一步确定p和q。

在实际数据中，图形往往不会像教科书那么完美。可能存在季节性相关（在滞后12、24等位置出现峰值），或者拖尾/截断的界限模糊。这时，ACF/PACF图主要提供初步的线索，我们需要结合更客观的信息准则来最终定夺。

4. 构建与评估：statsmodels建模全流程

有了对p, d, q的初步猜测，我们就可以开始构建模型了。statsmodels提供了ARIMA类（注意：在较新版本中，ARIMA是ARIMA模型的接口，它整合了之前的ARMA和ARIMA模型），使用起来非常直观。

4.1 模型拟合与参数估计

假设我们初步判断d=1，并从ACF/PACF图中猜测p=1,q=1，即尝试ARIMA(1,1,1)模型。

# 拟合 ARIMA(1,1,1) 模型 model_111 = ARIMA(sales_series, order=(1, 1, 1)) # order参数为 (p, d, q) model_fit_111 = model_111.fit() # 打印详细的模型摘要 print(model_fit_111.summary())

模型摘要会输出海量信息，我们需要重点关注以下几部分：

• 系数表 (coef)：显示ar.L1（AR系数）、ma.L1（MA系数）的估计值及其标准误、z统计量和p值。p值（P>|z|）小于0.05通常认为该系数显著不为零，应保留在模型中。
• 信息准则：AIC(Akaike Information Criterion) 和BIC(Bayesian Information Criterion)。它们的值越小，通常表示模型在拟合优度和复杂度之间取得了更好的平衡。用于比较不同模型的优劣。
• 残差检验：Ljung-Box (Q)统计量及其p值。用于检验残差是否为白噪声（即是否还有未被模型提取的信息）。理想的残差应该是白噪声，此时Q检验的p值应大于0.05。

4.2 模型诊断：残差分析

一个“好”的模型，其残差（实际值与模型预测值之差）应该看起来像随机的白噪声，没有自相关，也符合正态分布。我们可以对拟合模型的残差进行系统诊断。

# 获取模型残差 residuals = model_fit_111.resid # 绘制残差诊断图 fig, axes = plt.subplots(2, 2, figsize=(13, 8)) fig.suptitle('ARIMA(1,1,1) 模型残差诊断', fontsize=16, fontweight='bold') # 1. 残差序列图 axes[0, 0].plot(residuals, marker='o', markersize=4, linewidth=0.8) axes[0, 0].axhline(y=0, color='r', linestyle='--', alpha=0.7) axes[0, 0].set_title('残差序列') axes[0, 0].set_xlabel('时间') axes[0, 0].set_ylabel('残差') # 2. 残差直方图与正态分布曲线 from scipy.stats import norm axes[0, 1].hist(residuals, bins=15, density=True, alpha=0.7, edgecolor='black', color='skyblue') xmin, xmax = axes[0, 1].get_xlim() x = np.linspace(xmin, xmax, 100) p = norm.pdf(x, residuals.mean(), residuals.std()) axes[0, 1].plot(x, p, 'k', linewidth=2, label='正态分布') axes[0, 1].set_title('残差分布 (vs 正态分布)') axes[0, 1].legend() # 3. 残差ACF图 plot_acf(residuals, lags=30, ax=axes[1, 0], title='残差自相关函数') axes[1, 0].set_xlabel('滞后阶数') # 4. Q-Q图 (分位数-分位数图) from statsmodels.graphics.gofplots import qqplot qqplot(residuals, line='s', ax=axes[1, 1], marker='o', markersize=4, alpha=0.7) axes[1, 1].set_title('残差Q-Q图') plt.tight_layout(rect=[0, 0, 1, 0.96]) # 调整布局给总标题留空间 plt.show() # 进行Ljung-Box检验（statsmodels摘要中已包含，这里可单独再检验更长期滞后） from statsmodels.stats.diagnostic import acorr_ljungbox lb_test = acorr_ljungbox(residuals, lags=[10, 20], return_df=True) # 检验滞后10和20阶 print("\n残差Ljung-Box检验 (H0: 残差是白噪声):") print(lb_test)

理想的诊断结果应该是：

• 残差序列图：围绕0随机波动，无明显趋势或周期性。
• 残差ACF图：所有滞后阶数的自相关系数均落在置信区间内，无显著相关。
• Q-Q图：散点大致分布在红色参考线（理论正态分布线）附近。
• Ljung-Box检验：p值大于0.05，不拒绝“残差是白噪声”的原假设。

如果诊断结果不理想（如残差ACF仍有显著峰值，或Q-Q图偏离严重），说明当前(p,d,q)组合的模型未能充分捕捉数据特征，需要调整模型。

4.3 自动化寻优：网格搜索与信息准则

手动尝试不同的(p,d,q)组合效率低下。我们可以通过网格搜索（Grid Search），遍历一个合理的参数空间，并选择AIC或BIC最小的模型作为最优模型。这是一种基于信息准则的自动化模型选择方法。

# 定义参数搜索范围 p_range = range(0, 4) # AR阶数尝试0到3 d_range = [1] # 差分阶数，根据前期分析已确定为1 q_range = range(0, 4) # MA阶数尝试0到3 best_aic = np.inf best_bic = np.inf best_order = None best_model = None results_list = [] # 网格搜索 for p in p_range: for d in d_range: for q in q_range: try: model = ARIMA(sales_series, order=(p, d, q)) model_fit = model.fit() current_aic = model_fit.aic current_bic = model_fit.bic results_list.append({'order': (p,d,q), 'aic': current_aic, 'bic': current_bic}) # 寻找AIC最小的模型 if current_aic < best_aic: best_aic = current_aic best_order_aic = (p, d, q) best_model_aic = model_fit # 寻找BIC最小的模型（BIC对模型复杂度惩罚更重） if current_bic < best_bic: best_bic = current_bic best_order_bic = (p, d, q) best_model_bic = model_fit except Exception as e: # 某些参数组合可能导致估计失败，跳过 print(f"参数组合 ({p},{d},{q}) 拟合失败: {e}") continue # 将结果转为DataFrame便于查看 results_df = pd.DataFrame(results_list).sort_values('aic') print("所有尝试模型的信息准则排序 (按AIC升序):") print(results_df.to_string(index=False)) print(f"\n基于AIC的最优模型阶数: {best_order_aic}, AIC = {best_aic:.2f}") print(f"基于BIC的最优模型阶数: {best_order_bic}, BIC = {best_bic:.2f}") # 通常我们选择AIC或BIC更小的模型，这里以AIC为例 final_model = best_model_aic print(f"\n最终选定模型: ARIMA{best_order_aic}") print(final_model.summary())

网格搜索会输出一个按AIC排序的模型列表。AIC和BIC都权衡了模型的拟合优度与参数数量，倾向于选择“足够好但又简洁”的模型。BIC对参数数量的惩罚比AIC更重，因此在样本量较大时，BIC倾向于选择更简单的模型。你可以根据实际情况选择AIC或BIC准则。

5. 预测未来与模型部署

模型通过诊断检验并被选定后，就可以用来进行预测了。statsmodels的get_forecast和get_prediction方法非常方便。

5.1 样本内拟合与样本外预测

# 获取样本内拟合值（历史数据的拟合） fitted_values = final_model.fittedvalues # 预测未来12个时间点（步长） forecast_steps = 12 forecast_result = final_model.get_forecast(steps=forecast_steps) forecast_mean = forecast_result.predicted_mean # 点预测值 forecast_ci = forecast_result.conf_int(alpha=0.05) # 95%置信区间 # 创建预测时间索引 last_date = sales_series.index[-1] if isinstance(last_date, pd.Timestamp): future_index = pd.date_range(start=last_date + pd.tseries.offsets.MonthEnd(1), periods=forecast_steps, freq='M') else: future_index = range(len(sales_series), len(sales_series) + forecast_steps) # 绘制历史数据、拟合值及未来预测 fig, ax = plt.subplots(figsize=(14, 7)) ax.plot(sales_series.index, sales_series, label='历史实际值', marker='o', linewidth=2) ax.plot(sales_series.index[1:], fitted_values[1:], label='样本内拟合值', linestyle='--', linewidth=2, color='orange') # 注意差分后拟合值起始点 ax.plot(future_index, forecast_mean, label='未来预测值', marker='s', linewidth=2, color='red') ax.fill_between(future_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='red', alpha=0.2, label='95% 预测区间') ax.set_title('ARIMA模型拟合与预测', fontsize=16, fontweight='bold') ax.set_xlabel('时间') ax.set_ylabel('销售额') ax.legend(loc='best') ax.grid(True, alpha=0.3) plt.tight_layout() plt.show() # 打印预测结果表格 forecast_df = pd.DataFrame({ '预测日期': future_index, '点预测值': forecast_mean.values, '预测下限(95%)': forecast_ci.iloc[:, 0].values, '预测上限(95%)': forecast_ci.iloc[:, 1].values }) print("未来12期预测结果:") print(forecast_df.to_string(index=False))

图表将清晰地展示模型对历史数据的拟合情况以及对未来趋势的预测，并给出预测的不确定性范围（置信区间）。预测区间随着预测步长的增加而变宽，这符合我们对未来不确定性增加的认知。

5.2 模型评估与持续迭代

在将模型投入实际生产环境前，还需要进行严格的样本外评估。一种常见的方法是时间序列交叉验证或滚动预测。

# 简单的滚动预测评估示例（训练集/测试集分割） train_size = int(len(sales_series) * 0.8) train, test = sales_series[:train_size], sales_series[train_size:] predictions = [] history = list(train) for t in range(len(test)): model = ARIMA(history, order=best_order_aic) # 使用之前网格搜索得到的最优阶数 model_fit = model.fit() yhat = model_fit.forecast()[0] # 预测下一步 predictions.append(yhat) history.append(test.iloc[t]) # 将真实值加入历史，模拟实时更新 # 计算评估指标 from sklearn.metrics import mean_absolute_error, mean_squared_error mae = mean_absolute_error(test, predictions) rmse = np.sqrt(mean_squared_error(test, predictions)) mape = np.mean(np.abs((test.values - np.array(predictions)) / test.values)) * 100 print(f"测试集评估指标 (滚动预测):") print(f" 平均绝对误差 (MAE): {mae:.2f}") print(f" 均方根误差 (RMSE): {rmse:.2f}") print(f" 平均绝对百分比误差 (MAPE): {mape:.2f}%")

通过计算MAE、RMSE、MAPE等指标，可以量化模型的预测精度。如果指标不理想，可能需要回到之前的步骤：重新检查数据平稳性、尝试更复杂的模型（如SARIMA处理季节性）、或引入外部变量（走向ARIMAX等模型）。

最后，将训练好的最终模型（final_model）用pickle或joblib库保存下来，就可以集成到你的数据分析管道或应用系统中，定期用新数据更新模型参数，实现持续、自动化的时间序列预测。记住，没有一劳永逸的模型，定期用新数据重新评估和调整是保证预测效果的关键。