Seasonal-Trend decomposition using LOESS (STL)¶

This note book illustrates the use of STL to decompose a time series into three components: trend, season(al) and residual. STL uses LOESS (locally estimated scatterplot smoothing) to extract smooths estimates of the three components. The key inputs into STL are:

season - The length of the seasonal smoother. Must be odd.
trend - The length of the trend smoother, usually around 150% of season. Must be odd and larger than season.
low_pass - The length of the low-pass estimation window, usually the smallest odd number larger than the periodicity of the data.

First we import the required packages, prepare the graphics environment, and prepare the data.

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

Atmospheric CO2¶

The example in Cleveland, Cleveland, McRae, and Terpenning (1990) uses CO2 data, which is in the list below. This monthly data (January 1959 to December 1987) has a clear trend and seasonality across the sample.

[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

The decomposition requires 1 input, the data series. If the data series does not have a frequency, then you must also specify period. The default value for seasonal is 7, and so should also be changed in most applications.

[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 Fitting¶

Setting robust uses a data-dependent weighting function that re-weights data when estimating the LOESS (and so is using LOWESS). Using robust estimation allows the model to tolerate larger errors that are visible on the bottom plot.

Here we use a series the measures the production of electrical equipment in the EU.

[5]:

from statsmodels.datasets import elec_equip as ds

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

Next, we estimate the model with and without robust weighting. The difference is minor and is most pronounced during the financial crisis of 2008. The non-robust estimate places equal weights on all observations and so produces smaller errors, on average. The weights vary between 0 and 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 degree¶

The default configuration estimates the LOESS model with both a constant and a trend. This can be changed to only include a constant by setting COMPONENT_deg to 0. Here the degree makes little difference except in the trend around the financial crisis of 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

Performance¶

Three options can be used to reduce the computational cost of the STL decomposition:

seasonal_jump
trend_jump
low_pass_jump

When these are non-zero, the LOESS for component COMPONENT is only estimated ever COMPONENT_jump observations, and linear interpolation is used between points. These values should not normally be more than 10-20% of the size of seasonal, trend or low_pass, respectively.

The example below shows how these can reduce the computational cost by a factor of 15 using simulated data with both a low-frequency cosinusoidal trend and a sinusoidal seasonal pattern.

[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

First, the base line model is estimated with all jumps equal to 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

The jumps are all set to 15% of their window length. Limited linear interpolation makes little difference to the fit of the model.

[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

Forecasting with STL¶

STLForecast simplifies the process of using STL to remove seasonalities and then using a standard time-series model to forecast the trend and cyclical components.

Here we use STL to handle the seasonality and then an ARIMA(1,1,0) to model the deseasonalized data. The seasonal component is forecast from the find full cycle where

\[E[S_{T+h}|\mathcal{F}_T]=\hat{S}_{T-k}\]

where \(k= m - h + m \lfloor \frac{h-1}{m} \rfloor\). The forecast automatically adds the seasonal component forecast to the ARIMA forecast.

[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 contains information about both the time-series model and the STL decomposition.

[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).

Last update: Oct 03, 2024