Décomposition des Séries Temporelles#

Une série temporelle peut être décomposée en trois composants : tendance, saisonnalité et résidu. Plusieurs algorithmes peuvent être employés pour extraire des lissages des estimateurs de ces trois composants d’un jeu de données constituant la réalisation d’une série chronologique, au premier rang duquel \(\rm LOESS\) (Locally Estimated Scatterplot Smoothing).

Cette technique de régression non-paramètrique met en oeuvre une fonction pour capturer la variation dans un nuage de point tout en réduisant le bruit et avec des hypothèses minimales sur les relations entre les variables.

N.B.

LOESS est très flexible, ce qui le rend idéal pour modéliser des processus complexes pour lesquels aucun modèle théorique n’existe. Cependant, elle nécessite des jeux de données assez grands et densément échantillonnés afin de produire de bons modèles.

L’outil offert par la bibliothèque statsmodels pour la décomposition des séries temporelles est la méthode STL, implémentant la régression \(\rm LOESS\). Les paramètres clés de cette fonction sont :

seasonal : longueur du lissage de la composante saisonnière. Doit être obligatoirement impair.
trend : longueur du lissage de la tendance, généralement environ \(150\%\) de seasonal. Doit être obligatoirement impair et supérieur à la longueur de la saisonnalité.
low_pass : longueur de la fenêtre d’estimation passe-bas. Correspond généralement au plus petit nombre impair plus grand que la périodicité de la série.

On commence d’emblée par importer les paquets requis, notamment matplotlib.pyplot, pandas et numpy, puis préparer l’environnement graphique et les données.

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

register_matplotlib_converters()
sns.set_style("darkgrid")
plt.rc("figure", figsize=(20,12))
plt.rc("font", size=14)

Le jeu de données en question est une série des émissions quotidiennes du gaz atmosphérique CO2 à partir de janvier 1959 à décembre 1987, enregistrée dans l’observatoire de Mauna Loa à Hawaii, inspiré de l’ouvrage de Cleveland, McRae, and Terpenning (1990) sur \(\rm LOESS\). Il présente une tendance et une saisonnalité flagrantes sur l’ensemble de l’échantillon.

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

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

La décomposition nécessite un seul argument, la série de données. Si la série n’a pas de fréquence, il faut spécifier la période period. La valeur par défaut de seasonal est 7.

from statsmodels.tsa.seasonal import STL
stl = STL(co2, seasonal=13)
res = stl.fit()

Ajustement Robuste#

Spécifier le paramètre robuste permet d’utiliser une fonction de pondération dépendant des données qui repondère les données lors de l’estimation de la \(\rm LOESS\) (donnant ainsi \(\rm LOWESS\)). L’utilisation de l’estimation robuste permet au modèle de tolérer des erreurs plus importantes comme celles visibles sur le graphique du bas ci-dessus.

Nous utilisons ici une série qui mesure la production d’équipements électriques dans l’UE.

from statsmodels.datasets import elec_equip as ds

elec_equip = ds.load().data

Ensuite, nous estimons le modèle avec et sans pondération robuste. La différence est mineure et est plus prononcée pendant la crise financière de 2008. L’estimation non robuste attribue des poids égaux à toutes les observations et produit donc des erreurs plus faibles, en moyenne. Les pondérations varient entre 0 et 1.

def add_stl_plot(fig, res, legend):
    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, ["Robuste", "Non-Robuste"])

Conclusion

L’ajustement robuste est clairement plus performant, dès lors qu’il reflète mieux la tendance de la série temporelle étant donnée une anomalie dans les observations.

Degré du \(\rm LOESS\)#

La configuration par défaut estime le modèle \(\rm LOESS\) avec à la fois une constante et une tendance. Ceci peut être modifié pour n’inclure qu’une constante en spécifiant COMPONENT_deg à 0. Ici, le degré a peu d’influence, sauf pour la tendance autour de la crise financière de 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, ["Degré 1", "Degré 0"])

Complexité & Performance#

Trois options peuvent être utilisées pour réduire le coût calculatoire de la décomposition STL : seasonal_jump, trend_jump et low_pass_jump.

Lorsque ces valeurs sont différentes de zéro, \(\rm LOESS\) de la composante COMPONENT est estimé uniquement à partir des observations COMPONENT_jump, et une interpolation linéaire est effectuée entre les points. Ces valeurs ne devraient normalement pas être supérieures à \(10-20\%\) de la taille de seasonal, trend ou low_pass, respectivement.

L’exemple ci-dessous montre comment ceux-ci peuvent réduire le coût de calcul d’un facteur de 15 fois en utilisant des données générées avec une tendance cosinusoïdale à basse fréquence et un modèle saisonnier sinusoïdal.

rs = np.random.RandomState(0xA4FD94BC)
tau = 2000
t = np.arange(tau)
period = int(0.05 * tau)
seasonal = period + ((period % 2) == 0)  # Obligatoirement impair!
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("Données Simulées")
xlim = plt.gca().set_xlim(0, tau)

Tout d’abord, le modèle de base est estimé avec tous les sauts égaux à 1.

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

238 ms ± 577 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Les sauts sont tous fixés à \(15\%\) de la longueur de leur fenêtre. L’interpolation linéaire limitée fait peu de différence dans l’ajustement du modèle.

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)

16.7 ms ± 37.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Bilan

Les options susmentionnées améliorent significativement le temps d’exécution de la décomposition, et ce même pour des sauts très faibles.

Prédiction au moyen de `STL`#

STLForecast simplifie le processus de décomposition de la STL pour supprimer les saisonnalités et ensuite utiliser un modèle de série temporelle standard pour prévoir la tendance et les composantes cycliques.

Dans le cas envisagé, on emploie STL pour traiter la saisonnalité, puis un modèle \(\mathrm{ARIMA}(1,1,0)\) pour modéliser les données désaisonnalisées.

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.title("Prédiction de la production des équipements électriques dans l'UE")
plt.show()

summary renvoie des informations à la fois sur le modèle de séries temporelles et la décomposition 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, 05 Jan 2023   AIC                           1050.868
Time:                        23:24:41   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).

Introduction à Statsmodels

Décomposition des Séries Temporelles

Contents

Décomposition des Séries Temporelles#

Ajustement Robuste#

Degré du \(\rm LOESS\)#

Complexité & Performance#

Prédiction au moyen de `STL`#

Introduction à Statsmodels

Décomposition des Séries Temporelles

Contents

Décomposition des Séries Temporelles#

Ajustement Robuste#

Degré du \(\rm LOESS\)#

Complexité & Performance#

Prédiction au moyen de STL#

Prédiction au moyen de `STL`#