时间序列建模实战：从ARIMA到SARIMA，手把手预测山猫数量

1. 项目概述：从山猫数量预测看时间序列建模的实战价值

刚接触数学建模的朋友，常常会困惑于如何将课本上的理论转化为一个能解决实际问题的完整项目。我当年也是从一堆抽象的公式和算法里摸爬滚打过来的，深知“纸上得来终觉浅”的道理。今天，我们就以一个经典的入门级案例——山猫数量预测，来手把手拆解时间序列分析的全过程。这个案例之所以经典，不仅因为它数据公开、问题明确，更因为它几乎涵盖了时间序列建模的所有核心环节：数据探索、平稳性检验、模型识别、参数估计、诊断检验以及最终的预测。通过这个项目，你不仅能学会如何使用ARIMA这类经典模型，更能建立起一套处理时序数据的完整思维框架。无论你是参加数学建模竞赛的学生，还是希望用数据洞察业务趋势的从业者，这套方法都能让你在面对一串随时间变化的数据时，不再无从下手。

山猫数量数据，通常指的是加拿大哈德逊湾公司记录的1821年至1934年山猫毛皮收购量的年度时间序列。这组数据在统计学和生态学领域被广泛引用，其价值在于它呈现了一个清晰的、具有周期波动的生态种群数量变化，是练习预测的绝佳素材。我们的目标就是利用这百年的历史数据，构建一个可靠的模型，来预测未来若干年山猫的数量变化趋势。这个过程，远比单纯调用一个ARIMA()函数要丰富和深刻得多。

2. 核心思路与建模流程总览

在动手写代码之前，我们必须把整个建模的“作战地图”画清楚。时间序列分析不是一蹴而就的，它遵循一个严谨的迭代流程，通常被称为“Box-Jenkins方法”。我们的山猫预测项目将严格遵循这一流程，这能最大程度保证模型的可靠性和预测的准确性。

整个流程可以概括为四个核心阶段，它们环环相扣，上一步的输出往往是下一步的输入。第一阶段是数据准备与探索性分析。我们拿到原始的年度山猫数量数据，首先要做的就是将其导入分析环境（如Python的Pandas），并将其正确设置为时间序列索引。紧接着，不是急着建模，而是花时间“观察”数据：绘制时序图，直观感受数据的整体趋势、季节性（周期性）和波动情况；计算基本统计量，了解数据的集中和离散程度。对于山猫数据，我们预期会看到一个大约10年左右的波动周期，这是由猎物（雪鞋兔）数量周期驱动的经典生态现象。

第二阶段是序列的平稳化处理，这是时间序列建模的基石。绝大多数经典时间序列模型（如ARIMA）都要求数据是“平稳”的，即数据的统计特性（如均值、方差）不随时间推移而改变。显然，具有明显周期波动的山猫数据不满足这一点。因此，我们需要通过“差分”运算来消除趋势和周期。一阶差分可以消除线性趋势，季节性差分则可以消除固定周期的波动。我们会通过绘制差分后的序列图，并结合ADF单位根检验等统计方法，来科学判断序列是否已变得平稳。

第三阶段是模型识别与定阶。当序列平稳后，我们需要确定使用哪种模型，以及模型的参数（p, d, q）是什么。这里主要依赖两个工具：自相关函数图和偏自相关函数图。通过观察ACF和PACF图的截尾和拖尾特征，我们可以初步判断适合的模型类型（AR模型、MA模型还是ARMA模型）以及阶数的大致范围。例如，如果PACF在滞后p阶后突然截尾（落入置信区间），而ACF拖尾，则可能适合AR(p)模型。对于山猫数据，由于我们进行了差分，实际上是在构建ARIMA(p,d,q)模型，其中d就是差分的阶数。

第四阶段是模型估计、检验与预测。确定了模型形式和阶数后，我们使用最大似然估计等方法对模型参数进行估计，并得到具体的模型方程。然后，必须对模型的残差进行诊断检验：残差序列应该是白噪声（均值为零、方差恒定、无自相关）。我们可以通过绘制残差图、进行Ljung-Box检验来判断模型是否充分提取了原始序列中的信息。只有通过检验的模型，才能用于最终的预测。我们会利用拟合好的模型，对未来若干年的山猫数量进行点预测和区间预测，并直观地绘制在图表上，评估预测效果。

注意：这个流程不是线性的，而是一个循环。如果在模型检验阶段发现残差不是白噪声，说明模型拟合不充分，我们需要返回第三阶段，重新调整模型阶数或形式，直至找到一个满意的模型。这种迭代是建模工作的常态。

3. 数据准备与探索性深度解析

