Séries temporelles#

Une série temporelle est une donnée qui se souvient de son passé.

Introduction#

Une série temporelle est une suite d’observations ordonnées dans le temps : \(y_1, y_2, \ldots, y_T\). On la rencontre partout — ventes mensuelles, cours boursiers, température quotidienne, trafic web par heure. L’ordre temporel crée une dépendance entre observations qui invalide l’hypothèse d’indépendance classique et exige des méthodes spécifiques.

Les objectifs sont multiples : décrire la structure (tendance, saisonnalité), comprendre les mécanismes générateurs, et surtout prévoir les valeurs futures.

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import pandas as pd
from scipy import stats
from statsmodels.tsa.stattools import adfuller, kpss, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings('ignore')

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)

# Génération d'une série temporelle réaliste : ventes mensuelles
np.random.seed(2024)
T = 120  # 10 ans de données mensuelles
t = np.arange(T)
dates = pd.date_range('2015-01', periods=T, freq='ME')

# Composantes
tendance = 100 + 1.5 * t + 0.005 * t**2
saisonnalite = 30 * np.sin(2 * np.pi * t / 12 - np.pi/2)  # creux en hiver
bruit = np.random.normal(0, 12, T)
ventes = tendance + saisonnalite + bruit
ventes = np.maximum(ventes, 10)  # pas de ventes négatives

serie = pd.Series(ventes, index=dates, name='ventes')

# Vue d'ensemble
fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)

axes[0].plot(dates, ventes, 'steelblue', lw=1.5)
axes[0].set_title('Série temporelle : Ventes mensuelles simulées')
axes[0].set_ylabel('Ventes (unités)')

axes[1].plot(dates, tendance, 'tomato', lw=2, label='Tendance')
axes[1].set_ylabel('Tendance'); axes[1].legend(fontsize=9)

axes[2].plot(dates, saisonnalite, 'darkgreen', lw=1.5, label='Saisonnalité')
axes[2].axhline(0, color='gray', lw=0.8, ls='--')
axes[2].set_ylabel('Saisonnalité'); axes[2].legend(fontsize=9)

axes[3].plot(dates, bruit, 'gray', lw=1, label='Résidu (bruit)')
axes[3].axhline(0, color='gray', lw=0.8, ls='--')
axes[3].set_ylabel('Résidu'); axes[3].legend(fontsize=9)
axes[3].set_xlabel('Date')

plt.suptitle('Décomposition additive d\'une série temporelle', fontsize=12,
             fontweight='bold', y=1.01)
plt.tight_layout()
plt.savefig('_static/15_decomposition_intro.png', dpi=120, bbox_inches='tight')
plt.show()
_images/f6e21788ebfc80a215f6bb7da9cc865789df396dc09de2eb6c8c8c07a6ab3074.png

Décomposition d’une série temporelle#

Toute série peut être vue comme la superposition de composantes :

Modèle additif : \(y_t = T_t + S_t + R_t\)

Modèle multiplicatif : \(y_t = T_t \times S_t \times R_t\)

Le modèle multiplicatif est adapté quand l’amplitude de la saisonnalité croît avec le niveau (typique des séries financières ou de ventes en croissance forte).

Hide code cell source

# Décomposition STL (Seasonal and Trend decomposition using Loess)
stl = STL(serie, period=12, robust=True)
result_stl = stl.fit()

fig = result_stl.plot()
fig.set_size_inches(14, 8)
fig.axes[0].set_title('Décomposition STL — Ventes mensuelles', fontsize=12)
plt.tight_layout()
plt.savefig('_static/15_stl.png', dpi=120, bbox_inches='tight')
plt.show()

print("Statistiques de la décomposition STL :")
print(f"  Force de la tendance    : {1 - result_stl.resid.var()/(result_stl.resid + result_stl.trend).var():.3f}")
print(f"  Force de la saisonnalité: {1 - result_stl.resid.var()/(result_stl.resid + result_stl.seasonal).var():.3f}")
print("  (0=absent, 1=composante dominante)")
_images/da6eb8fecb6d1a41f52b41203be8ec10e6936b39cb38f42914ff1870ba175833.png
Statistiques de la décomposition STL :
  Force de la tendance    : 0.977
  Force de la saisonnalité: 0.798
  (0=absent, 1=composante dominante)

Stationnarité#

