使用 LOESS 进行季节与趋势分解(STL)

本笔记本演示用 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_jump
  • trend_jump
  • low_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。

© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发

请登录后发表评论

    暂无评论内容