本笔记本演示用 STL 将时间序列分解为趋势、季节和残差三个部分。STL 使用 LOESS(局部估计散点平滑)获得这三个部分的平滑估计。主要输入参数如下:
season:季节平滑器长度,必须为奇数。trend:趋势平滑器长度,通常约为season的 150%,必须为奇数且大于season。low_pass:低通估计窗口长度,通常取大于数据周期的最小奇数。
先导入所需包,配置绘图环境并准备数据。
import matplotlib.pyplot as plt
import pandas as pd
from pandas.plotting import register_matplotlib_converters
import seaborn as sns
register_matplotlib_converters()
sns.set_style("darkgrid")
plt.rc("figure", figsize=(16, 12))
plt.rc("font", size=13)
大气中的 CO₂
Cleveland、Cleveland、McRae 和 Terpenning(1990)的示例使用下面的 CO₂ 数据。这组月度数据覆盖 1959 年 1 月至 1987 年 12 月,整个样本中都有明显的趋势和季节性。
co2 = pd.read_csv(
"https://raw.githubusercontent.com/statsmodels/smdatasets/refs/heads/main/data/stl-decomposition/co2.csv",
parse_dates=True,
index_col=0,
).iloc[:, 0]
co2.describe()
原文数据描述输出:
count 348.000000
mean 330.123879
std 10.059747
min 313.550000
25% 321.302500
50% 328.820000
75% 338.002500
max 351.340000
Name: co2, dtype: float64
分解只需要一个输入:数据序列。如果序列没有频率信息,还必须指定 period。seasonal 的默认值是 7,多数应用中也应调整它。
from statsmodels.tsa.seasonal import STL
stl = STL(co2, seasonal=13)
res = stl.fit()
fig = res.plot()
稳健拟合
设置 robust 会使用依赖数据的权重函数,在估计 LOESS 时重新赋权,因此使用的是 LOWESS。稳健估计允许模型容忍较大的误差;这些误差体现在原文分解图最下方的残差部分。
这里使用欧盟电气设备生产量序列。
from statsmodels.datasets import elec_equip as ds
elec_equip = ds.load().data.iloc[:, 0]
接下来分别使用和不使用稳健权重估计模型。两者差异较小,在 2008 年金融危机期间最明显。非稳健估计对所有观测使用相同权重,因此平均误差较小。稳健权重介于 0 和 1 之间。
def add_stl_plot(fig, res, legend):
"""Add 3 plots from a second STL fit"""
axs = fig.get_axes()
comps = ["trend", "seasonal", "resid"]
for ax, comp in zip(axs[1:], comps, strict=False):
series = getattr(res, comp)
if comp == "resid":
ax.plot(series, marker="o", linestyle="none")
else:
ax.plot(series)
if comp == "trend":
ax.legend(legend, frameon=False)
stl = STL(elec_equip, period=12, robust=True)
res_robust = stl.fit()
fig = res_robust.plot()
res_non_robust = STL(elec_equip, period=12, robust=False).fit()
add_stl_plot(fig, res_non_robust, ["Robust", "Non-robust"])
fig = plt.figure(figsize=(16, 5))
lines = plt.plot(res_robust.weights, marker="o", linestyle="none")
ax = plt.gca()
xlim = ax.set_xlim(elec_equip.index[0], elec_equip.index[-1])
LOESS 的阶数
默认配置的 LOESS 模型同时包含常数和趋势。将 COMPONENT_deg 设为 0,可以改为只包含常数项。在本例中,除 2008 年金融危机附近的趋势外,阶数选择带来的差异不大。
stl = STL(
elec_equip, period=12, seasonal_deg=0, trend_deg=0, low_pass_deg=0, robust=True
)
res_deg_0 = stl.fit()
fig = res_robust.plot()
add_stl_plot(fig, res_deg_0, ["Degree 1", "Degree 0"])
性能
以下三个选项可以降低 STL 分解的计算开销:
seasonal_jumptrend_jumplow_pass_jump
这些参数非零时,某个 COMPONENT 的 LOESS 只会每隔 COMPONENT_jump 个观测估计一次,中间位置用线性插值填充。通常不应超过对应 seasonal、trend 或 low_pass 窗口长度的 10%~20%。
下例采用同时包含低频余弦趋势与正弦季节模式的模拟数据,展示这些选项如何降低计算开销。原文描述约 15 倍加速;实际计时结果随环境和版本变化。
import numpy as np
rs = np.random.RandomState(0xA4FD94BC)
tau = 2000
t = np.arange(tau)
period = int(0.05 * tau)
seasonal = period + ((period % 2) == 0) # Ensure odd
e = 0.25 * rs.standard_normal(tau)
y = np.cos(t / tau * 2 * np.pi) + 0.25 * np.sin(t / period * 2 * np.pi) + e
plt.plot(y)
plt.title("Simulated Data")
xlim = plt.gca().set_xlim(0, tau)
先拟合基准模型,所有 jump 都等于 1。
mod = STL(y, period=period, seasonal=seasonal)
%timeit mod.fit()
res = mod.fit()
fig = res.plot(observed=False, resid=False)
原文运行记录:
342 ms ± 27.6 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
接着把 jump 设为相应窗口长度的 15%。有限的线性插值对模型拟合影响很小。
low_pass_jump = seasonal_jump = int(0.15 * (period + 1))
trend_jump = int(0.15 * 1.5 * (period + 1))
mod = STL(
y,
period=period,
seasonal=seasonal,
seasonal_jump=seasonal_jump,
trend_jump=trend_jump,
low_pass_jump=low_pass_jump,
)
%timeit mod.fit()
res = mod.fit()
fig = res.plot(observed=False, resid=False)
原文运行记录:
32.6 ms ± 371 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
使用 STL 预测
STLForecast 简化了先通过 STL 消除季节性,再使用常规时间序列模型预测趋势与周期成分的过程。
这里用 STL 处理季节性,然后通过 ARIMA(1,1,0) 建模去季节化数据。季节部分根据最后一个完整周期预测:
E[S_{T+h}|\mathcal{F}_T]=\hat{S}_{T-k}
其中 k= m – h + m \lfloor \frac{h-1}{m} \rfloor。最终预测会自动把季节成分预测加回 ARIMA 预测。
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.forecasting.stl import STLForecast
elec_equip.index.freq = "MS"
stlf = STLForecast(elec_equip, ARIMA, model_kwargs=dict(order=(1, 1, 0), trend="t"))
stlf_res = stlf.fit()
forecast = stlf_res.forecast(24)
plt.plot(elec_equip)
plt.plot(forecast)
plt.show()
summary 同时包含时间序列模型与 STL 分解的信息。
print(stlf_res.summary())
原文模型报告:
STL Decomposition and SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 257
Model: ARIMA(1, 1, 0) Log Likelihood -522.434
Date: Thu, 27 Aug 2026 AIC 1050.868
Time: 07:06:57 BIC 1061.504
Sample: 01-01-1995 HQIC 1055.146
- 05-01-2016
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
x1 0.1171 0.118 0.995 0.320 -0.113 0.348
ar.L1 -0.0435 0.049 -0.880 0.379 -0.140 0.053
sigma2 3.4682 0.188 18.406 0.000 3.099 3.837
===================================================================================
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 223.01
Prob(Q): 0.92 Prob(JB): 0.00
Heteroskedasticity (H): 0.33 Skew: -0.26
Prob(H) (two-sided): 0.00 Kurtosis: 7.54
STL Configuration
=================================================================================
Period: 12 Trend Length: 23
Seasonal: 7 Trend deg: 1
Seasonal deg: 1 Trend jump: 1
Seasonal jump: 1 Low pass: 13
Robust: False Low pass deg: 1
Low pass jump: 1
---------------------------------------------------------------------------------
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
原文:Seasonal-Trend decomposition using LOESS (STL),statsmodels 文档与示例贡献者。正文与代码来自同名官方 notebook 源文件,说明文字译为中文,代码保留不变。原文输出是网页已有的运行记录,本次未运行示例;配图仅链接到原文,未作视觉或 HTTP 检查。原文使用 season 描述参数,调用示例使用 seasonal;此处保留原文写法。
Copyright (C) 2006, Jonathan E. Taylor;Copyright (c) 2006–2008 Scipy Developers;Copyright (c) 2009–2018 statsmodels Developers。依据 BSD 三条款许可证改编,完整版权、条款和免责声明见随附 LICENSE-source.txt。文档页版权标注为 © 2009–2025 Josef Perktold、Skipper Seabold、Jonathan Taylor、statsmodels-developers。











暂无评论内容