让我们进入实战环节。首先，我们需要获取并加载数据。山猫数量数据集在很多统计软件和开源数据包中都有收录，例如在Python的statsmodels库或R语言中都可以轻松找到。这里我们以Python环境为例。

import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.datasets import get_rdataset # 加载著名的lynx数据集（山猫） lynx = get_rdataset('lynx', 'datasets').data # 查看数据前几行 print(lynx.head()) # 通常数据是一列，列名为'value'，索引可能是数字，我们需要将其设置为时间索引 # 假设数据是1821-1934年的年度数据 years = pd.date_range(start='1821', end='1934', freq='A') # 'A'表示年末 lynx_ts = pd.Series(lynx['value'].values, index=years) lynx_ts.name = 'Lynx Trappings'

数据加载后，第一步永远是可视化。绘制时序图能给我们最直接的洞察。

plt.figure(figsize=(12, 6)) plt.plot(lynx_ts) plt.title('Annual Number of Lynx Trappings (1821-1934)') plt.xlabel('Year') plt.ylabel('Number') plt.grid(True) plt.show()

观察这张图，你可以清晰地看到几个特点：第一，序列没有明显的长期上升或下降趋势，整体围绕一个平均水平上下波动。第二，存在非常显著的周期性波动，波峰和波谷规律性地出现，周期大约在9-11年之间，这完美印证了生态学中捕食者-猎物系统的周期震荡理论。第三，波动的幅度（方差）似乎在整个时间范围内相对稳定，没有出现前期波动小、后期波动剧烈的情况，这初步暗示序列可能具有“弱平稳性”。

除了看图，我们还需要用统计量量化这些观察。计算序列的基本描述统计（均值、标准差、最小值、最大值）和分位数。更重要的是，我们可以绘制年度子图，将每一年的数据（本例中一年只有一个点，不适用）或按周期分段查看，但对于年度数据，更有效的方法是计算滚动统计量，比如10年滚动均值和滚动标准差，来观察局部趋势和波动是否稳定。

# 计算10年滚动均值和标准差 rolling_mean = lynx_ts.rolling(window=10).mean() rolling_std = lynx_ts.rolling(window=10).std() plt.figure(figsize=(12, 8)) plt.subplot(2,1,1) plt.plot(lynx_ts, label='Original') plt.plot(rolling_mean, label='10-Year Rolling Mean', color='red') plt.legend() plt.title('Original Series with Rolling Mean') plt.grid(True) plt.subplot(2,1,2) plt.plot(rolling_std, label='10-Year Rolling Std', color='green') plt.legend() plt.title('Rolling Standard Deviation') plt.grid(True) plt.tight_layout() plt.show()

如果滚动均值线大致水平，滚动标准差线也大致水平，那么序列是平稳的有力证据。对于山猫数据，滚动均值线可能会有小幅波动，但整体无趋势；滚动标准差可能会在波峰波谷处有所变化，这是周期序列的特点，我们需要通过差分来进一步处理。

实操心得：在探索性分析阶段，一定要“不厌其烦”地看图。除了时序图，还可以绘制分布直方图和Q-Q图来检查数据是否服从正态分布。许多时间序列模型假设误差项服从正态分布，虽然不是绝对必须，但了解这一点对后续模型诊断有帮助。对于山猫数据，其分布通常是有偏的（非正态），这提示我们在解释预测区间时需要谨慎。

4. 平稳性检验与差分处理实战

经过探索性分析，我们怀疑原始序列是非平稳的（因为有周期性）。现在需要用更严格的统计检验来验证，并对其进行平稳化处理。最常用的检验是增强迪基-富勒检验。

from statsmodels.tsa.stattools import adfuller result = adfuller(lynx_ts) print('ADF Statistic: %f' % result[0]) print('p-value: %f' % result[1]) print('Critical Values:') for key, value in result[4].items(): print('\t%s: %.3f' % (key, value))

ADF检验的原假设是“序列存在单位根，即非平稳”。如果p值小于显著性水平（如0.05），我们拒绝原假设，认为序列平稳。对于原始山猫序列，p值很可能大于0.05，说明它是非平稳的。

接下来进行差分处理。由于序列有周期性（约10年），我们首先考虑进行季节性差分，阶数为周期长度。由于是年度数据且周期约为10，我们尝试周期S=10。

# 季节性差分 (周期为10) lynx_diff_s = lynx_ts.diff(periods=10).dropna() # 对季节性差分后的序列再次进行ADF检验 result_diff_s = adfuller(lynx_diff_s) print('Seasonally Differenced Series p-value: ', result_diff_s[1])

