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

此笔记本演示了 STL 的使用，它将时间序列分解为三个部分：趋势、季节性和残差。STL 使用 LOESS（局部加权回归）来提取三个部分的平滑估计。 STL 的关键输入是

season - 季节性平滑器的长度。必须为奇数。
trend - 趋势平滑器的长度，通常约为 season 的 150%。必须为奇数且大于 season。
low_pass - 低通估计窗口的长度，通常是大于数据周期性的最小奇数。

首先，我们导入所需的包，准备图形环境并准备数据。

[1]:

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from pandas.plotting import register_matplotlib_converters

register_matplotlib_converters()
sns.set_style("darkgrid")

[2]:

plt.rc("figure", figsize=(16, 12))
plt.rc("font", size=13)

大气 CO2¶

Cleveland、Cleveland、McRae 和 Terpenning（1990 年）的示例使用 CO2 数据，这些数据在下面的列表中。此月度数据（1959 年 1 月至 1987 年 12 月）在整个样本中具有明显的趋势和季节性。

[3]:

co2 = [
58,
39,
79,
82,
39,
22,
68,
01,
02,
55,
02,
75,
52,
10,
79,
22,
08,
70,
27,
99,
24,
05,
05,
23,
92,
76,
54,
49,
64,
85,
70,
96,
17,
47,
19,
17,
12,
72,
79,
68,
28,
89,
79,
56,
46,
59,
85,
87,
87,
25,
13,
49,
34,
62,
85,
87,
36,
24,
13,
46,
57,
23,
89,
54,
20,
90,
42,
60,
73,
15,
94,
91,
73,
78,
23,
49,
59,
35,
61,
24,
23,
76,
36,
50,
35,
40,
22,
45,
80,
50,
16,
09,
26,
66,
47,
70,
06,
23,
78,
10,
63,
79,
34,
73,
00,
99,
41,
68,
30,
89,
59,
65,
30,
15,
88,
80,
99,
86,
88,
36,
59,
23,
34,
33,
03,
24,
39,
16,
87,
31,
34,
74,
61,
58,
55,
81,
82,
53,
29,
66,
12,
09,
01,
10,
12,
62,
16,
94,
15,
79,
53,
65,
60,
78,
13,
26,
93,
84,
96,
93,
25,
24,
13,
42,
97,
29,
56,
73,
73,
70,
46,
70,
66,
22,
02,
39,
58,
27,
30,
81,
44,
89,
62,
85,
29,
44,
35,
58,
58,
55,
56,
73,
45,
98,
63,
88,
63,
53,
90,
08,
59,
31,
44,
64,
62,
45,
36,
46,
84,
29,
04,
88,
23,
83,
18,
50,
80,
22,
54,
82,
45,
97,
65,
40,
28,
73,
05,
54,
65,
06,
32,
39,
66,
56,
24,
39,
43,
22,
61,
78,
88,
43,
61,
53,
06,
92,
39,
72,
64,
65,
07,
53,
82,
19,
89,
56,
22,
92,
26,
27,
66,
54,
71,
79,
79,
06,
93,
02,
65,
80,
01,
94,
17,
28,
76,
05,
18,
04,
16,
01,
64,
91,
72,
52,
75,
68,
14,
37,
32,
45,
05,
91,
77,
30,
98,
41,
89,
03,
19,
87,
74,
55,
28,
00,
37,
74,
36,
19,
97,
20,
76,
96,
82,
82,
24,
09,
66,
90,
27,
21,
88,
58,
99,
31,
98,
72,
63,
24,
83,
10,
52,
43,
48,
89,
29,
54,
66,
07,
12,
55,
34,
80,
10,
54,
20,
20,
44,
67,
]
co2 = pd.Series(
    co2, index=pd.date_range("1-1-1959", periods=len(co2), freq="M"), name="CO2"
)
co2.describe()

/tmp/ipykernel_3881/1071343913.py:352: FutureWarning: 'M' is deprecated and will be removed in a future version, please use 'ME' instead.
  co2, index=pd.date_range("1-1-1959", periods=len(co2), freq="M"), name="CO2"

[3]:

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

分解需要 1 个输入，即数据序列。如果数据序列没有频率，那么你还必须指定 period。 seasonal 的默认值为 7，因此在大多数应用程序中也应该更改。

[4]:

from statsmodels.tsa.seasonal import STL

stl = STL(co2, seasonal=13)
res = stl.fit()
fig = res.plot()

../../../_images/examples_notebooks_generated_stl_decomposition_6_0.png

稳健拟合¶

设置 robust 将使用数据依赖的加权函数，该函数在估计 LOESS 时重新加权数据（因此使用 LOWESS）。使用稳健估计允许模型容忍在底部图中可见的较大误差。

这里我们使用一个测量欧盟电气设备生产量的序列。

[5]:

from statsmodels.datasets import elec_equip as ds

elec_equip = ds.load().data.iloc[:, 0]

接下来，我们使用和不使用稳健加权估计模型。差异很小，在 2008 年金融危机期间最为明显。非稳健估计对所有观察结果赋予相同的权重，因此平均产生较小的误差。权重在 0 到 1 之间变化。

[6]:

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):
        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"])

../../../_images/examples_notebooks_generated_stl_decomposition_10_0.png

[7]:

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])

../../../_images/examples_notebooks_generated_stl_decomposition_11_0.png

LOESS 度数¶

默认配置使用常数和趋势来估计 LOESS 模型。可以通过将 COMPONENT_deg 设置为 0 来更改为仅包含常数。这里，度数几乎没有区别，除了在 2008 年金融危机期间的趋势之外。

[8]:

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"])

../../../_images/examples_notebooks_generated_stl_decomposition_13_0.png

性能¶

可以使用三种选项来降低 STL 分解的计算成本

seasonal_jump
trend_jump
low_pass_jump

当这些值不为零时，组件 COMPONENT 的 LOESS 仅在每 COMPONENT_jump 个观察值进行估计，并且在点之间使用线性插值。这些值通常不应超过 seasonal、trend 或 low_pass 大小的 10-20%。

下面的示例展示了如何使用模拟数据（具有低频余弦趋势和正弦季节模式）将计算成本降低 15 倍。

[9]:

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)

../../../_images/examples_notebooks_generated_stl_decomposition_15_0.png

首先，使用所有跳跃都等于 1 来估计基线模型。

[10]:

mod = STL(y, period=period, seasonal=seasonal)
%timeit mod.fit()
res = mod.fit()
fig = res.plot(observed=False, resid=False)

284 ms ± 38.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

../../../_images/examples_notebooks_generated_stl_decomposition_17_1.png

所有跳跃都设置为其窗口长度的 15%。有限的线性插值对模型的拟合几乎没有影响。

[11]:

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)

22.9 ms ± 1.04 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

../../../_images/examples_notebooks_generated_stl_decomposition_19_1.png

使用 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 预测中。

[12]:

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.forecasting.stl import STLForecast

elec_equip.index.freq = elec_equip.index.inferred_freq
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()

../../../_images/examples_notebooks_generated_stl_decomposition_21_0.png

summary 包含有关时间序列模型和 STL 分解的信息。

[13]:

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, 03 Oct 2024   AIC                           1050.868
Time:                        15:45: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
---------------------------------------------------------------------------------

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

最后更新：2024 年 10 月 3 日