---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.0
kernelspec:
  name: python3
  display_name: Python 3
---

# 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.

```{code-cell} python
:tags: [hide-input]

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

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

```{code-cell} python
:tags: [hide-input]

# 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)")
```

## 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é

```{code-cell} python
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)')
```

```{code-cell} python
:tags: [hide-input]

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

## 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 |

```{code-cell} python
:tags: [hide-input]

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

## 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.

```{code-cell} python
:tags: [hide-input]

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

## 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.

```{code-cell} python
# 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))
```

```{code-cell} python
# 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())
```

## 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$$

```{code-cell} python
# 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])
```

```{code-cell} python
:tags: [hide-input]

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

## Prévision avec intervalles de prédiction

```{code-cell} python
# 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}%")
```

```{code-cell} python
:tags: [hide-input]

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

## 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é.

```{code-cell} python
# 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}")
```

```{code-cell} python
:tags: [hide-input]

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

## Métriques d'évaluation

```{code-cell} python
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")
```

```{admonition} Choix entre ARIMA et ETS
:class: note

- **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

```{code-cell} python
:tags: [hide-input]

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}")
```

## Résumé

```{admonition} Points clés — Séries temporelles
:class: note

**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.
```