Une série est stationnaire au sens faible si :

  1. \(E[y_t] = \mu\) (constante)

  2. \(\text{Var}(y_t) = \sigma^2\) (constante)

  3. \(\text{Cov}(y_t, y_{t+k}) = \gamma_k\) (dépend seulement du lag \(k\), pas de \(t\))

La stationnarité est une condition nécessaire pour la plupart des modèles ARIMA.

Tests de stationnarité#

def test_stationnarite(serie, nom='Série'):
    """Applique les tests ADF et KPSS."""
    print(f"\n=== Tests de stationnarité : {nom} ===")

    # Test ADF : H0 = racine unitaire (non stationnaire)
    adf_result = adfuller(serie.dropna(), autolag='AIC')
    print(f"\nTest ADF (Augmented Dickey-Fuller) :")
    print(f"  Statistique : {adf_result[0]:.4f}")
    print(f"  p-valeur    : {adf_result[1]:.4f}")
    print(f"  Valeurs critiques : {adf_result[4]}")
    if adf_result[1] < 0.05:
        print("  → Rejet H0 : la série est STATIONNAIRE (ADF)")
    else:
        print("  → Non-rejet H0 : la série est probablement NON-STATIONNAIRE (ADF)")

    # Test KPSS : H0 = stationnarité
    kpss_result = kpss(serie.dropna(), regression='c', nlags='auto')
    print(f"\nTest KPSS :")
    print(f"  Statistique : {kpss_result[0]:.4f}")
    print(f"  p-valeur    : {kpss_result[1]:.4f}")
    if kpss_result[1] > 0.05:
        print("  → Non-rejet H0 : la série est STATIONNAIRE (KPSS)")
    else:
        print("  → Rejet H0 : la série est NON-STATIONNAIRE (KPSS)")

# Série originale (avec tendance)
test_stationnarite(serie, 'Ventes originales')

# Après différenciation d'ordre 1
serie_diff = serie.diff().dropna()
test_stationnarite(serie_diff, 'Ventes différenciées (d=1)')
=== Tests de stationnarité : Ventes originales ===

Test ADF (Augmented Dickey-Fuller) :
  Statistique : 2.3637
  p-valeur    : 0.9990
  Valeurs critiques : {'1%': np.float64(-3.4924012594942333), '5%': np.float64(-2.8886968193364835), '10%': np.float64(-2.5812552709190673)}
  → Non-rejet H0 : la série est probablement NON-STATIONNAIRE (ADF)

Test KPSS :
  Statistique : 1.7617
  p-valeur    : 0.0100
  → Rejet H0 : la série est NON-STATIONNAIRE (KPSS)

=== Tests de stationnarité : Ventes différenciées (d=1) ===

Test ADF (Augmented Dickey-Fuller) :
  Statistique : -8.4726
  p-valeur    : 0.0000
  Valeurs critiques : {'1%': np.float64(-3.4924012594942333), '5%': np.float64(-2.8886968193364835), '10%': np.float64(-2.5812552709190673)}
  → Rejet H0 : la série est STATIONNAIRE (ADF)

Test KPSS :
  Statistique : 0.0137
  p-valeur    : 0.1000
  → Non-rejet H0 : la série est STATIONNAIRE (KPSS)

Hide code cell source

fig, axes = plt.subplots(2, 2, figsize=(14, 7))

# Série originale
ax = axes[0, 0]
ax.plot(dates, ventes, 'steelblue', lw=1.5)
ax.set_title('Série originale\n(non stationnaire : tendance)')
ax.set_ylabel('Ventes')

# Après différenciation d'ordre 1
ax = axes[0, 1]
ax.plot(dates[1:], serie_diff, 'tomato', lw=1.2)
ax.axhline(0, color='gray', lw=0.8, ls='--')
ax.set_title('Série différenciée (d=1)\n(stationnaire)')
ax.set_ylabel('Δ Ventes')

# Série à transformer : log pour variance croissante
ventes_log = np.log(ventes)
ventes_log_diff = np.diff(ventes_log)
ax = axes[1, 0]
ax.plot(dates, ventes_log, 'darkgreen', lw=1.5)
ax.set_title('log(Ventes)\n(stabilise la variance)')
ax.set_ylabel('log(Ventes)')

# Log-différencié
ax = axes[1, 1]
ax.plot(dates[1:], ventes_log_diff, 'darkgreen', lw=1.2)
ax.axhline(0, color='gray', lw=0.8, ls='--')
ax.set_title('Δlog(Ventes) = taux de croissance\n(stationnaire)')
ax.set_ylabel('Δlog(Ventes)')

plt.tight_layout()
plt.savefig('_static/15_stationnarite.png', dpi=120, bbox_inches='tight')
plt.show()
_images/30e2b57e62de55b3273c5cd8ce9e72d11aeeed0a5fd7f9129b05de89e9a3fcbc.png

ACF et PACF : identifier l’ordre du modèle#

  • ACF (Autocorrelation Function) : corrélation entre \(y_t\) et \(y_{t-k}\).

  • PACF (Partial ACF) : corrélation entre \(y_t\) et \(y_{t-k}\) après avoir éliminé l’effet des lags intermédiaires.

Règles d’identification#

Modèle

ACF

PACF

AR(p)

Décroissance géométrique ou oscillante

Coupure nette après le lag p

MA(q)

Coupure nette après le lag q

Décroissance géométrique ou oscillante

ARMA(p,q)

Décroissance après lag q

Décroissance après lag p

Hide code cell source

# ACF et PACF sur la série différenciée (série stationnaire)
fig, axes = plt.subplots(2, 2, figsize=(14, 8))

# ACF série originale
plot_acf(serie, lags=40, ax=axes[0, 0], color='steelblue', alpha=0.05)
axes[0, 0].set_title('ACF — Ventes originales\n(autocorrélation lente = non stationnaire)')

# PACF série originale
plot_pacf(serie, lags=40, ax=axes[0, 1], color='steelblue', alpha=0.05, method='ywm')
axes[0, 1].set_title('PACF — Ventes originales')

# ACF série différenciée
plot_acf(serie_diff, lags=40, ax=axes[1, 0], color='tomato', alpha=0.05)
axes[1, 0].set_title('ACF — Ventes différenciées\n(pics aux lags 12, 24 → saisonnalité)')

# PACF série différenciée
plot_pacf(serie_diff, lags=40, ax=axes[1, 1], color='tomato', alpha=0.05, method='ywm')
axes[1, 1].set_title('PACF — Ventes différenciées')

plt.tight_layout()
plt.savefig('_static/15_acf_pacf.png', dpi=120, bbox_inches='tight')
plt.show()
_images/aa500da11197c28d7fc49a28da020ef6e8d620e1dac6e8d59ef29607fffa6219.png

Modèles AR, MA et ARMA#

Processus AR(p) — Autorégressif#

\[y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \varepsilon_t\]

Processus MA(q) — Moyenne mobile#

\[y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}\]

Processus ARMA(p, q)#

Combinaison des deux.

Hide code cell source

# Simulation d'AR(1), AR(2), MA(1), MA(2)
np.random.seed(42)
n_sim = 200
eps = np.random.normal(0, 1, n_sim + 100)

def simulate_ar(phi_list, n=200):
    p = len(phi_list)
    y = np.zeros(n + 100)
    eps_sim = np.random.normal(0, 1, n + 100)
    for t in range(p, n + 100):
        y[t] = sum(phi_list[k] * y[t-k-1] for k in range(p)) + eps_sim[t]
    return y[100:]

def simulate_ma(theta_list, n=200):
    q = len(theta_list)
    y = np.zeros(n + 100)
    eps_sim = np.random.normal(0, 1, n + 100)
    for t in range(q, n + 100):
        y[t] = eps_sim[t] + sum(theta_list[k] * eps_sim[t-k-1] for k in range(q))
    return y[100:]

t_range = range(n_sim)
series_sim = {
    'AR(1) φ=0.8': simulate_ar([0.8]),
    'AR(2) φ₁=1.2, φ₂=-0.6': simulate_ar([1.2, -0.6]),
    'MA(1) θ=0.8': simulate_ma([0.8]),
    'MA(2) θ₁=0.8, θ₂=-0.5': simulate_ma([0.8, -0.5]),
}

fig, axes = plt.subplots(4, 3, figsize=(15, 12))

