一、时间序列基础概念
时间序列（Time Series）是按时间顺序排列的一系列观测值。与普通回归数据不同，时间序列的核心特征是观测值之间存在时间依赖结构，即当前值往往与过去的值相关。理解这种依赖关系是进行预测的基础。
1.1 时间序列的三个基本成分
任何时间序列理论上都可以分解为以下三个基本成分的组合，识别这些成分是建立预测模型的第一步：
- 趋势（Trend）： 序列长期单调上升或下降的走向。例如气温上升趋势、用户数增长趋势。趋势可以是线性的，也可以是非线性的（如指数型、饱和型）。
- 季节性（Seasonality）： 以固定周期重复出现的波动模式。例如零售业的年末销售高峰、电力消耗的日内双峰形态。季节性周期通常为年、季、月、周、日。
- 残差（Residual）： 去除趋势和季节性后的随机波动部分，也称噪声。理想情况下残差应为白噪声（均值为零、方差恒定的随机序列）。
在经典时间序列理论中，这三个成分通过特定的数学模型组合起来，形成了加法模型和乘法模型两种基本框架。
1.2 加法模型与乘法模型
加法模型（Additive Model）假设各成分独立叠加，乘法模型（Multiplicative Model）假设成分之间存在比例关系。选择哪种模型取决于序列的波动幅度是否随趋势水平变化：
加法模型
Y[t] = Trend[t] + Seasonal[t] + Residual[t]
适用于季节性波动幅度不随趋势变化的场景。例如月平均气温序列，无论整体温度上升还是下降，季节性波动的绝对幅度基本不变。
Python 实现：
from statsmodels.tsa.seasonal import seasonal_decompose
result = seasonal_decompose(
series, model='additive', period=12
)
trend = result.trend
seasonal = result.seasonal
residual = result.resid
# 可视化
result.plot()
乘法模型
Y[t] = Trend[t] x Seasonal[t] x Residual[t]
适用于季节性波幅随趋势水平同步变化（比例恒定）的场景。例如电商销售额，当销售额从100万增长到1000万时，双十一的促销峰值也会同比放大。
Python 实现：
from statsmodels.tsa.seasonal import seasonal_decompose
result = seasonal_decompose(
series, model='multiplicative',
period=12
)
trend = result.trend
seasonal = result.seasonal
residual = result.resid
如何判断使用加法还是乘法？
- 画出序列图，观察季节波动的振幅是否随时间（趋势水平）变化
- 振幅基本不变 → 加法模型；振幅与趋势水平成正比 → 乘法模型
- 不确定时，两种都试，比较残差的白噪声程度
- 乘法模型取对数后可转化为加法模型：log(Y) = log(Trend) + log(Seasonal) + log(Residual)
二、时间序列分解方法
时间序列分解是将原始序列拆分为趋势、季节性和残差的过程。分解是理解数据内在结构的关键步骤，也是后续平稳化处理的前置工作。除了经典方法外，还有更鲁棒的STL分解。
2.1 经典分解（seasonal_decompose）
statsmodels.tsa.seasonal.seasonal_decompose 是最常用的分解函数，底层使用滑动平均法提取趋势，再通过周期平均提取季节成分。它的计算步骤可以归纳为：
- 使用中心化滑动平均（CMA）估计趋势分量 Trend[t]
- 从原始序列中减去趋势（加法）或除以趋势（乘法），得到 detrended 序列
- 按周期对 detrended 序列取平均，得到季节性分量 Seasonal[t]
- 从 detrended 中减去（或除以）季节性，得到残差 Residual[t]
经典分解的优势在于简单直观，但存在两个明显缺陷：一是要求周期严格固定（如 period=12 表示每年12个月），二是对异常值非常敏感，一个离群点就会污染整个滑动窗口。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
# 生成模拟数据：趋势 + 季节 + 噪声
np.random.seed(42)
t = np.arange(120)
trend = 10 + 0.05 * t
seasonal = 2 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 0.5, 120)
series = trend + seasonal + noise
# 转换为 pandas Series
ts = pd.Series(series, index=pd.date_range(
'2020-01-01', periods=120, freq='M'
))
# 加法分解
decomp = seasonal_decompose(ts, model='additive', period=12)
print("趋势分量前5个值:", decomp.trend.head().values)
print("季节分量前12个值:", decomp.seasonal[:12].values)
print("残差标准差:", np.nanstd(decomp.resid))
2.2 STL 分解（Seasonal-Trend decomposition using LOESS）
STL 分解由 Cleveland 等人于1990年提出，使用局部加权回归（LOESS）进行平滑和分解，是业界公认最稳健的分解方法。与经典分解相比，STL 具有以下优势：
- 鲁棒性强： 对异常值不敏感，内置鲁棒性权重迭代
- 季节成分可随时间变化： 允许季节性模式的渐变（如零售业中节假日促销效应的逐年演变）
- 灵活控制平滑度： 通过 seasonal、trend 等参数独立控制各成分的平滑程度
- 缺失值处理： 对数据中的小缺口有一定容忍度
from statsmodels.tsa.seasonal import STL
# STL分解 -- 推荐用于大多数实际场景
stl = STL(ts, period=12,
seasonal=13, # 季节窗口（奇数）
trend=None, # None 时自动选择
seasonal_deg=1, # 季节LOESS多项式阶数
trend_deg=1, # 趋势LOESS多项式阶数
robust=True) # 启用鲁棒性
result = stl.fit()
trend_stl = result.trend
seasonal_stl = result.seasonal
residual_stl = result.resid
# 查看残差统计
print("STL残差均值:", residual_stl.mean())
print("STL残差标准差:", residual_stl.std())
# 提取后的数据可直接用于后续建模
detrended = ts - result.trend
deseasonalized = ts - result.seasonal
STL 参数调优建议
seasonal（季节平滑窗口）：必须为奇数。越大则季节成分越平滑，默认值为 period+1 取下一个奇数。trend（趋势平滑窗口）：约为 period 的1.5倍。通过 STL(..., seasonal=13, trend=19) 分别控制。low_pass：低频滤波窗口，通常取 season 的最大奇数。设置 robust=True 可以显著提高对异常脉冲击波的抗干扰能力。
2.3 季节调整与成分提取
分解完成后，可以提取后续建模所需的数据版本。在实际项目中，分解提纯的常见做法包括：
- 季节调整（Seasonal Adjustment）： adjusted = ts - seasonal，得到去除季节效应后的"经季节调整"序列，常用于官方经济数据发布（如CPI、GDP）。
- 去趋势（Detrending）： detrended = ts - trend，得到平稳化的序列用于进一步建模。
- 提取残差做异常检测： 残差中超出3倍标准差的点往往是异常事件。
- 成分重组： 各成分可独立预测后再相加，即"分解-预测-重组"范式。
# 成分提取与重组：分别预测各成分再相加的范例
# 1. 提取各成分
trend_part = result.trend.dropna()
seasonal_part = result.seasonal.dropna()
resid_part = result.resid.dropna()
# 2. 对趋势做简单的线性外推预测（实践中可用ARIMA等）
from sklearn.linear_model import LinearRegression
X_trend = np.arange(len(trend_part)).reshape(-1, 1)
y_trend = trend_part.values
lr = LinearRegression().fit(X_trend, y_trend)
# 预测未来12期趋势
future_idx = np.arange(len(trend_part), len(trend_part)+12).reshape(-1, 1)
trend_forecast = lr.predict(future_idx)
# 3. 季节成分直接使用已知周期（或做季节预测）
seasonal_forecast = seasonal_part[:12].values # 简单复用
# 4. 组合预测
forecast = trend_forecast + seasonal_forecast
print("未来12期预测值:", forecast)
三、平稳性检验与差分
平稳性（Stationarity）是经典时间序列建模的前提条件。一个平稳时间序列的统计性质（均值、方差、自协方差）不随时间改变。大多数时间序列模型（如ARIMA）要求输入序列是平稳的，否则会产生"伪回归"问题。
3.1 平稳性的直观理解与检验必要性
如果一个序列的均值随时间上升（如股价长期上涨），或者方差不恒定（如波动率聚集），它就是非平稳的。非平稳序列的统计分析会得出误导性结论：两列不相关的随机游走可能算出显著的相关性。因此，建模前必须进行平稳性检验。
3.2 ADF 检验（Augmented Dickey-Fuller Test）
ADF 检验是检验单位根（Unit Root）是否存在的最常用方法。原假设 H0 为"序列存在单位根（非平稳）"，备择假设 H1 为"序列平稳"。当 p-value 小于显著性水平（通常 0.05）时，拒绝 H0，认为序列平稳。
from statsmodels.tsa.stattools import adfuller
def adf_test(series, title=''):
"""
对序列执行ADF检验并打印详细结果
"""
result = adfuller(series, autolag='AIC')
print(f'=== ADF检验: {title} ===')
print(f'ADF统计量: {result[0]:.6f}')
print(f'p-value: {result[1]:.6f}')
print(f'临界值:')
for key, val in result[4].items():
print(f' {key}: {val:.6f}')
print(f'结论: {"序列平稳 ✓" if result[1] < 0.05 else "序列非平稳 ✗"}')
return result[1] < 0.05
# 对原始序列检验
adf_test(ts, '原始数据')
# 对一阶差分检验
ts_diff = ts.diff().dropna()
adf_test(ts_diff, '一阶差分')
3.3 KPSS 检验
KPSS 检验是 ADF 的互补方法。它的原假设与 ADF 相反：H0 为"序列平稳"，H1 为"存在单位根"。同时使用 ADF 和 KPSS 可以减少误判风险：
- ADF 显著 + KPSS 不显著 → 序列平稳
- ADF 不显著 + KPSS 显著 → 序列非平稳
- 两者都显著 → 序列可能有趋势平稳特征
- 两者都不显著 → 数据量不足，需要更多观察
from statsmodels.tsa.stattools import kpss
def kpss_test(series, title=''):
"""
对序列执行KPSS检验
"""
result = kpss(series, regression='c', nlags='auto')
print(f'=== KPSS检验: {title} ===')
print(f'KPSS统计量: {result[0]:.6f}')
print(f'p-value: {result[1]:.6f}')
print(f'结论: {"序列非平稳" if result[1] < 0.05 else "序列平稳 ✓"}')
# 联合检验
adf_result = adfuller(ts, autolag='AIC')
kpss_result = kpss(ts, regression='c', nlags='auto')
if adf_result[1] < 0.05 and kpss_result[1] > 0.05:
print("判定结论：数据是平稳的")
elif adf_result[1] > 0.05 and kpss_result[1] < 0.05:
print("判定结论：数据是非平稳的，需要进行差分")
else:
print("判定结论：结果不一致，建议进一步分析")
3.4 差分与季节差分
差分是最常用的平稳化手段。一阶差分 diff(1) 可以消除线性趋势，二阶差分消除二次趋势。对于周期为 m 的季节性数据，使用季节差分（Seasonal Differencing）：
- 普通差分： diff(1) = Y[t] - Y[t-1]，消除趋势
- 季节差分： diff(m) = Y[t] - Y[t-m]，消除季节性（m 为周期长度）
- 双重差分： 先季节差分再一阶差分，即 diff(m).diff(1)
import pandas as pd
import numpy as np
# 模拟带有趋势和季节性的数据
np.random.seed(123)
t = np.arange(200)
trend = 0.1 * t
seasonal = 3 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 0.8, 200)
data = trend + seasonal + noise
ts = pd.Series(data, index=pd.date_range('2020-01-01', periods=200, freq='M'))
# 普通一阶差分 -- 消除趋势
ts_diff1 = ts.diff(1).dropna()
# 季节差分（周期12） -- 消除季节性
ts_diff12 = ts.diff(12).dropna()
# 季节差分后再次一阶差分 -- 完全平稳化
ts_diff12_1 = ts_diff12.diff(1).dropna()
# 检验效果
print("原始序列ADF p-value:", adfuller(ts, autolag='AIC')[1])
print("一阶差分ADF p-value:", adfuller(ts_diff1, autolag='AIC')[1])
print("季节差分ADF p-value:", adfuller(ts_diff12.dropna(), autolag='AIC')[1])
print("双重差分ADF p-value:", adfuller(ts_diff12_1, autolag='AIC')[1])
差分的注意事项
- 不要过度差分：每差分一次就会丢失一个观测值，且增加模型复杂度
- 大多数经济金融序列一阶差分后就已平稳，很少需要 d > 2
- 差分阶数 d 确定了 ARIMA(p,d,q) 模型的核心参数之一
- 季节差分应用于有明确周期性的数据（周期 m 已知）
四、自相关与偏自相关分析
自相关函数（ACF）和偏自相关函数（PACF）是识别时间序列模型结构（确定 AR 和 MA 的阶数）的核心可视化工具，也是连接平稳性分析与模型选择的桥梁。
4.1 ACF（自相关函数）
ACF 衡量序列与其自身滞后版本之间的线性相关性。滞后 k 的自相关 r[k] 是 Y[t] 和 Y[t-k] 之间的相关系数。如果 ACF 在若干滞后后截尾（突然降为不显著），或者缓慢衰减，这提供了判断序列结构的重要线索。
4.2 PACF（偏自相关函数）
PACF 衡量在剔除中间滞后项的影响后，Y[t] 和 Y[t-k] 之间的净相关性。例如，PACF 在滞后3处显著，说明 Y[t] 和 Y[t-3] 之间存在不能被 Y[t-1]、Y[t-2] 解释的额外相关关系，这在识别 AR 模型的阶数时尤为关键。
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
# 确定画布大小
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
# 原始序列的ACF和PACF
plot_acf(ts, lags=30, ax=axes[0, 0], title='原始序列 ACF')
plot_pacf(ts, lags=30, ax=axes[0, 1], title='原始序列 PACF')
# 一阶差分后序列的ACF和PACF
ts_diff = ts.diff().dropna()
plot_acf(ts_diff, lags=30, ax=axes[1, 0], title='一阶差分后 ACF')
plot_pacf(ts_diff, lags=30, ax=axes[1, 1], title='一阶差分后 PACF')
plt.tight_layout()
plt.show()
4.3 利用 ACF/PACF 识别 ARMA 阶数
ACF 和 PACF 的形态可以指导 ARIMA(p,d,q) 模型的选择：
| 模型 |
ACF 形态 |
PACF 形态 |
| AR(p) -- 自回归 |
拖尾：指数衰减或震荡衰减 |
p步后截尾（突然降为0） |
| MA(q) -- 移动平均 |
q步后截尾 |
拖尾：指数衰减 |
| ARMA(p,q) -- 混合 |
拖尾（p步后开始衰减） |
拖尾（q步后开始衰减） |
| 随机游走 |
缓慢衰减，不截尾 |
滞后1处显著，之后截尾 |
实际识别步骤
- 观察 PACF：看哪个滞后之后突然变得不显著（蓝色置信带内）→ 确定为 AR 阶数 p
- 观察 ACF：看哪个滞后之后突然变得不显著 → 确定为 MA 阶数 q
- 如果两者都拖尾，尝试 ARMA(1,1) 并迭代调整
- 始终使用 AIC/BIC 辅助判断，避免仅靠肉眼识别造成的过拟合
五、ARIMA 模型
ARIMA（自回归积分滑动平均模型）是时间序列分析中最经典、应用最广泛的预测模型。它将三个核心组件统一在一个框架中：AR（自回归）捕捉序列的"记忆"效应，I（差分）实现平稳化，MA（移动平均）捕捉冲击的衰减效应。
5.1 ARIMA 模型的数学表达与直觉理解
ARIMA(p,d,q): (1 - φ₁B - ... - φₚBᵖ)(1-B)ᵈY[t] = (1 + θ₁B + ... + θₙBᵠ)ε[t]
其中 B 是滞后算子（B Y[t] = Y[t-1]）。用直觉来理解这三个部分：
- AR (p) -- "惯性"部分： 当前值是过去 p 期值的线性组合。例如 AR(1) 意味着 Y[t] = φ₁Y[t-1] + ε[t]，今天的值部分地由昨天决定。φ₁ 接近1表示强惯性（持久性），φ₁ 接近0表示弱记忆。
- I (d) -- "变化"部分： 对序列进行 d 次差分以消除趋势。d=1 意味着用"变化量"（增量）而非"水平值"建模。
- MA (q) -- "冲击"部分： 当前值是过去 q 期预测误差的加权和。MA(1) 表示 Y[t] = ε[t] + θ₁ε[t-1]，即今天的值受到昨天"意外冲击"的影响。
5.2 ARIMA 建模完整流程
import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
# === 步骤1: 准备数据 ===
np.random.seed(42)
dates = pd.date_range('2020-01-01', periods=200, freq='M')
data = (10 + 0.05 * np.arange(200) +
3 * np.sin(2 * np.pi * np.arange(200) / 12) +
np.random.normal(0, 1, 200))
ts = pd.Series(data, index=dates)
# === 步骤2: 平稳性检验 ===
from statsmodels.tsa.stattools import adfuller
adf_pval = adfuller(ts, autolag='AIC')[1]
print(f"ADF p-value = {adf_pval:.4f}")
d = 0 if adf_pval < 0.05 else 1
print(f"确定 d = {d}")
# === 步骤3: 差分并查看ACF/PACF（略，见前文）===
# === 步骤4: 拟合ARIMA模型 ===
model = ARIMA(ts, order=(2, d, 2))
fitted = model.fit()
# === 步骤5: 模型诊断 ===
print(fitted.summary())
# 残差检验：残差应为白噪声
residuals = fitted.resid
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
residuals.plot(ax=axes[0], title='残差序列')
axes[0].axhline(y=0, linestyle='--', color='gray')
plot_acf(residuals, lags=20, ax=axes[1], title='残差ACF')
residuals.hist(bins=30, ax=axes[2], alpha=0.7)
axes[2].set_title('残差直方图')
plt.tight_layout()
plt.show()
# === 步骤6: 预测 ===
forecast = fitted.forecast(steps=12)
print("未来12期预测值:")
print(forecast)
5.3 auto_arima 自动定阶
手动识别 ACF/PACF 并结合 AIC 选取参数是一个费时费力的过程。pmdarima 库提供的 auto_arima 函数可以自动搜索最优的 (p,d,q) 组合，极大提升建模效率。它基于 AIC 或 BIC 在指定参数空间中进行网格搜索：
# 需要安装: pip install pmdarima
from pmdarima import auto_arima
# 自动搜索最优ARIMA参数
auto_model = auto_arima(
ts,
start_p=0, max_p=5, # AR阶数搜索范围
start_q=0, max_q=5, # MA阶数搜索范围
d=None, # 自动确定差分阶数
seasonal=False, # 是否考虑季节性（SARIMA见下一章）
trace=True, # 打印搜索过程
error_action='ignore',
suppress_warnings=True,
stepwise=True, # 使用逐步搜索（更快）
information_criterion='aic', # 使用AIC选择最佳模型
max_order=10, # p+q 的最大值限制
)
print("\n最优模型参数:")
print(f"ARIMA{auto_model.order}")
print(f"AIC: {auto_model.aic():.2f}")
# 使用最佳模型预测
preds = auto_model.predict(n_periods=12)
print("\n自动ARIMA预测结果:", preds.values)
auto_arima 参数详解
stepwise=True 使用 Hyndman-Khandakar 算法（逐步搜索），比穷举搜索快数倍。d=None 让 auto_arima 自动执行 KPSS 检验来确定差分阶数。seasonal=False 表示非季节性模型。information_criterion='aic' 使用 AIC 作为选择标准（另一个选项是 'bic'，更倾向于简约模型）。max_order=10 限制 p+q 的总和，防止模型过于复杂。
六、季节性 ARIMA（SARIMA）
当数据存在显著的季节性模式时，仅靠 ARIMA(p,d,q) 是不够的。SARIMA（Seasonal ARIMA）在 ARIMA 的基础上增加了季节性分量，用大写字母 P,D,Q,m 表示：P 为季节自回归阶数，D 为季节差分阶数，Q 为季节移动平均阶数，m 为季节周期长度。
6.1 SARIMA 模型结构
SARIMA(p,d,q)(P,D,Q)[m]
其中 (p,d,q) 是非季节部分，(P,D,Q)[m] 是季节部分，m 是每个季节周期的观测数。例如 SARIMA(1,1,1)(1,1,1)[12] 表示：
- 非季节部分：AR阶数1，差分阶数1，MA阶数1
- 季节部分（月数据，年周期=12）：季节AR阶数1，季节差分阶数1，季节MA阶数1
- 等价于对序列做 diff(1).diff(12) 后，再拟合 ARMA(1,1)x(1,1)[12] 模型
6.2 SARIMA 建模示例
from statsmodels.tsa.statespace.sarimax import SARIMAX
# 构建SARIMA模型: ARIMA(1,1,1)(1,1,1)[12]
model = SARIMAX(
ts,
order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False
)
fitted = model.fit(disp=False)
print(fitted.summary())
# 模型诊断
fitted.plot_diagnostics(figsize=(15, 10))
plt.show()
# 预测未来24期
forecast = fitted.get_forecast(steps=24)
pred_mean = forecast.predicted_mean
pred_ci = forecast.conf_int()
print("未来24期预测均值:")
print(pred_mean)
# 获取置信区间
lower = pred_ci.iloc[:, 0]
upper = pred_ci.iloc[:, 1]
6.3 自动选择季节阶数
auto_arima 同样支持 SARIMA 的自动定阶，只需设置 seasonal=True 并指定 m：
from pmdarima import auto_arima
# SARIMA自动定阶
sarima_auto = auto_arima(
ts,
start_p=0, max_p=3,
start_q=0, max_q=3,
d=None,
seasonal=True,
m=12, # 季节周期长度
start_P=0, max_P=2,
start_Q=0, max_Q=2,
D=None, # 自动确定季节差分阶数
trace=True,
stepwise=True,
error_action='ignore',
suppress_warnings=True,
max_order=8,
)
print(f"\n最优SARIMA模型: {sarima_auto.order} x {sarima_auto.seasonal_order}")
print(f"AIC: {sarima_auto.aic():.2f}")
如何确定季节周期 m？
- 月度数据：m=12（一年内12个月的周期）
- 季度数据：m=4（一年4个季度）
- 周数据（年度周期）：m=52（一年52周）
- 日数据（年度周期）：m=365（一年365天），或 m=7（周周期）
- 小时数据（日周期）：m=24（一天24小时）
- 当 m 很大时（如 m=365），计算量剧增，此时改用 Prophet 可能更合适
6.4 SARIMAX -- 加入外生变量
SARIMAX 是 SARIMA 的扩展，允许加入外生解释变量（Exogenous Variables）。例如在预测销售额时，可以将促销活动标志、节假日、天气数据等作为外生变量加入模型，提高预测准确性：
import pandas as pd
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX
# 模拟外生变量（如：促销活动标志）
np.random.seed(42)
exog = pd.DataFrame({
'promotion': np.random.binomial(1, 0.15, 200), # 15%的时间有促销
'holiday': np.random.binomial(1, 0.05, 200) # 5%的时间是节假日
}, index=ts.index)
# SARIMAX模型：加入外生变量
sarimax_model = SARIMAX(
ts,
exog=exog,
order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False
)
sarimax_fitted = sarimax_model.fit(disp=False)
print(sarimax_fitted.summary())
# 预测时需要同时提供未来的外生变量
future_exog = pd.DataFrame({
'promotion': np.random.binomial(1, 0.15, 12),
'holiday': np.random.binomial(1, 0.05, 12)
})
forecast_sarimax = sarimax_fitted.forecast(
steps=12, exog=future_exog
)
print("带外生变量的预测结果:", forecast_sarimax.values)
七、波动率分析与异常检测
在许多实际问题中，不仅序列的水平值重要，序列的波动性（Volatility）同样携带关键信息。金融市场的风险度量、工业设备的状态监测、服务器负载的异常告警，都依赖波动率分析。
7.1 滚动标准差（Rolling Standard Deviation）
最简单的波动率估计是在固定窗口内计算标准差。窗口越小对短期波动的追踪越快，但也越容易被噪声干扰；窗口越大越平滑，但对波动的响应延迟越大。
import pandas as pd
import numpy as np
# 模拟有波动率聚集效应的序列
np.random.seed(42)
returns = np.random.normal(0, 0.01, 500)
volatility_clusters = np.where(
np.abs(returns) > 0.02,
returns * 1.5, # 大波动后放大波动
returns
)
price = 100 + np.cumsum(volatility_clusters) # 模拟价格
ts_price = pd.Series(price, index=pd.date_range(
'2022-01-01', periods=500, freq='D'
))
# 计算滚动标准差（30天窗口）
rolling_std_30 = ts_price.pct_change().rolling(30).std()
rolling_std_60 = ts_price.pct_change().rolling(60).std()
# 将波动率超过阈值（90%分位数）的点标记为异常高波动
threshold = rolling_std_30.quantile(0.90)
high_volatility = rolling_std_30 > threshold
print(f"高波动阈值: {threshold:.4f}")
print(f"高波动天数: {high_volatility.sum()} / {len(rolling_std_30)}")
# 绘制日收益率与滚动标准差
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(14, 8))
ts_price.pct_change().plot(ax=axes[0], title='日收益率', alpha=0.6)
rolling_std_30.plot(ax=axes[1], label='30日滚动标准差', color='red')
rolling_std_60.plot(ax=axes[1], label='60日滚动标准差', color='blue')
axes[1].axhline(y=threshold, linestyle='--', color='gray', label='90%阈值')
axes[1].legend()
plt.tight_layout()
plt.show()
7.2 波动率聚类（Volatility Clustering）
波动率聚类是时间序列中广泛存在的现象：大幅波动往往跟随着大幅波动，小幅波动也倾向于聚集。这是因为外部冲击的传导、市场参与者的情绪传染、风险管理的连锁反应等机制会导致波动具有自相关性。ARCH/GARCH 模型族可用于正式建模波动率聚类，但实践中滚动标准差+阈值法已经能在很多场景下提供可靠的异常检测信号。
7.3 CUSUM 变化点检测
CUSUM（Cumulative Sum，累计和控制图）是一种经典的变化点检测算法。它通过累积偏差的累加来放大微小的持续性偏移，非常适合检测序列均值的结构性变化。其核心思想是：如果序列均值发生偏移，累积和会快速偏离零。
def cusum_detection(series, threshold=5, drift=0.5):
"""
CUSUM变化点检测算法
参数:
- series: 输入序列
- threshold: 检测阈值（决定报警灵敏度）
- drift: 允许的微小漂移量（过滤掉不必要的报警）
返回:
- change_points: 检测到的变化点索引列表
"""
mean_val = np.mean(series)
deviations = series - mean_val
cum_sum_upper = np.zeros(len(series))
cum_sum_lower = np.zeros(len(series))
change_points = []
for i in range(1, len(series)):
# 上侧累积和：检测均值上升
cum_sum_upper[i] = max(
0, cum_sum_upper[i-1] + deviations.iloc[i] - drift
)
# 下侧累积和：检测均值下降
cum_sum_lower[i] = max(
0, cum_sum_lower[i-1] - deviations.iloc[i] - drift
)
# 超过阈值则标记为变化点并重置
if cum_sum_upper[i] > threshold:
change_points.append(i)
cum_sum_upper[i] = 0
if cum_sum_lower[i] > threshold:
change_points.append(i)
cum_sum_lower[i] = 0
return change_points, cum_sum_upper, cum_sum_lower
# 模拟带有均值偏移的序列
np.random.seed(42)
normal = np.random.normal(0, 1, 200)
shifted = np.random.normal(1.5, 1, 100) # 均值从0跳到1.5
combined = np.concatenate([normal, shifted])
ts_cusum = pd.Series(combined, index=pd.date_range(
'2022-01-01', periods=300, freq='D'
))
# 应用CUSUM检测
cps, upper, lower = cusum_detection(
ts_cusum, threshold=5, drift=0.3
)
print(f"检测到变化点位置: {cps}")
# 预期在第200个点左右检测到均值偏移
CUSUM 参数调优
threshold（阈值 h）：决定了检测的灵敏度。较小的值能更快发现变化但容易误报；较大的值减少误报但可能滞后检测。工程实践中，建议以交叉验证或历史数据回测来确定。drift（漂移参数 k）：允许序列在不触发报警的情况下自然漂移的幅度，通常设为可接受偏移量的一半。
八、Prophet 时间序列预测
Prophet 是 Facebook（现 Meta）在2017年开源的时间序列预测工具，由 Sean Taylor 和 Ben Letham 主导开发。它基于可加模型（Additive Model），将时间序列分解为趋势、季节性、节假日效应三个分量，特别适用于具有强季节效应和多周期的日常业务数据。
8.1 Prophet 的核心思想
y(t) = g(t) + s(t) + h(t) + ε[t]
其中 g(t) 是趋势函数（支持线性增长和逻辑斯蒂增长两种），s(t) 是季节性分量（用傅里叶级数逼近），h(t) 是节假日效应（用户自定义的哑变量），ε[t] 是误差项。Prophet 的几个关键创新点在于：
- 自动变化点检测（Changepoint）： 趋势可以在预设的时间点发生斜率变化，自动适应结构性突变
- 傅里叶季节逼近： 不使用固定周期平均，而是用傅里叶级数拟合季节性模式，支持多周期叠加
- 区间预测（Uncertainty Intervals）： 通过马尔可夫链蒙特卡洛（MCMC）或朴素模拟给出预测区间
- 缺失值和异常值鲁棒： 对数据质量和缺失值不太敏感
8.2 Prophet 完整建模流程
# 需要安装: pip install prophet
import pandas as pd
from prophet import Prophet
import matplotlib.pyplot as plt
# === 数据准备：Prophet要求两列 ===
# ds: 日期列（日期格式）
# y: 数值列（目标变量）
np.random.seed(42)
dates = pd.date_range('2020-01-01', periods=365*3, freq='D')
trend = 10 + 0.01 * np.arange(len(dates))
seasonal_yearly = 2 * np.sin(2 * np.pi * np.arange(len(dates)) / 365.25)
seasonal_weekly = 0.5 * np.sin(2 * np.pi * np.arange(len(dates)) / 7)
noise = np.random.normal(0, 0.5, len(dates))
data = trend + seasonal_yearly + seasonal_weekly + noise
df = pd.DataFrame({
'ds': dates,
'y': data
})
# === 创建并拟合Prophet模型 ===
model = Prophet(
growth='linear', # 增长模式：linear / logistic
seasonality_mode='additive', # 季节模式：additive / multiplicative
yearly_seasonality=True, # 年季节性（自动开启）
weekly_seasonality=True, # 周季节性（自动开启）
daily_seasonality=False, # 日季节性（数据为日频时开启）
changepoint_prior_scale=0.05, # 变化点灵敏度（越大越容易检测变化）
seasonality_prior_scale=10.0, # 季节性强度（越大季节效应越强）
interval_width=0.80, # 80%预测区间
uncertainty_samples=1000, # 不确定性模拟次数
)
model.add_country_holidays('CN') # 加入中国节假日效应
model.fit(df)
# === 构建未来日期 ===
future = model.make_future_dataframe(periods=365) # 预测未来1年
# === 预测 ===
forecast = model.predict(future)
# 查看关键预测列
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail())
# === 可视化 ===
fig1 = model.plot(forecast)
fig2 = model.plot_components(forecast) # 趋势、季节、节假日分解图
plt.show()
8.3 自定义季节性与节假日
Prophet 的强大之处在于可以灵活添加自定义的季节性和节假日效应。例如，某个产品每周一至周五销量高而周末低（周季节），同时双十一期间有显著的促销脉冲（自定义节假日）：
from prophet import Prophet
import pandas as pd
# === 自定义季节性 ===
model = Prophet()
# 添加周中自定义季节性
model.add_seasonality(
name='monthly_cycle',
period=30.5, # 周期长度（月周期约30.5天）
fourier_order=5, # 傅里叶阶数（越大越灵活）
prior_scale=10.0
)
# === 自定义节假日 ===
# 准备节假日数据框
promotions = pd.DataFrame({
'holiday': 'double_11', # 节假日名称
'ds': pd.to_datetime([
'2020-11-11', '2021-11-11',
'2022-11-11', '2023-11-11'
]),
'lower_window': -3, # 节前3天开始影响
'upper_window': 3, # 节后3天结束影响
})
# 可以加入多个节假日
holidays = pd.concat([
promotions,
pd.DataFrame({
'holiday': 'new_year_sale',
'ds': pd.to_datetime(['2021-01-01', '2022-01-01', '2023-01-01']),
'lower_window': -7,
'upper_window': 7,
})
])
# 将自定义节假日传入模型
model = Prophet(holidays=holidays)
model.add_seasonality('weekly', period=7, fourier_order=3)
model.add_seasonality('yearly', period=365.25, fourier_order=10)
model.fit(df)
# 未来预测
future = model.make_future_dataframe(periods=90)
forecast = model.predict(future)
# 查看节假日效应
holiday_effects = forecast[[
'ds', 'double_11', 'new_year_sale'
]]
print("节假日效应示例:")
print(holiday_effects.head(10))
8.4 Changepoint 检测与调整
Prophet 会自动检测趋势中的变化点（Changepoints）。当时间序列在某个时间点发生结构性变化（如业务策略调整、外部政策变化、竞争环境突变）时，变化点检测可以帮助模型更好地拟合历史数据。通过调整 changepoint_prior_scale 参数可以控制检测的敏感度：
- changepoint_prior_scale=0.001：几乎不检测变化点，趋势平滑（适合长期稳定序列）
- changepoint_prior_scale=0.05：默认值，适中的变化点检测
- changepoint_prior_scale=0.5：非常灵敏，可能过拟合（适合多变的业务数据）
# 手动指定变化点位置
changepoints = pd.to_datetime([
'2021-06-01', '2022-01-01', '2022-09-01'
])
model = Prophet(
changepoints=changepoints, # 手动指定变化点
changepoint_prior_scale=0.1, # 变化点先验规模
changepoint_range=0.9, # 在历史数据的前90%中检测变化点
)
model.fit(df)
future = model.make_future_dataframe(periods=180)
forecast = model.predict(future)
# 查看趋势中的变化点
from prophet.plot import add_changepoints_to_plot
fig = model.plot(forecast)
add_changepoints_to_plot(
fig.gca(), model, forecast
)
plt.show()
Prophet vs SARIMA 选型指南
- 数据量： Prophet 在较少的数据上（如半年到一年）也能获得合理结果；SARIMA 通常需要至少2-3个完整周期
- 多季节性： Prophet 原生支持年、周、日等任意多周期；SARIMA 只能处理单一固定周期
- 节假日效应： Prophet 支持自定义窗口长度的节假日；SARIMA 需要借助外生变量（SARIMAX）
- 可解释性： Prophet 的趋势/季节/节假日分解非常直观；SARIMA 的参数（p,d,q）需要专业知识来解释
- 计算速度： SARIMA（特别是 auto_arima）通常更快；Prophet 随着数据量增大计算变慢
- 模型精度： 对于有固定周期的干净数据，SARIMA 往往更精确；对于多周期、有节假日效应的业务数据，Prophet 更胜一筹
九、时间序列模型评估
模型评估是时间序列预测中至关重要的环节。与普通机器学习不同，时间序列评估必须尊重时间顺序，不能使用随机交叉验证。选择恰当的评估指标和验证策略直接决定了模型在真实场景中的表现。
9.1 核心评估指标
| 指标 |
公式 |
说明 |
优点 |
缺点 |
| RMSE |
sqrt(mean((y - ŷ)²)) |
均方根误差 |
对大误差惩罚大，量纲与y相同 |
对异常值敏感 |
| MAE |
mean(|y - ŷ|) |
平均绝对误差 |
直观，鲁棒 |
无法区分欠/过预测 |
| MAPE |
mean(|(y-ŷ)/y|) x 100% |
平均绝对百分比误差 |
无量纲，可跨序列比较 |
y接近0时无穷大 |
| MASE |
MAE / MAE_naive |
平均绝对缩放误差 |
可比性强，不受量纲影响 |
需要基准模型 |
import numpy as np
import pandas as pd
from sklearn.metrics import mean_absolute_error, mean_squared_error
def evaluate_forecast(y_true, y_pred, y_train=None):
"""
全面评估时间序列预测结果
参数:
- y_true: 真实值
- y_pred: 预测值
- y_train: 训练数据（用于计算MASE）
返回: 包含所有评估指标的字典
"""
# RMSE
rmse = np.sqrt(mean_squared_error(y_true, y_pred))
# MAE
mae = mean_absolute_error(y_true, y_pred)
# MAPE（避免除以零）
nonzero_mask = y_true != 0
mape = np.mean(
np.abs((y_true[nonzero_mask] - y_pred[nonzero_mask])
/ y_true[nonzero_mask])
) * 100
# MASE（需要训练数据的naive预测MAE作为基准）
if y_train is not None:
naive_forecast = y_train[:-1] # naive: 用上一个值预测下一个
naive_actual = y_train[1:]
mae_naive = mean_absolute_error(naive_actual, naive_forecast)
mase = mae / mae_naive if mae_naive > 0 else np.inf
else:
mase = None
return {
'RMSE': round(rmse, 4),
'MAE': round(mae, 4),
'MAPE': round(mape, 2),
'MASE': round(mase, 4) if mase is not None else 'N/A'
}
# 使用示例
np.random.seed(42)
y_true = np.array([100, 102, 105, 103, 108, 110])
y_pred = np.array([101, 103, 104, 104, 107, 109])
y_train = np.random.normal(100, 5, 50)
metrics = evaluate_forecast(y_true, y_pred, y_train)
for k, v in metrics.items():
print(f"{k}: {v}")
9.2 滚动预测验证（Time Series Cross-Validation）
标准的 k 折交叉验证会破坏时间序列的时间顺序（使用未来数据预测过去），导致评估结果过于乐观。时间序列的交叉验证必须使用前向逐步验证：训练集只包含历史数据，测试集总是来自更晚的时间点。
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error
import pandas as pd
import numpy as np
def time_series_cv(series, model_order, n_splits=5, test_size=12):
"""
时间序列滚动交叉验证
参数:
- series: 完整时间序列
- model_order: ARIMA参数 (p,d,q)
- n_splits: 验证轮数
- test_size: 每轮测试集大小
返回: 各轮评估结果
"""
n = len(series)
results = []
for i in range(n_splits):
# 计算训练集和测试集的分割点
train_end = n - (n_splits - i) * test_size
test_end = train_end + test_size
if train_end <= 0 or test_end > n:
continue
train = series[:train_end]
test = series[train_end:test_end]
# 训练模型
model = ARIMA(train, order=model_order)
fitted = model.fit()
# 预测
forecast = fitted.forecast(steps=test_size)
# 评估
mae = mean_absolute_error(test, forecast[:len(test)])
rmse = np.sqrt(np.mean((test.values - forecast[:len(test)])**2))
results.append({
'fold': i + 1,
'train': f'[{train.index[0]}, {train.index[-1]}]',
'test': f'[{test.index[0]}, {test.index[-1]}]',
'MAE': round(mae, 4),
'RMSE': round(rmse, 4)
})
return pd.DataFrame(results)
# 使用示例
np.random.seed(42)
ts_data = pd.Series(
100 + np.cumsum(np.random.normal(0, 1, 200)),
index=pd.date_range('2020-01-01', periods=200, freq='M')
)
cv_results = time_series_cv(ts_data, (1, 1, 1), n_splits=5, test_size=12)
print("滚动交叉验证结果:")
print(cv_results)
print(f"\n平均MAE: {cv_results['MAE'].mean():.4f}")
print(f"平均RMSE: {cv_results['RMSE'].mean():.4f}")
9.3 扩展滚动预测（Expanding Window）
与固定窗口滚动不同，扩展窗口验证中训练集始终从起点开始，随着时间推进不断加入更多历史数据。这种方式更贴近真实的应用场景：模型始终使用"迄今为止的全部数据"来预测未来。
def expanding_window_cv(series, model_order, min_train=50, step=12):
"""
扩展窗口验证
参数:
- series: 时间序列
- model_order: ARIMA (p,d,q)
- min_train: 最小训练样本数
- step: 每次扩展的步长
返回: 各轮评估结果
"""
n = len(series)
results = []
for train_end in range(min_train, n - step + 1, step):
train = series[:train_end]
test = series[train_end:train_end + step]
model = ARIMA(train, order=model_order)
fitted = model.fit()
forecast = fitted.forecast(steps=len(test))
mae = mean_absolute_error(test, forecast)
rmse = np.sqrt(np.mean((test.values - forecast)**2))
results.append({
'train_size': train_end,
'MAE': round(mae, 4),
'RMSE': round(rmse, 4)
})
return pd.DataFrame(results)
# 扩展窗口评估
exp_results = expanding_window_cv(
ts_data, (1, 1, 1),
min_train=50, step=24
)
print("扩展窗口验证结果:")
print(exp_results)
print(f"\n最终MAE: {exp_results['MAE'].iloc[-1]:.4f} （基于更多数据的评估）")
9.4 残差诊断与模型比较
评估模型时，除了预测精度，还需要检查残差是否满足白噪声假设。好的模型残差应该是均值为零、无自相关、近似正态分布的。以下是常用的残差诊断方法：
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy import stats
def residual_diagnostics(residuals, title=''):
"""
残差诊断综合函数
检查: 白噪声检验、正态性检验、自相关
"""
print(f"=== 残差诊断: {title} ===")
# 1. Ljung-Box白噪声检验
lb_test = acorr_ljungbox(
residuals, lags=[10, 20, 30],
return_df=True
)
print("\nLjung-Box白噪声检验:")
print(lb_test)
# 2. 正态性检验（Shapiro-Wilk）
stat, p_value = stats.shapiro(residuals.dropna())
print(f"\nShapiro-Wilk正态性检验:")
print(f" 统计量: {stat:.4f}, p-value: {p_value:.4f}")
print(f" 结论: {'服从正态分布' if p_value > 0.05 else '不服从正态分布'}")
# 3. 残差基本统计
print(f"\n残差基本统计:")
print(f" 均值: {residuals.mean():.4f}")
print(f" 标准差: {residuals.std():.4f}")
print(f" 偏度: {residuals.skew():.4f}")
print(f" 峰度: {residuals.kurtosis():.4f}")
# 4. 如果残差异方差，说明模型未完全捕捉到结构
squared_resid = residuals**2
arch_test = acorr_ljungbox(
squared_resid.dropna(), lags=[10],
return_df=True
)
print(f"\nARCH效应检验（残差平方的自相关检验）:")
if arch_test['lb_pvalue'].iloc[0] < 0.05:
print(" 存在ARCH效应（波动率聚类），考虑GARCH模型")
else:
print(" 无ARCH效应，模型充分捕捉了波动结构")
return lb_test
# 对ARIMA残差进行诊断
model = ARIMA(ts_data, order=(1, 1, 1))
fitted = model.fit()
residual_diagnostics(fitted.resid, 'ARIMA(1,1,1)')
十、核心要点总结
- 分解为先： 时间序列分析的第一步始终是分解（经典分解或STL），理解趋势、季节性、残差的结构是后续所有建模的基础。
- 平稳性是门槛： 使用 ADF + KPSS 联合检验判断平稳性，通过差分（普通差分 + 季节差分）使序列满足建模前提。d 和 D 的确定是 ARIMA/SARIMA 建模的关键。
- ACF/PACF 是罗盘： ACF 截尾定 MA 阶数 q，PACF 截尾定 AR 阶数 p。两者都拖尾时尝试低阶 ARMA 并用 AIC 辅助选择。
- auto_arima 是效率利器： 手动定阶耗时且有主观性，使用 pmdarima 的 auto_arima 自动搜索，同时用逐步搜索（stepwise=True）加速。
- SARIMA 处理固定周期： 月度数据 m=12，季度数据 m=4。SARIMAX 可加入外生变量。当周期很长（如 m=365）或有多重周期时，慎用 SARIMA，改选 Prophet。
- 异常检测从波动率入手： 滚动标准差是最直观的波动率估计，CUSUM 可有效检测均值偏移，两者结合适用于大多数工业异常检测场景。
- Prophet 专为业务预测设计： 自动处理多周期性、节假日效应、缺失值，趋势变化点自适应。在业务预测中比 SARIMA 更易用，尤其适合非专业用户。
- 评估必须尊重时间顺序： 使用滚动预测验证（固定窗口或扩展窗口）代替随机交叉验证。RMSE 和 MAE 是最常用的指标，MASE 提供跨模型可比性。每次评估后务必做残差诊断。
十一、进一步思考与实践
时间序列分析是一个理论与实践紧密结合的领域。以下是一些值得进一步探索的方向和实践建议：
11.1 经典方法 vs. 深度学习方法
近年来，基于深度学习的时序模型（如 LSTM、Transformer、TimesNet、PatchTST）在长期预测任务上表现出色。然而，深度学习并非万能：在数据量较小（少于数千个时间点）时，ARIMA/SARIMA 通常优于神经网络；当预测周期较短（1-3步）时，简单的滚动平均或指数平滑往往已经足够。实践中建议从简单模型开始（如 Naive、指数平滑、ARIMA），只有当性能不满足要求时才尝试更复杂的模型。
11.2 实战建议
- 始终将数据分为训练集、验证集、测试集（按时间顺序），避免数据泄露
- 对原始数据进行可视化探索：趋势是线性还是非线性？季节模式是固定的还是随时间变化的？是否有异常点？
- 当序列受外部因素影响强烈（如促销、政策、突发事件）时，优先考虑 Prophet 或 SARIMAX 这类可引入外部信息的模型
- 不要追求在测试集上绝对精确，关注模型的稳定性和泛化能力。一个好的模型在滚动验证中表现应保持稳定
- 考虑模型集成（Ensemble）：将 ARIMA 的线性捕捉能力与 Prophet 的多周期适应性结合起来，可以有效降低预测方差
- 建立预测的监控和回测机制，定期评估模型在实际环境中的表现，当预测误差持续放大时及时重新训练
11.3 常用工具速查
# ========== 工具速查表 ==========
# 1. 分解
from statsmodels.tsa.seasonal import seasonal_decompose, STL
# 2. 平稳性检验
from statsmodels.tsa.stattools import adfuller, kpss
# 3. ACF/PACF
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# 4. ARIMA
from statsmodels.tsa.arima.model import ARIMA
# 5. SARIMA/SARIMAX
from statsmodels.tsa.statespace.sarimax import SARIMAX
# 6. 自动定阶（需安装pmdarima）
from pmdarima import auto_arima
# 7. Prophet（需安装prophet）
from prophet import Prophet
# 8. 残差诊断
from statsmodels.stats.diagnostic import acorr_ljungbox
# 9. 评估
from sklearn.metrics import mean_absolute_error, mean_squared_error
# 10. 数据清洗与可视化
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt