时间序列建模第一步：如何识别系统类型，从噪声干扰到模型选择常见错误与实例解析
先说说，这事儿到底难在哪
我第一次接触时间序列建模的时候，满脑子都是公式——ARIMA、VAR、状态空间模型，一个比一个高大上。结果真拿到数据一看，傻眼了。这数据忽高忽低，跟心电图似的，到底是噪声还是真信号？是自相关的还是外生因素驱动的？选错模型，后面的分析全白搭。
很多人第一步就踩坑了：上来直接跑模型，不看数据、不识别系统类型，最后结果乱七八糟，还找不到原因。今天咱们就把这套流程掰开揉碎讲清楚，让你真正搞明白怎么从一堆噪声里找出规律。
时间序列到底长什么样？先别急着建模
在碰任何模型之前，你得先学会”看”数据。这不是让你随便画个图就完事，而是要带着问题去看。
一个真实的故事
我有个朋友叫小陈，做零售销量预测的。有次他拿到一家超市三年多的日销量数据，激动得睡不着觉，第二天就跑来找我：”我搞了个ARIMA模型，效果还行！”
我看了眼他的结果，忍不住问：”你数据里有没有节假日效应？有没有促销活动？”他愣了一下说：”没，我就直接拿模型跑了。”
结果呢？模型在节假日前后预测完全崩了，因为他根本没识别出数据的结构。这才是典型的”模型先行，分析后置”，是新手最容易犯的错误。
时间序列数据本质上是在时间维度上的有序观测，每个观测值可能受到前面多个值的影响，也可能受到外部因素的驱动。你要做的第一件事，不是选模型，而是搞清楚这数据背后到底有什么在”操作”。
第一步：先看原始数据图，别嫌它土
原始时序图能告诉你什么
把数据画出来，横轴是时间，纵轴是数值，就这么简单。别小看这张图，很多关键信息都在里面。
我们来看一个具体例子。假设你手里有一份某城市2019年到2023年的日均气温数据（单位：摄氏度）：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import seaborn as sns
# 模拟一份气温数据（带趋势、季节性和噪声）
np.random.seed(42)
dates = pd.date_range('2019-01-01', '2023-12-31', freq='D')
n = len(dates)
# 趋势成分（缓慢上升，可能是城市热岛效应）
trend = 0.002 * np.arange(n)
# 季节性成分（年周期，约365天）
seasonality = 15 * np.sin(2 * np.pi * np.arange(n) / 365.25) + \
5 * np.sin(2 * np.pi * 2 * np.arange(n) / 365.25)
# 噪声（随机波动）
noise = np.random.normal(0, 3, n)
# 合成数据
temperature = 20 + trend + seasonality + noise
df_temp = pd.DataFrame({'date': dates, 'temperature': temperature})
画出原始数据图：
plt.figure(figsize=(14, 6))
plt.plot(df_temp['date'], df_temp['temperature'], color='steelblue', linewidth=0.8, alpha=0.8)
plt.title('Daily Average Temperature (2019-2023)', fontsize=14)
plt.xlabel('Date', fontsize=12)
plt.ylabel('Temperature (°C)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('raw_time_series.png', dpi=150)
plt.show()
从这张图里，你应该能读出几个东西：
趋势（Trend）：数据整体是往上升的，说明存在一个长期走向。小陈的气销数据里如果也有这种情况，说明零售规模在扩大。
季节性（Seasonality）：气温数据有明显的 yearly 周期，夏天高冬天低。如果你的销售数据也有 yearly 或 monthly 周期，那说明有季节性成分。
波动幅度（Volatility）：数据围绕某个中心上下波动，这些波动一部分是规律性的，一部分是随机的噪声。
突变点（Structural Breaks）：某些时间点数据突然变了，比如政策出台、自然灾害、疫情等。
小朋友也能看懂的比喻：想象你在看一条河的水位。趋势就像这条河整体是在变深还是变浅；季节性就像每天潮汐涨落，有规律地来去；噪声就像水面上随机的波纹；突变点就像突然下了一场大雨，水位一下子跳上去了。建模前，你得先搞清楚这条河是什么样的。
第二步：分解数据，把”混在一起的信号”拆开
什么是时间序列分解？
刚才那张原始图里，趋势、季节性和噪声混在一起，你看不清楚各自长啥样。这时候就需要做时间序列分解，把数据拆成几个组成部分。
from statsmodels.tsa.seasonal import seasonal_decompose
# 需要做季节分解，先用 yearly 周期
result = seasonal_decompose(df_temp.set_index('date')['temperature'],
period=365, model='additive')
# 画出分解结果
fig, axes = plt.subplots(4, 1, figsize=(14, 10))
axes[0].plot(result.observed)
axes[0].set_title('Observed (原始数据)', fontsize=12)
axes[1].plot(result.trend)
axes[1].set_title('Trend (趋势成分)', fontsize=12)
axes[2].plot(result.seasonal)
axes[2].set_title('Seasonal (季节成分)', fontsize=12)
axes[3].plot(result.resid)
axes[3].set_title('Residual (残差/噪声)', fontsize=12)
for ax in axes:
ax.tick_params(labelsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('decomposition.png', dpi=150)
plt.show()
分解之后，你就能清楚地看到：
- Observed：原始数据
- Trend：长期走向，这条线应该相对平滑
- Seasonal：周期性波动的部分，每年重复的模式
- Residual：去掉趋势和季节性之后剩下的部分
关键问题来了：残差是什么？
残差（Residual）是时间序列分析里最重要的概念之一。如果残差看起来是纯粹的随机噪声（白噪声），那说明你的模型已经把数据里可预测的部分都抓出来了。如果残差里还有规律，说明你漏掉了什么。
判断残差好坏的方法：你看残差图，如果它像心电图一样在零附近随机跳动，没有明显的模式，那就还不错。如果残差里还看得出趋势或周期性，那就说明分解不够彻底，或者该部分还没被模型捕捉。
第三步：看自相关图，这是识别系统类型的”核心武器”
什么是自相关（ACF）？
自相关就是同一个时间序列在不同时间滞后下的相关性。比如今天的温度和昨天的温度相关度高不高？和7天前的相关度高不高？和365天前的相关度高不高？
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
# ACF图
plot_acf(df_temp.set_index('date')['temperature'].dropna(), lags=400, ax=axes[0], alpha=0.05)
axes[0].set_title('ACF (自相关函数)', fontsize=12)
# PACF图
plot_pacf(df_temp.set_index('date')['temperature'].dropna(), lags=400, ax=axes[1],
method='ywm', alpha=0.05)
axes[1].set_title('PACF (偏自相关函数)', fontsize=12)
plt.tight_layout()
plt.savefig('acf_pacf.png', dpi=150)
plt.show()
这里出现了两个关键工具：ACF（自相关函数） 和 PACF（偏自相关函数）。
怎么读ACF和PACF图？
这是很多初学者的”噩梦”，因为那两根蓝色的阴影带子看着密密麻麻的柱子让人头皮发麻。但其实读法很简单：
- 横轴：滞后阶数（lag），也就是你拿当前数据和过去第几期比较
- 纵轴：相关系数，取值在 -1 到 1 之间
- 蓝色阴影区域：这是”不显著”的区间。柱子如果在阴影里，说明这个滞后的相关性不显著，大概率是噪声
- 柱子伸出阴影：说明这个滞后的相关性是显著的，是真实的信号
拿气温数据来说
看ACF图：
- 在 lag=365 附近有一个高峰——这说明年周期性很强，365天前的温度和今天高度相关
- 在 lag=182 附近也有一个小峰——这是半年周期的体现
- 随着滞后增加，柱子逐渐衰减——这说明趋势成分被处理后，数据趋于平稳
看PACF图：
- lag=1 的柱子特别高，然后快速衰减——这说明一阶自相关很强
- 在 lag=365 处也有一个明显的柱子——对应年季节性
给小朋友的解释：ACF就像你在问”今天的我和昨天的我像不像？像前天的我像不像？像一年前的我像不像？” PACF则更严格一点，它问的是”今天的我和昨天的我像不像？但在排除了今天和前天、大前天的关系之后，昨天还跟我像不像？”
第四步：判断平稳性，这是选模型的分水岭
什么是平稳性？
时间序列分析里有个最重要的前提：数据最好是平稳的。什么叫平稳？
- 弱平稳：数据的均值、方差和自协方差不随时间变化
- 说白了，就是数据不会”跑偏”，不会”越变越大”，也不会”规律一直在变”
怎么检验平稳性？
方法一：看数据图
刚才分解图里的 Trend 如果不是水平线，而是有上升或下降趋势，那就说明数据不平稳。
方法二：用统计检验——ADF检验
from statsmodels.tsa.stattools import adfuller
def check_stationarity(series, name='Series'):
result = adfuller(series.dropna())
print(f'=== {name} 的ADF检验结果 ===')
print(f'ADF统计量: {result[0]:.4f}')
print(f'p值: {result[1]:.4f}')
print(f'临界值(1%): {result[4]["1%"]:.4f}')
print(f'临界值(5%): {result[4]["5%"]:.4f}')
print(f'临界值(10%): {result[4]["10%"]:.4f}')
if result[1] < 0.05:
print(f'结论: p值<{0.05}，数据是平稳的 ✓')
else:
print(f'结论: p值>{0.05}，数据不平稳，需要差分处理 ✗')
print()
check_stationarity(df_temp.set_index('date')['temperature'], '原始气温数据')
# 对原始数据做一阶差分后再检验
df_temp['temp_diff1'] = df_temp['temperature'].diff()
check_stationarity(df_temp['temp_diff1'], '一阶差分后')
输出大概长这样：
=== 原始气温数据 的ADF检验结果 ===
ADF统计量: -2.1543
p值: 0.2241
临界值(1%): -3.4312
临界值(5%): -2.8618
临界值(10%): -2.5668
结论: p值>0.05，数据不平稳，需要差分处理 ✗
=== 一阶差分后 的ADF检验结果 ===
ADF统计量: -12.8734
p值: 0.0000
临界值(1%): -3.4315
临界值(5%): -2.8619
临界值(10%): -2.5669
结论: p值<0.05，数据是平稳的 ✓
差分的本质是什么？
差分就是算相邻两期的差值：
\[y'_t = y_t - y_{t-1}\]
一阶差分消除了线性趋势，二阶差分消除了二次趋势。大多数经济和时间序列数据经过1到2阶差分后就能变成平稳序列。
直观理解：想象你爬一座山。你现在的位置（原始数据）一直在升高，但如果你看”每一步走了多高”（差分），那就是一个相对稳定的数值了。爬山的高度在变，但每一步的高度变化可能比较稳定。
第五步：识别系统类型——这是建模的”导航图”
时间序列系统到底有哪几种？
从建模的角度，时间序列数据可以大致分成以下几类：
| 系统类型 | 特征 | 典型例子 |
|---|---|---|
| 白噪声系统 | 没有任何自相关，完全是随机的 | 抛硬币的结果序列 |
| AR系统（自回归） | 当前值依赖过去几个值 | 气温的日间波动 |
| MA系统（滑动平均） | 当前值依赖过去的噪声项 | 某些传感器数据 |
| ARMA系统 | 既有AR又有MA成分 | 大多数平稳经济数据 |
| ARIMA系统 | 非平稳数据，需要差分后才是ARMA | 股价、GDP、气温 |
| SARIMA系统 | ARIMA + 季节性 | 有年/月周期的数据 |
| 状态空间系统 | 数据由隐藏的”状态”驱动 | 经济周期、库存变化 |
怎么用ACF和PACF来识别？
这是过程识别的核心技巧。对于平稳的ARMA类数据，ACF和PACF的”截尾”和”拖尾”特征是识别模型阶数的关键：
# 做一个对比示例：AR(1)、MA(1)和ARMA(1,1)的ACF/PACF特征
np.random.seed(123)
n = 500
# AR(1): y_t = 0.7*y_{t-1} + ε_t
ar1 = np.zeros(n)
ar1[0] = np.random.normal(0, 1)
for i in range(1, n):
ar1[i] = 0.7 * ar1[i-1] + np.random.normal(0, 1)
# MA(1): y_t = ε_t + 0.5*ε_{t-1}
epsilon = np.random.normal(0, 1, n)
ma1 = np.zeros(n)
for i in range(1, n):
ma1[i] = epsilon[i] + 0.5 * epsilon[i-1]
# ARMA(1,1): y_t = 0.5*y_{t-1} + ε_t + 0.3*ε_{t-1}
arma11 = np.zeros(n)
arma11[0] = np.random.normal(0, 1)
for i in range(1, n):
arma11[i] = 0.5 * arma11[i-1] + epsilon[i] + 0.3 * epsilon[i-1]
fig, axes = plt.subplots(3, 2, figsize=(14, 10))
# AR(1)
plot_acf(ar1, lags=40, ax=axes[0,0], title='AR(1) - ACF', alpha=0.05)
plot_pacf(ar1, lags=40, ax=axes[0,1], method='ywm', title='AR(1) - PACF', alpha=0.05)
# MA(1)
plot_acf(ma1, lags=40, ax=axes[1,0], title='MA(1) - ACF', alpha=0.05)
plot_pacf(ma1, lags=40, ax=axes[1,1], method='ywm', title='MA(1) - PACF', alpha=0.05)
# ARMA(1,1)
plot_acf(arma11, lags=40, ax=axes[2,0], title='ARMA(1,1) - ACF', alpha=0.05)
plot_pacf(arma11, lags=40, ax=axes[2,1], method='ywm', title='ARMA(1,1) - PACF', alpha=0.05)
for ax in axes.flatten():
ax.grid(True, alpha=0.2)
plt.tight_layout()
plt.savefig('acf_pacf_patterns.png', dpi=150)
plt.show()
看完这些图，你应该能总结出规律：
AR(p) 模型的特征：
- ACF：拖尾（逐渐衰减，不突然截断）
- PACF：在 lag=p 处截尾（p阶之后突然变得不显著）
MA(q) 模型的特征：
- ACF：在 lag=q 处截尾
- PACF：拖尾
ARMA(p,q) 模型的特征：
- ACF和PACF都拖尾
记忆口诀：AR看PACF，MA看ACF，ARMA两边看。AR的PACF在p阶”刹车”，MA的ACF在q阶”刹车”。
但现实中的数据不会这么听话
上面那些教科书式的例子，是理想情况下的纯ARMA数据。实际数据往往更复杂：
- 有趋势 → 先差分
- 有季节性 → 考虑SARIMA或季节性分解
- 有外生变量 → 考虑ARIMAX或状态空间模型
- 有结构性突变 → 需要在模型中加入虚拟变量或分段处理
常见错误：这些坑我踩过，希望你别踩
错误一：看到数据就急着跑模型
这是最常见的错误。很多人拿到数据，打开R或Python，直接 auto.arima() 就跑了。结果模型效果不好，还找不到原因。
正确做法：先看数据图 → 做分解 → 看ACF/PACF → 检验平稳性 → 再选模型。这个过程不能省，每一步都是为了给模型选择提供依据。
错误二：不理解”白噪声”的意义
白噪声就是完全没有规律的纯随机数据。如果你的数据残差是白噪声，说明模型已经提取了所有可预测的信息，这是好模型的目标。
但很多人分不清”原始数据是白噪声”和”残差是白噪声”。前者意味着你这数据根本没得建模——它本来就是随机的；后者才是建模成功的标志。
用代码验证残差是否为白噪声：
from statsmodels.stats.diagnostic import acorr_ljungbox
# 假设你已经拟合了模型，得到残差
# residuals = model.resid
# Ljung-Box检验：检验残差是否为白噪声
lb_test = acorr_ljungbox(residuals, lags=[10, 20], return_df=True)
print(lb_test)
# 如果p值都大于0.05，说明残差是白噪声，模型拟合良好
错误三：混淆相关性和因果性
ACF高，不代表一个值”导致”了另一个值。它只说明存在统计依赖关系。比如气温的ACF在lag=1很高，不代表”昨天的气温导致了今天的气温”，而是说明两者有共同的物理机制驱动。
建模时要结合领域知识来判断因果关系，不能光靠统计数据。
错误四：阶数选择不当
很多人用 auto.arima() 自动选择阶数，觉得这样最科学。但自动选择也有问题：
- 它在小样本上容易选过头（过拟合）
- 它可能忽略季节性
- 它不保证你得到的模型在业务上有意义
更好的做法：先用ACF/PACF图确定候选阶数范围，再用AIC/BIC等准则辅助选择，最后结合业务理解做判断。
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')
# 手动尝试几个候选模型，比较AIC
candidate_models = []
for p in range(0, 4):
for d in range(0, 3):
for q in range(0, 4):
try:
model = ARIMA(df_temp['temperature'], order=(p, d, q))
fitted = model.fit()
candidate_models.append({
'order': (p, d, q),
'AIC': fitted.aic,
'BIC': fitted.bic,
'loglik': fitted.loglike
})
except:
continue
candidates_df = pd.DataFrame(candidate_models)
# 按AIC排序，取前5个
top5 = candidates_df.nsmallest(5, 'AIC')
print(top5)
错误五：忽略结构性突变
时间序列中偶尔会有”突变点”——某个时间点数据的行为模式突然变了。比如2020年3月疫情爆发，全球航空销量断崖式下跌；2022年俄乌冲突导致欧洲天然气价格飙升。
如果不识别这些突变点，直接建模，模型会被”带偏”，预测效果很差。
from ruptures import PELT # 突变点检测库
# 找到突变点
signal = df_temp['temperature'].values.reshape(-1, 1)
# 注意：这里只是示例，气温数据可能不太需要找突变点
# 实际应用中，比如电价、销量等数据更常见
# 用BIC惩罚的PELT算法
algo = PELT(model="rbf").fit(signal)
result = algo.predict(pen=10)
print(f'突变点位置: {result}')
完整流程总结：从噪声到模型的实战路径
我给自己和团队成员总结了一个”五步法”，每次拿到新数据都按这个流程走，基本没踩过大的坑：
第一步：可视化
- 画原始时序图
- 画分解图（趋势、季节、残差）
- 目的：建立对数据的”手感”
第二步：平稳性检验
- ADF检验或KPSS检验
- 如果不平稳，做差分
- 目的：确认数据是否满足模型前提
第三步：ACF/PACF分析
- 看ACF和PACF的截尾/拖尾特征
- 确定ARIMA的p和q候选范围
- 目的：为模型阶数选择提供依据
第四步：模型拟合与诊断
- 尝试多个候选模型
- 比较AIC/BIC
- 做残差诊断（Ljung-Box检验、正态性检验）
- 目的：找到最优模型并验证其可靠性
第五步：业务验证
- 把模型结果和领域知识对比
- 做样本外预测，看实际效果
- 目的：确保模型不只是统计上好看，业务上也合理
一个完整的实战案例
我们用一份真实的销售数据来演示整个流程。假设你是某电商平台的数据分析师，公司想预测未来30天的某品类销量。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings
warnings.filterwarnings('ignore')
# ========== 1. 生成模拟销售数据 ==========
np.random.seed(2024)
dates = pd.date_range('2022-01-01', '2024-06-30', freq='D')
n = len(dates)
# 趋势：整体缓慢增长（电商发展）
trend = 1000 + 0.5 * np.arange(n)
# 年度季节性：双峰（618和双11）
seasonal = (
300 * np.exp(-((np.arange(n) % 365) - 170) ** 2 / (2 * 30 ** 2)) + # 618
500 * np.exp(-((np.arange(n) % 365) - 335) ** 2 / (2 * 25 ** 2)) + # 双11
100 * np.sin(2 * np.pi * np.arange(n) / 30) # 月度周期
)
# 噪声
noise = np.random.normal(0, 80, n)
# 销量数据
sales = trend + seasonal + noise
sales = np.maximum(sales, 100) # 销量不能为负
df_sales = pd.DataFrame({'date': dates, 'sales': sales})
# ========== 2. 可视化原始数据 ==========
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
axes[0].plot(df_sales['date'], df_sales['sales'], color='darkblue', linewidth=0.8)
axes[0].set_title('E-commerce Daily Sales (2022-2024)', fontsize=13)
axes[0].set_ylabel('Sales', fontsize=11)
axes[0].grid(True, alpha=0.3)
# 放大最近一年的数据，看得更清楚
recent = df_sales[df_sales['date'] >= '2023-07-01']
axes[1].plot(recent['date'], recent['sales'], color='darkblue', linewidth=0.6)
axes[1].set_title('Sales (Last 12 Months)', fontsize=13)
axes[1].set_ylabel('Sales', fontsize=11)
axes[1].set_xlabel('Date', fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('sales_raw.png', dpi=150)
plt.show()
从这张图中你可以清楚看到：
- 整体呈上升趋势（趋势成分）
- 有明显的高峰（618和双11促销活动）
- 日常有随机波动（噪声）
# ========== 3. 时间序列分解 ==========
df_sales.set_index('date', inplace=True)
# 尝试季节性分解（周期30天，捕捉月度规律）
result = seasonal_decompose(df_sales['sales'], period=30, model='additive')
fig, axes = plt.subplots(4, 1, figsize=(14, 10))
axes[0].plot(result.observed)
axes[0].set_title('Observed (原始销售数据)', fontsize=12)
axes[1].plot(result.trend)
axes[1].set_title('Trend (趋势)', fontsize=12)
axes[2].plot(result.seasonal)
axes[2].set_title('Seasonal (季节/周期成分)', fontsize=12)
axes[3].plot(result.resid)
axes[3].set_title('Residual (残差)', fontsize=12)
for ax in axes:
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('sales_decompose.png', dpi=150)
plt.show()
分解后注意到：30天周期只能捕捉到月度规律，双峰（618和双11）这种半年级的季节性需要更长的周期。可以尝试 period=365 来分解年度季节性。
# ========== 4. ADF检验平稳性 ==========
def adf_test(series, name):
result = adfuller(series.dropna())
print(f'=== {name} ADF检验 ===')
print(f'ADF统计量: {result[0]:.4f}, p值: {result[1]:.6f}')
print(f'5%临界值: {result[4]["5%"]:.4f}')
if result[1] < 0.05:
print('结论: 平稳 ✓')
else:
print('结论: 不平稳，需要差分 ✗')
return result[1] < 0.05
# 检验原始数据
is_stationary = adf_test(df_sales['sales'], '原始数据')
# 检验一阶差分
df_sales['sales_diff1'] = df_sales['sales'].diff()
is_stationary_diff1 = adf_test(df_sales['sales_diff1'], '一阶差分后')
# 检验二阶差分
df_sales['sales_diff2'] = df_sales['sales_diff1'].diff()
is_stationary_diff2 = adf_test(df_sales['sales_diff2'], '二阶差分后')
输出结果大致是：
=== 原始数据 ADF检验 ===
ADF统计量: -1.8234, p值: 0.385621
5%临界值: -2.8618
结论: 不平稳，需要差分 ✗
=== 一阶差分后 ADF检验 ===
ADF统计量: -8.2341, p值: 0.000001
5%临界值: -2.8619
结论: 平稳 ✓
# ========== 5. ACF和PACF分析（对平稳后的数据） ==========
# 用一阶差分后的数据
diff_series = df_sales['sales_diff1'].dropna()
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plot_acf(diff_series, lags=60, ax=axes[0], alpha=0.05)
axes[0].set_title('ACF of 1st Difference (一阶差分的自相关)', fontsize=12)
plot_pacf(diff_series, lags=60, ax=axes[1], method='ywm', alpha=0.05)
axes[1].set_title('PACF of 1st Difference (一阶差分的偏自相关)', fontsize=12)
for ax in axes:
ax.grid(True, alpha=0.2)
plt.tight_layout()
plt.savefig('sales_acf_pacf.png', dpi=150)
plt.show()
从ACF/PACF图中，你可能观察到：
- ACF在 lag=1 处显著，然后逐渐衰减 → 提示MA(1)成分
- PACF在 lag=1 处显著，然后逐渐衰减 → 提示AR(1)成分
- 也可能在 lag=7, 14, 21 等有周期性峰值 → 提示周季节性
根据这些信息，候选模型可能是 ARIMA(1,1,1) 或 ARIMA(1,1,1)(1,1,1)7（带周季节性）。
# ========== 6. 模型拟合与比较 ==========
# 准备数据
train = df_sales['sales'][:'2024-03-01']
test = df_sales['sales']['2024-03-02':]
# 尝试几个候选模型
models = {
'ARIMA(1,1,1)': (1, 1, 1),
'ARIMA(2,1,1)': (2, 1, 1),
'ARIMA(1,1,2)': (1, 1, 2),
'ARIMA(1,1,1)(1,1,1)7': (1, 1, 1, (1, 1, 1, 7)),
}
results = []
for name, order in models.items():
try:
model = ARIMA(train, order=order)
fitted = model.fit()
# 在测试集上预测
forecast = fitted.forecast(steps=len(test))
# 计算RMSE
rmse = np.sqrt(np.mean((forecast - test.values) ** 2))
results.append({
'model': name,
'AIC': fitted.aic,
'BIC': fitted.bic,
'RMSE': rmse
})
print(f'{name}: AIC={fitted.aic:.2f}, RMSE={rmse:.2f}')
except Exception as e:
print(f'{name}: 拟合失败 - {e}')
results_df = pd.DataFrame(results)
print('\n模型比较结果:')
print(results_df.sort_values('RMSE'))
# ========== 7. 残差诊断 ==========
best_model_name = 'ARIMA(1,1,1)(1,1,1)7'
best_order = (1, 1, 1, (1, 1, 1, 7))
model_best = ARIMA(train, order=best_order)
fitted_best = model_best.fit()
# 残差
residuals = fitted_best.resid
# Ljung-Box检验
lb_result = acorr_ljungbox(residuals, lags=[10, 20], return_df=True)
print('Ljung-Box检验:')
print(lb_result)
# 残差图
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes[0,0].plot(residuals)
axes[0,0].set_title('Residuals (残差)', fontsize=11)
axes[0,0].axhline(0, color='red', linestyle='--')
plot_acf(residuals, lags=30, ax=axes[0,1], alpha=0.05)
axes[0,1].set_title('Residual ACF (残差自相关)', fontsize=11)
axes[1,0].hist(residuals, bins=30, density=True, alpha=0.6, color='steelblue')
axes[1,0].set_title('Residual Distribution (残差分布)', fontsize=11)
# QQ图
from scipy import stats
stats.probplot(residuals, dist="norm", plot=axes[1,1])
axes[1,1].set_title('Residual QQ-Plot (残差正态Q-Q图)', fontsize=11)
plt.tight_layout()
plt.savefig('residual_diagnostic.png', dpi=150)
plt.show()
如果Ljung-Box检验的p值都大于0.05，说明残差是白噪声，模型拟合良好。QQ图接近对角线，说明残差近似正态分布。
# ========== 8. 预测可视化 ==========
forecast = fitted_best.forecast(steps=len(test))
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(train.index, train.values, label='Training Data (训练数据)', color='steelblue')
ax.plot(test.index, test.values, label='Actual (实际值)', color='green', linewidth=1.5)
ax.plot(test.index, forecast, label='Forecast (预测值)', color='red', linewidth=1.5, alpha=0.8)
ax.fill_between(forecast.index,
forecast - 1.96 * 80,
forecast + 1.96 * 80,
color='red', alpha=0.1, label='95% Confidence Interval')
ax.set_title('Sales Forecast (销量预测)', fontsize=14)
ax.set_xlabel('Date', fontsize=12)
ax.set_ylabel('Sales', fontsize=12)
ax.legend(loc='upper left', fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('sales_forecast.png', dpi=150)
plt.show()
最后几句大实话
过程识别这一步，说难也难，说不难也不难。难在于你需要同时看懂图、理解统计检验、结合业务背景做判断；不难在于只要按流程走，每一步都有明确的操作和判断标准，不会迷路。
我见过太多人跳过过程识别直接建模，然后对着结果发愁。其实时间序列建模就像医生看病：你不能一上来就开刀，得先问诊（看数据图）、做检查（平稳性检验、ACF/PACF）、确诊（识别系统类型），最后再开药（选模型）。
记住一句话：好的模型不是跑出来的，是”看”出来的。 多花时间在前期的探索性分析上，后面的建模和预测才会顺理成章。
如果你现在手里有一份新的时间序列数据，不妨按照上面的流程走一遍。看完原始图、做完分解、看了ACF/PACF、做了ADF检验之后，你自然就知道该用什么样的模型了。这个过程不需要天赋，只需要一点耐心和正确的方法。
好了，今天就聊到这里。如果还有具体问题，随时来问我。