for row, (label, y_sim) in enumerate(series_sim.items()):
    # Série
    axes[row, 0].plot(t_range, y_sim, 'steelblue', lw=0.8)
    axes[row, 0].set_title(label); axes[row, 0].set_ylabel('y_t')

    # ACF
    acf_vals = acf(y_sim, nlags=25, fft=True)
    axes[row, 1].bar(range(26), acf_vals, color='steelblue', alpha=0.8, width=0.7)
    ci = 1.96 / np.sqrt(n_sim)
    axes[row, 1].axhline(ci, color='tomato', ls='--', lw=1)
    axes[row, 1].axhline(-ci, color='tomato', ls='--', lw=1)
    axes[row, 1].axhline(0, color='black', lw=0.5)
    axes[row, 1].set_title(f'ACF — {label}'); axes[row, 1].set_xlabel('Lag')

    # PACF
    pacf_vals = pacf(y_sim, nlags=25, method='ywm')
    axes[row, 2].bar(range(26), pacf_vals, color='tomato', alpha=0.8, width=0.7)
    axes[row, 2].axhline(ci, color='steelblue', ls='--', lw=1)
    axes[row, 2].axhline(-ci, color='steelblue', ls='--', lw=1)
    axes[row, 2].axhline(0, color='black', lw=0.5)
    axes[row, 2].set_title(f'PACF — {label}'); axes[row, 2].set_xlabel('Lag')

plt.suptitle('Reconnaissance des modèles AR et MA\npar les profils ACF/PACF',
             fontsize=13, fontweight='bold', y=1.01)
plt.tight_layout()
plt.savefig('_static/15_ar_ma_patterns.png', dpi=120, bbox_inches='tight')
plt.show()
_images/4d5895f3fc6364924637334cdf3e5453a0b2bbe1b21ced87aed874643d18c3f2.png

ARIMA : identification et estimation#

Le modèle ARIMA(p, d, q) (AutoRegressive Integrated Moving Average) généralise ARMA en permettant la différenciation d’ordre \(d\) pour rendre la série stationnaire.

# Sélection automatique d'ordre par AIC sur la série de ventes différenciée
# On teste quelques modèles ARIMA manuellement

print("Sélection du modèle ARIMA par AIC :")
print("="*50)
aic_results = []
for p in range(0, 4):
    for q in range(0, 4):
        try:
            mod = ARIMA(serie, order=(p, 1, q)).fit()
            aic_results.append({'p': p, 'd': 1, 'q': q, 'AIC': mod.aic, 'BIC': mod.bic})
        except Exception:
            pass

df_aic = pd.DataFrame(aic_results).sort_values('AIC').head(8)
print(df_aic.to_string(index=False))
Sélection du modèle ARIMA par AIC :
==================================================
 p  d  q         AIC         BIC
 2  1  3  995.483412 1012.158153
 2  1  2 1035.355548 1049.251165
 3  1  3 1049.481567 1068.935432
 3  1  2 1060.427522 1077.102263
 0  1  0 1060.733293 1063.512417
 1  1  3 1061.164003 1075.059620
 0  1  3 1061.912076 1073.028570
 2  1  1 1062.093608 1073.210102
# Ajustement du meilleur modèle ARIMA
best_order = tuple(df_aic.iloc[0][['p', 'd', 'q']].astype(int))
print(f"Meilleur modèle : ARIMA{best_order}")

model_arima = ARIMA(serie, order=best_order).fit()
print(model_arima.summary())
Meilleur modèle : ARIMA(2, 1, 3)
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                 ventes   No. Observations:                  120
Model:                 ARIMA(2, 1, 3)   Log Likelihood                -491.742
Date:                Wed, 01 Apr 2026   AIC                            995.483
Time:                        22:30:59   BIC                           1012.158
Sample:                    01-31-2015   HQIC                          1002.254
                         - 12-31-2024                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          1.7314      0.004    424.047      0.000       1.723       1.739
ar.L2         -0.9999      0.002   -539.558      0.000      -1.004      -0.996
ma.L1         -2.3287      4.741     -0.491      0.623     -11.622       6.964
ma.L2          2.0450      8.387      0.244      0.807     -14.393      18.483
ma.L3         -0.6071      3.340     -0.182      0.856      -7.154       5.940
sigma2       209.4804   1155.390      0.181      0.856   -2055.043    2474.003
===================================================================================
Ljung-Box (L1) (Q):                   5.71   Jarque-Bera (JB):                 3.34
Prob(Q):                              0.02   Prob(JB):                         0.19
Heteroskedasticity (H):               0.82   Skew:                            -0.40
Prob(H) (two-sided):                  0.54   Kurtosis:                         3.16
===================================================================================

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