如果季节性差分后序列仍然不平稳，或者为了更彻底地消除可能存在的长期依赖，我们可以在季节性差分的基础上再进行一阶普通差分。更常见的做法是，对于有明显周期的序列，直接使用非季节性差分（d=1）有时也能在一定程度上削弱周期性，我们可以比较不同差分方案的效果。

# 方案一：先季节性差分，再一阶差分（如果必要） # lynx_diff_s1 = lynx_diff_s.diff().dropna() # 方案二：直接进行一阶差分（更常用，作为ARIMA模型中的d参数） lynx_diff_1 = lynx_ts.diff().dropna() # 绘制差分后序列 fig, axes = plt.subplots(2, 1, figsize=(12, 8)) axes[0].plot(lynx_diff_s) axes[0].set_title('Seasonally Differenced (s=10) Series') axes[0].grid(True) axes[1].plot(lynx_diff_1) axes[1].set_title('First-Order Differenced Series') axes[1].grid(True) plt.tight_layout() plt.show() # 分别检验平稳性 print('ADF test for first-order differenced series p-value: ', adfuller(lynx_diff_1)[1])

通过看图，你会发现一阶差分后的序列，虽然周期性依然存在（因为一阶差分主要消除线性趋势），但序列整体围绕0值波动，ADF检验的p值通常会变得非常小（<<0.05），从统计意义上可以认为它是平稳的。在ARIMA模型中，我们正是通过设置差分阶数d=1来达成这一目的。对于季节性，我们将通过ARIMA模型中的季节性参数或后续的SARIMA模型来捕捉。

注意事项：差分是一把双刃剑。每做一次差分，就意味着我们丢失一个数据点，并且可能会过度差分，导致序列的方差变大或引入不必要的相关性。一个实用的原则是：使用能达到平稳性的最小差分阶数。通常，d=0,1,2就够了，很少需要更高阶。对于山猫数据，d=1通常是足够的。

5. 模型识别：解读ACF与PACF图的密码

获得平稳序列（这里我们以一阶差分序列lynx_diff_1为例）后，下一步就是为其选择合适的ARIMA(p,d,q)模型。其中d我们已经确定为1。现在的任务是确定自回归阶数p和移动平均阶数q。这就需要借助自相关函数图和偏自相关函数图这两把钥匙。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf fig, axes = plt.subplots(2, 1, figsize=(12, 8)) plot_acf(lynx_diff_1, lags=40, ax=axes[0]) # 查看40个滞后的自相关 axes[0].set_title('ACF Plot for First-Differenced Lynx Series') plot_pacf(lynx_diff_1, lags=40, ax=axes[1], method='ywm') # 使用Yule-Walker方法计算PACF axes[1].set_title('PACF Plot for First-Differenced Lynx Series') plt.tight_layout() plt.show()

如何解读这两张图？ACF图描述的是当前观测值与过去各期观测值之间的简单相关系数。PACF图则是在排除了中间滞后项影响后，当前观测值与过去某期观测值之间的“纯”相关系数。

观察ACF图：你会发现自相关系数在滞后10、20、30等处出现明显的峰值，并且缓慢衰减，呈现一种“拖尾”特征。这强烈暗示序列中存在季节性自相关，周期约为10。在非季节性滞后处（如滞后1、2），ACF值可能显著不为零，然后快速衰减，这提示我们可能需要非季节性的MA成分。

观察PACF图：偏自相关函数在滞后1、2处可能有显著峰值，然后迅速截尾（后面的值落在蓝色置信区间内），这提示我们可能需要一个低阶的AR成分（如p=1或2）。

基于以上观察，对于非季节性部分，一个初步的候选模型可能是ARIMA(1,1,1)、ARIMA(2,1,0)或ARIMA(0,1,2)。但是，由于明显的季节性特征，标准的ARIMA可能不够，我们需要考虑季节性ARIMA模型，即SARIMA。SARIMA模型表示为SARIMA(p,d,q)(P,D,Q,s)，其中(P,D,Q,s)是季节性部分的参数，s是周期长度（这里s=10）。

在季节性部分，观察ACF/PACF在滞后10、20处的特征：如果ACF在滞后10s处截尾，PACF拖尾，则季节性部分可能是MA(Q)；反之则可能是AR(P)。对于山猫数据，ACF在滞后10处有一个正尖峰，在滞后20处有一个负尖峰，这通常暗示需要包含季节性MA(1)成分，即Q=1。季节性差分D通常取1，以消除季节性非平稳性。

因此，一个合理的候选模型是SARIMA(1,1,1)(0,1,1,10)。这意味着：

• 非季节性部分：AR(1), 差分1阶, MA(1)
• 季节性部分：无AR，季节性差分1阶，季节性MA(1)，周期s=10。

实操心得：模型识别没有唯一正确答案，更像是一门艺术。ACF/PACF图只能给出初步指引。一个更稳健的方法是“网格搜索”：在合理的范围内（如p, q, P, Q从0到2），尝试所有可能的模型组合，然后根据信息准则（如AIC或BIC）来选择最优模型。AIC/BIC值越小，说明模型在拟合优度和复杂度之间取得了更好的平衡。我们可以在下一步模型估计中自动化这个过程。

6. 模型拟合、诊断与预测全流程

确定了候选模型结构后，我们使用统计软件来拟合模型，即估计模型中的所有参数（如AR系数、MA系数等），并对模型进行严格的诊断检验。

import warnings warnings.filterwarnings('ignore') # 忽略一些不影响结果的警告 from statsmodels.tsa.statespace.sarimax import SARIMAX import itertools # 定义参数搜索范围（为了演示，我们缩小范围。实际可扩大搜索） p = d = q = range(0, 2) # 非季节性p,d,q P = D = Q = range(0, 2) # 季节性P,D,Q s = 10 # 周期 # 生成所有参数组合 pdq = list(itertools.product(p, [1], q)) # d固定为1 seasonal_pdq = list(itertools.product(P, [1], Q, [s])) # D固定为1 best_aic = np.inf best_order = None best_seasonal_order = None print('开始网格搜索...') for param in pdq: for param_seasonal in seasonal_pdq: try: mod = SARIMAX(lynx_ts, order=param, seasonal_order=param_seasonal, enforce_stationarity=False, enforce_invertibility=False) results = mod.fit(disp=False) # disp=False不显示迭代日志 current_aic = results.aic if current_aic < best_aic: best_aic = current_aic best_order = param best_seasonal_order = param_seasonal # print(f'SARIMA{param}x{param_seasonal} - AIC:{current_aic:.2f}') except Exception as e: continue print(f'\n最优模型: SARIMA{best_order}x{best_seasonal_order}') print(f'最优AIC值: {best_aic:.2f}')

假设网格搜索得出的最优模型是SARIMA(1,1,1)(0,1,1,10)。我们用全部数据拟合这个模型。

# 拟合最优模型 best_model = SARIMAX(lynx_ts, order=best_order, # 例如 (1,1,1) seasonal_order=best_seasonal_order, # 例如 (0,1,1,10) enforce_stationarity=False, enforce_invertibility=False) best_results = best_model.fit() print(best_results.summary())

在模型摘要中，重点关注几点：1. 系数显著性：查看P>|z|列，通常小于0.05认为该系数显著不为零。2. 模型诊断：摘要底部提供了对标准化残差的一系列检验。更直观的方法是绘制诊断图。

best_results.plot_diagnostics(figsize=(12, 8)) plt.tight_layout() plt.show()

诊断图包含四个子图：

1. 标准化残差时序图：残差应该像白噪声一样随机分布在0附近，没有明显的趋势或周期。如果有，说明模型未充分提取信息。
2. 残差直方图+核密度估计与正态分布曲线对比：理想情况下，残差应近似服从正态分布。
3. 正态Q-Q图：点应大致分布在45度参考线附近，如果严重偏离，说明残差非正态。
4. 残差自相关图：所有滞后期的自相关系数都应落在置信区间内（图中蓝色区域），表明残差不存在自相关。

如果诊断图通过检验（特别是残差无自相关），说明模型是充分的。接下来就可以进行预测了。

# 预测未来20年 forecast_steps = 20 forecast_obj = best_results.get_forecast(steps=forecast_steps) forecast_mean = forecast_obj.predicted_mean forecast_ci = forecast_obj.conf_int() # 置信区间 # 创建预测时间索引 last_year = lynx_ts.index[-1] forecast_index = pd.date_range(start=last_year + pd.DateOffset(years=1), periods=forecast_steps, freq='A') # 绘制结果 plt.figure(figsize=(14, 7)) plt.plot(lynx_ts.index, lynx_ts, label='Observed (历史数据)') plt.plot(forecast_index, forecast_mean, label='Forecast', color='red') plt.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='red', alpha=0.2, label='95% Confidence Interval') plt.title('Lynx Trappings: Historical Data and 20-Year Forecast') plt.xlabel('Year') plt.ylabel('Number of Lynx Trappings') plt.legend() plt.grid(True) plt.show()

预测图会展示历史数据的拟合情况以及未来20年的预测值（红色实线）和95%的置信区间（红色阴影区域）。观察预测曲线，你应该能看到它延续了大约10年左右的周期性波动模式。

注意事项：时间序列预测的置信区间会随着预测步长的增加而迅速变宽，这意味着长期预测的不确定性非常大。对于山猫预测，20年后的预测值误差范围可能已经大到失去实际指导意义。因此，时间序列模型更擅长短期和中期预测。在报告中，务必强调这一点，避免对长期预测结果过度解读。

7. 常见问题、陷阱与实战调优技巧

在实际操作中，你几乎一定会遇到下面这些问题。这里我把自己踩过的坑和总结的技巧分享给你。

问题一：ADF检验结果与图形判断不一致怎么办？有时看图觉得序列已经平稳了，但ADF检验的p值还是大于0.05。这种情况很常见。我的建议是：图形判断优先，并结合多个检验。除了ADF，还可以做KPSS检验（原假设是平稳）。如果图形显示无明显趋势和周期，即使ADF结果稍显模糊，也可以尝试进行建模。过度差分会导致模型冗余和预测性能下降。

问题二：ACF/PACF图没有清晰的截尾或拖尾，难以定阶。这是处理真实数据，尤其是带有复杂季节性的数据时的常态。别慌，有几种策略：

1. 尝试简单的模型：从低阶开始，如ARIMA(1,1,1)，然后根据残差诊断逐步增加阶数。
2. 依赖信息准则：如前所述，使用网格搜索配合AIC/BIC来选择模型。这是更客观、自动化的方法。
3. 考虑其他模型：如果标准ARIMA/SARIMA拟合不佳，可以探索更复杂的模型，如带傅里叶项的ARIMA来处理非整数周期，或者指数平滑状态空间模型。

问题三：模型残差检验未通过，存在自相关。这说明当前模型没有完全捕捉数据中的依赖关系。你需要：

1. 增加模型阶数：回头检查ACF/PACF图，看是否在某个滞后期有显著相关被忽略了，相应增加p或q。
2. 检查是否遗漏了季节性成分：如果残差ACF在周期倍数处（如10，20）仍有峰值，说明季节性没处理好，需要调整季节性参数(P,D,Q)。
3. 添加外部变量：考虑是否有其他影响山猫数量的因素（如气候数据）可以作为外生变量加入模型（即SARIMAX模型）。

问题四：预测结果看起来“太平滑”，捕捉不到极端值。ARIMA类模型是线性模型，其预测本质上是历史数据的加权平均，因此对于极端波峰和波谷的预测往往会“收敛”向均值，显得保守。这是线性模型的固有局限。如果预测极端值至关重要，可能需要研究非线性时间序列模型，或者对数据进行变换（如对数变换）以稳定方差后再建模。

问题五：如何评估预测性能？我们不能只在训练集上自娱自乐。标准的做法是进行样本外预测评估。

1. 滚动预测：将历史数据分为训练集和测试集（如用前100年数据训练，预测后14年）。
2. 一步预测：用训练集拟合模型，预测下一期，将真实值加入训练集，重新拟合模型，再预测下一期，如此滚动进行。这模拟了实时预测场景。
3. 多步预测：直接用训练好的模型预测测试集的所有未来点。
4. 使用评估指标：计算测试集上的均方根误差、平均绝对百分比误差等，量化预测精度。RMSE对异常值敏感，MAPE是相对误差，更适合不同量级序列的比较。

# 示例：简单的样本外评估（将最后14年作为测试集） train = lynx_ts.iloc[:-14] test = lynx_ts.iloc[-14:] # 在训练集上拟合相同参数的模型 model_train = SARIMAX(train, order=best_order, seasonal_order=best_seasonal_order, enforce_stationarity=False, enforce_invertibility=False) results_train = model_train.fit(disp=False) # 预测未来14期 forecast_test = results_train.get_forecast(steps=14) predicted_values = forecast_test.predicted_mean # 计算RMSE from sklearn.metrics import mean_squared_error rmse = np.sqrt(mean_squared_error(test, predicted_values)) mape = np.mean(np.abs((test - predicted_values) / test)) * 100 print(f'测试集RMSE: {rmse:.2f}') print(f'测试集MAPE: {mape:.2f}%')

最后的建议：时间序列建模是一个需要耐心和反复迭代的过程。从山猫这个案例出发，掌握好数据探索、平稳化、模型识别、拟合诊断和评估这一套组合拳。然后，你可以将这套方法应用到任何领域的时间序列数据上，无论是股票价格、月度销售额、每日气温还是服务器流量。记住，没有一个模型是万能的，最好的模型永远是在你对业务（或问题背景）的理解与数据本身特征之间找到的那个平衡点。多练、多思考、多调参，你会越来越得心应手。