SARIMA : modèle saisonnier#

Le modèle SARIMA(p,d,q)(P,D,Q)ₛ ajoute des composantes autorégressives et de moyennes mobiles saisonnières d’ordre \(s\) :

\[\Phi_P(B^s) \phi_p(B) \nabla^d \nabla_s^D y_t = \Theta_Q(B^s) \theta_q(B) \varepsilon_t\]
# Modèle SARIMA(1,1,1)(1,1,1)12 pour les ventes mensuelles
model_sarima = SARIMAX(serie,
                       order=(1, 1, 1),
                       seasonal_order=(1, 1, 1, 12),
                       enforce_stationarity=False,
                       enforce_invertibility=False).fit(disp=False)

print(f"SARIMA(1,1,1)(1,1,1)12")
print(f"  AIC : {model_sarima.aic:.1f}")
print(f"  BIC : {model_sarima.bic:.1f}")
print(model_sarima.summary().tables[1])
SARIMA(1,1,1)(1,1,1)12
  AIC : 747.6
  BIC : 760.2
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1232      0.142     -0.870      0.384      -0.401       0.154
ma.L1         -0.8796      0.079    -11.083      0.000      -1.035      -0.724
ar.S.L12       0.0127      0.008      1.514      0.130      -0.004       0.029
ma.S.L12      -1.0000      0.170     -5.870      0.000      -1.334      -0.666
sigma2       129.0331      0.001   9.77e+04      0.000     129.031     129.036
==============================================================================

Hide code cell source

# Diagnostic des résidus SARIMA
fig, axes = plt.subplots(2, 2, figsize=(14, 8))

residus = model_sarima.resid[13:]  # après initialisation

# Résidus dans le temps
ax = axes[0, 0]
ax.plot(dates[13:], residus, 'steelblue', lw=1)
ax.axhline(0, color='tomato', ls='--', lw=1.5)
ax.axhline(2*residus.std(), color='orange', ls=':', lw=1, label='±2σ')
ax.axhline(-2*residus.std(), color='orange', ls=':', lw=1)
ax.set_title('Résidus du modèle SARIMA dans le temps')
ax.set_ylabel('Résidu'); ax.legend(fontsize=8)

# Histogramme des résidus
ax = axes[0, 1]
ax.hist(residus, bins=40, density=True, color='steelblue', alpha=0.7, edgecolor='white')
x_range = np.linspace(residus.min(), residus.max(), 300)
ax.plot(x_range, stats.norm.pdf(x_range, residus.mean(), residus.std()),
        'tomato', lw=2.5, label='N(0, σ̂²)')
ax.set_title('Distribution des résidus')
ax.legend(fontsize=9)

# ACF des résidus (doit être proche de 0)
plot_acf(residus, lags=30, ax=axes[1, 0], color='steelblue', alpha=0.05)
axes[1, 0].set_title('ACF des résidus\n(bruit blanc si pas de structure résiduelle)')

# Q-Q plot
ax = axes[1, 1]
qq = stats.probplot(residus, dist='norm')
ax.scatter(qq[0][0], qq[0][1], color='steelblue', s=20, alpha=0.6)
ax.plot(qq[0][0], qq[1][0] * qq[0][0] + qq[1][1], 'tomato', lw=2)
ax.set_title('Q-Q Plot des résidus')
ax.set_xlabel('Quantiles théoriques'); ax.set_ylabel('Quantiles observés')

plt.suptitle('Diagnostic des résidus — SARIMA(1,1,1)(1,1,1)₁₂', fontsize=12,
             fontweight='bold', y=1.01)
plt.tight_layout()
plt.savefig('_static/15_sarima_diag.png', dpi=120, bbox_inches='tight')
plt.show()
_images/6d26a06e49ca51514c3bc59645662d0aad8d03f53b231000164f1fec2e5712a7.png

Prévision avec intervalles de prédiction#

# Prévision 24 mois en avant
horizon = 24
forecast = model_sarima.get_forecast(steps=horizon)
pred_mean = forecast.predicted_mean
pred_ci = forecast.conf_int(alpha=0.05)

# Dates de prévision
dates_forecast = pd.date_range(serie.index[-1] + pd.DateOffset(months=1),
                                periods=horizon, freq='ME')

# Évaluation sur la fin de la série (walk-forward)
# Utiliser les 96 premiers mois pour entraîner, tester sur les 24 derniers
n_train = 96
serie_train = serie.iloc[:n_train]
serie_test = serie.iloc[n_train:]

model_eval = SARIMAX(serie_train,
                     order=(1, 1, 1),
                     seasonal_order=(1, 1, 1, 12),
                     enforce_stationarity=False).fit(disp=False)

forecast_eval = model_eval.get_forecast(steps=len(serie_test))
y_pred_eval = forecast_eval.predicted_mean
y_true_eval = serie_test.values

# Métriques
mae = np.mean(np.abs(y_true_eval - y_pred_eval))
rmse = np.sqrt(np.mean((y_true_eval - y_pred_eval)**2))
mape = np.mean(np.abs((y_true_eval - y_pred_eval) / y_true_eval)) * 100

print("Évaluation de la prévision (test : 24 derniers mois) :")
print(f"  MAE  : {mae:.1f}")
print(f"  RMSE : {rmse:.1f}")
print(f"  MAPE : {mape:.1f}%")
Évaluation de la prévision (test : 24 derniers mois) :
  MAE  : 9.1
  RMSE : 10.6
  MAPE : 2.9%

Hide code cell source

fig, axes = plt.subplots(2, 1, figsize=(15, 10))

# Prévision sur le jeu de test
ax = axes[0]
ax.plot(serie_train.index, serie_train, 'steelblue', lw=1.5, label='Entraînement')
ax.plot(serie_test.index, serie_test, 'darkgreen', lw=1.5, label='Observé (test)')
ax.plot(serie_test.index, y_pred_eval, 'tomato', lw=2, ls='--', label='Prédit')
ci_eval = forecast_eval.conf_int(alpha=0.05)
ax.fill_between(serie_test.index,
                ci_eval.iloc[:, 0], ci_eval.iloc[:, 1],
                alpha=0.2, color='tomato', label='IC 95%')
ax.axvline(serie_train.index[-1], color='gray', ls=':', lw=2, label='Fin entraînement')
ax.set_title(f'Prévision SARIMA — Évaluation\nMAE={mae:.1f}, RMSE={rmse:.1f}, MAPE={mape:.1f}%')
ax.set_ylabel('Ventes'); ax.legend(fontsize=8)

# Prévision future (fan chart)
ax = axes[1]
ax.plot(serie.index[-36:], serie.iloc[-36:], 'steelblue', lw=1.5, label='Observé')
ax.plot(dates_forecast, pred_mean, 'tomato', lw=2, ls='--', label='Prévision')

# Fan chart avec plusieurs niveaux de confiance
for alpha_fc, color_fc in zip([0.2, 0.4, 0.6, 0.8],
                               ['#d73027', '#f46d43', '#fdae61', '#fee090']):
    ci_fc = forecast.conf_int(alpha=alpha_fc)
    label_fc = f'IC {(1-alpha_fc)*100:.0f}%' if alpha_fc == 0.2 else None
    ax.fill_between(dates_forecast,
                    ci_fc.iloc[:, 0], ci_fc.iloc[:, 1],
                    alpha=0.3, color=color_fc, label=label_fc)

ax.axvline(serie.index[-1], color='gray', ls=':', lw=2)
ax.set_title('Prévision future — Fan chart (intervalles de prédiction 20%/60%/80%/95%)')
ax.set_ylabel('Ventes'); ax.legend(fontsize=8)
ax.set_xlabel('Date')

plt.tight_layout()
plt.savefig('_static/15_forecast.png', dpi=120, bbox_inches='tight')
plt.show()
_images/7098908ed82a6d14a65a5438d3467dee0fefb2e89eaab55dcc60dc23865e83d0.png

Lissage exponentiel (ETS)#

Les modèles ETS (Error, Trend, Seasonality) offrent une alternative aux modèles ARIMA via une pondération exponentielle décroissante des observations passées.

Lissage simple (SES)#

\[\hat{y}_{t+1} = \alpha y_t + (1-\alpha) \hat{y}_t, \quad 0 < \alpha < 1\]

Méthode de Holt (double lissage)#

Prend en compte la tendance avec deux équations de mise à jour.

Holt-Winters (triple lissage)#

Prend en compte tendance et saisonnalité.

# Comparaison des méthodes de lissage exponentiel
models_ets = {
    'SES (α=0.3)': ExponentialSmoothing(serie_train, trend=None, seasonal=None).fit(
        smoothing_level=0.3),
    'Holt (tendance)': ExponentialSmoothing(serie_train, trend='add', seasonal=None,
                                             damped_trend=False).fit(optimized=True),
    'Holt-Winters (additif)': ExponentialSmoothing(serie_train, trend='add',
                                                    seasonal='add', seasonal_periods=12,
                                                    damped_trend=True).fit(optimized=True),
    'Holt-Winters (multiplicatif)': ExponentialSmoothing(serie_train, trend='add',
                                                          seasonal='mul', seasonal_periods=12,
                                                          damped_trend=True).fit(optimized=True),
}

print("Comparaison des modèles ETS :")
print(f"{'Modèle':<35} {'AIC':>10} {'RMSE test':>12}")
print('-' * 60)
for nom, mod in models_ets.items():
    preds = mod.forecast(len(serie_test))
    rmse_ets = np.sqrt(np.mean((serie_test.values - preds.values)**2))
    try:
        aic_ets = mod.aic
    except Exception:
        aic_ets = float('nan')
    print(f"{nom:<35} {aic_ets:>10.1f} {rmse_ets:>12.1f}")
Comparaison des modèles ETS :
Modèle                                     AIC    RMSE test
------------------------------------------------------------
SES (α=0.3)                              626.0         49.7
Holt (tendance)                          591.6         43.7
Holt-Winters (additif)                   515.9         10.0
Holt-Winters (multiplicatif)             532.5         14.5

Hide code cell source

fig, ax = plt.subplots(figsize=(15, 5))

ax.plot(serie.index, serie, 'steelblue', lw=1.5, alpha=0.7, label='Observé', zorder=3)

colors_ets = ['tomato', 'darkgreen', 'orange', 'purple']
for (nom, mod), color in zip(models_ets.items(), colors_ets):
    preds = mod.forecast(len(serie_test))
    # Fitted values sur entraînement
    ax.plot(serie_train.index, mod.fittedvalues,
            color=color, lw=1.2, alpha=0.6, ls=':')
    # Prévisions sur test
    ax.plot(serie_test.index, preds, color=color, lw=2, ls='--', label=nom)

ax.axvline(serie_train.index[-1], color='gray', ls=':', lw=1.5)
ax.set_title('Comparaison des modèles ETS — Prévision sur 24 mois')
ax.set_ylabel('Ventes'); ax.legend(fontsize=8, loc='upper left')
ax.set_xlabel('Date')

plt.tight_layout()
plt.savefig('_static/15_ets.png', dpi=120, bbox_inches='tight')
plt.show()
_images/cd4f9812e8f113b1c3121ecfbfc7beaca290c7bd6277e4ee2726db8a4aa33cd1.png

Métriques d’évaluation#

from sklearn.metrics import mean_absolute_error, mean_squared_error

def metriques_prevision(y_true, y_pred, nom=''):
    """Calcule MAE, RMSE, MAPE et MASE."""
    mae_val = mean_absolute_error(y_true, y_pred)
    rmse_val = np.sqrt(mean_squared_error(y_true, y_pred))
    mape_val = np.mean(np.abs((y_true - y_pred) / (y_true + 1e-9))) * 100
    return {'Modèle': nom, 'MAE': mae_val, 'RMSE': rmse_val, 'MAPE (%)': mape_val}

resultats_eval = []
# SARIMA
resultats_eval.append(metriques_prevision(serie_test.values, y_pred_eval.values, 'SARIMA(1,1,1)(1,1,1)₁₂'))

# ETS
for nom, mod in models_ets.items():
    preds = mod.forecast(len(serie_test))
    resultats_eval.append(metriques_prevision(serie_test.values, preds.values, nom))

df_resultats = pd.DataFrame(resultats_eval).set_index('Modèle')
print("Comparaison finale des modèles :")
print(df_resultats.round(2).to_string())
print("\n→ MAPE < 5% : excellente prévision")
print("→ MAPE 5-10% : bonne prévision")
print("→ MAPE > 20% : prévision médiocre")
Comparaison finale des modèles :
                                MAE   RMSE  MAPE (%)
Modèle                                              
SARIMA(1,1,1)(1,1,1)₁₂         9.07  10.60      2.93
SES (α=0.3)                   43.04  49.65     13.04
Holt (tendance)               37.25  43.73     11.21
Holt-Winters (additif)         8.12   9.98      2.70
Holt-Winters (multiplicatif)  11.53  14.49      3.60

→ MAPE < 5% : excellente prévision
→ MAPE 5-10% : bonne prévision
→ MAPE > 20% : prévision médiocre

Choix entre ARIMA et ETS

  • ARIMA/SARIMA : meilleur pour les séries avec des dépendances complexes, adapté quand ACF/PACF montrent une structure claire. Peut intégrer des variables exogènes (SARIMAX).

  • ETS/Holt-Winters : plus intuitif, robuste, souvent très compétitif en pratique. Plus rapide à estimer.

  • En production, tester les deux et sélectionner par validation croisée temporelle (time series cross-validation).

Validation croisée temporelle#

Hide code cell source

fig, ax = plt.subplots(figsize=(14, 5))

n_splits = 5
n_total = len(serie)
n_test_cv = 12  # 12 mois de test à chaque fold
n_min_train = 48  # minimum 4 ans d'entraînement

mae_folds = []
for fold in range(n_splits):
    n_end = n_total - fold * n_test_cv
    n_start = max(0, n_end - n_min_train - n_test_cv)
    train_cv = serie.iloc[n_start:n_end - n_test_cv]
    test_cv = serie.iloc[n_end - n_test_cv:n_end]

    if len(train_cv) < n_min_train:
        continue

    try:
        mod_cv = SARIMAX(train_cv, order=(1,1,1), seasonal_order=(1,1,1,12),
                         enforce_stationarity=False).fit(disp=False)
        pred_cv = mod_cv.forecast(len(test_cv))
        mae_folds.append(mean_absolute_error(test_cv.values, pred_cv.values))

        color_fold = plt.cm.viridis(fold / n_splits)
        ax.plot(train_cv.index, train_cv, color=color_fold, lw=1, alpha=0.5)
        ax.plot(test_cv.index, test_cv, 'steelblue', lw=1.5, alpha=0.8)
        ax.plot(test_cv.index, pred_cv, color=color_fold, lw=2, ls='--')
    except Exception:
        pass

ax.set_title(f'Validation croisée temporelle (Time Series CV)\n'
             f'MAE moyen : {np.mean(mae_folds):.1f} ± {np.std(mae_folds):.1f}')
ax.set_ylabel('Ventes'); ax.set_xlabel('Date')
ax.text(0.02, 0.95, '← Entraînement  | Prévision →',
        transform=ax.transAxes, fontsize=10, verticalalignment='top')

plt.tight_layout()
plt.savefig('_static/15_cv_ts.png', dpi=120, bbox_inches='tight')
plt.show()

print(f"MAE par fold : {[f'{m:.1f}' for m in mae_folds]}")
print(f"MAE moyen : {np.mean(mae_folds):.1f}")
_images/e3dd33afb68aa44045977f6eec9583ed03d8318ac8ca5d22919341fa26331906.png
MAE par fold : ['9.1', '13.5', '12.7', '9.7', '16.8']
MAE moyen : 12.4

Résumé#

Points clés — Séries temporelles

Prétraitement :

  • Vérifier la stationnarité (tests ADF et KPSS) et différencier si nécessaire.

  • Visualiser ACF/PACF pour identifier l’ordre (p, q) du modèle.

  • Décomposer avec STL pour séparer tendance, saisonnalité et résidu.

Modèles :

  • ARIMA(p,d,q) : modèle général pour séries stationnaires après différenciation.

  • SARIMA(p,d,q)(P,D,Q)ₛ : extension saisonnière indispensable pour les données périodiques.

  • ETS/Holt-Winters : approche de lissage, souvent compétitive et plus interprétable.

Validation :

  • Toujours évaluer sur des données futures (jamais d’inversion temporelle !).

  • Métriques : MAE, RMSE, MAPE selon le contexte.

  • Valider les résidus : bruit blanc (ACF ≈ 0, gaussien).

  • Utiliser la validation croisée temporelle pour estimer la généralisation.