Prévision Hebdo Linky : ETS, SARIMA, SARIMAX

Analyse et prévision de séries temporelles Linky avec ETS, SARIMA et SARIMAX.

Introduction

30 Second Summary

Every week, thousands of Linky smart meters are deployed across France for maintenance, network growth, and equipment replacement. Planning teams need accurate forecasts of these installation volumes to schedule crews, order materials, and manage budgets months in advance.

In this project, you will build a complete time series forecasting notebook that explores three weekly Linky deployment streams, compares three statistical models (ETS, SARIMA, SARIMAX), and identifies the best predictor for each stream using rigorous residual diagnostics.

What You'll Build

You'll run your notebook and watch it automatically test dozens of model configurations, diagnose their residuals with statistical tests, and produce weekly forecasts with confidence intervals stretching all the way to 2028.

By the end of this project, you'll have:

  • A full exploratory analysis with seasonal decomposition, stationarity tests (ADF/KPSS), and ACF/PACF plots revealing the hidden structure in each deployment flux.
  • An automated grid search across ETS, SARIMA, and SARIMAX configurations that selects the optimal model per flux based on AIC, RMSE, and residual diagnostics.
  • Production-ready forecasts with confidence intervals through 2028, powered by French calendar effects (public holidays and school vacations) as exogenous variables.
  • Secret Mission: Add temporal cross-validation with an expanding window to stress-test each model's real-world stability.

Are there any prerequisites?

You need Python 3.10+ with Jupyter Notebook, pandas, matplotlib, and statsmodels already installed. You also need your weekly Linky deployment data (Maintenance, Croissance, Resorption) in a CSV file with a WEEK column as the index.

Before We Start

Before we dive in, let's lock in what you're about to build and why it matters.

Set Up Your Environment and Load Data

Les librairies de base comme pandas et statsmodels ne suffisent pas pour ce projet.

Le grid search SARIMA necessite pmdarima, et la modelisation des effets calendaires francais necessite des librairies dediees aux jours feries et vacances scolaires. Dans cette etape, tu vas installer les trois librairies manquantes, importer tous les modules, et preparer les sous-series temporelles de chaque flux Linky.

In this step, get ready to:
  • Installer pmdarima, holidays et vacances-scolaires-france.
  • Importer tous les modules et charger les donnees CSV Linky.
  • Definir les parametres globaux et preparer les sous-series filtrees par flux.
Install the missing libraries

Ton environnement Python dispose deja de pandas, matplotlib et statsmodels. Il manque trois packages pour le grid search automatise et les effets calendaires.

  • Installe les librairies manquantes en executant cette commande dans une cellule de ton notebook :
pip install pmdarima==2.1.1 holidays vacances-scolaires-france

What do these libraries do?

  • pmdarima fournit auto_arima, un grid search automatise pour trouver les meilleurs parametres SARIMA.
  • holidays genere la liste des jours feries francais par annee.
  • vacances-scolaires-france identifie les semaines de vacances scolaires par zone (A, B, C).
  • Verifie que l'installation a fonctionne en executant cette commande dans une nouvelle cellule :
import pmdarima
print(pmdarima.__version__)

Tu devrais voir 2.1.1 s'afficher dans la sortie de la cellule.

Seeing an installation error?

  • Verifie que tu utilises Python 3.10 ou superieur avec python --version.
  • Si tu as un conflit de dependances, essaie de relancer le kernel Jupyter puis reinstaller.
  • help me debug a pip install error for pmdarima
Import all modules and load the data

Tous les modules dont tu auras besoin pour le projet complet sont importes en une seule cellule. Cela inclut les outils de decomposition, les tests statistiques, les modeles ETS/SARIMA/SARIMAX et les librairies calendaires.

  • Cree une nouvelle cellule dans ton notebook et colle le bloc d'imports suivant :
# =============================================================================
# Prevision Hebdo Linky - Projet 1 : EDA + Modeles Statistiques
# =============================================================================

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox, het_arch
from statsmodels.stats.stattools import jarque_bera
import pmdarima as pm
from pmdarima.arima import StepwiseContext
import holidays
from vacances_scolaires_france import SchoolHolidayDates
import datetime

What does this code do?

  • warnings.filterwarnings('ignore') supprime les warnings de convergence qui polluent la sortie pendant les grid searches.
  • Les imports statsmodels couvrent la decomposition saisonniere, les tests de stationnarite (ADF, KPSS), les correlogrammes (ACF/PACF), et les modeles ETS/SARIMAX.
  • StepwiseContext de pmdarima permet de limiter le temps de recherche lors du grid search SARIMA.
  • SchoolHolidayDates identifie les vacances scolaires par zone geographique francaise.
  • Execute la cellule. Si aucune erreur ne s'affiche, tous les modules sont correctement charges.

Maintenant, charge ton fichier CSV Linky dans une nouvelle cellule.

  • Charge les donnees en adaptant le chemin vers ton fichier CSV :
# Charger les donnees (adapter le chemin)
df = pd.read_csv('data_linky.csv', index_col='WEEK', parse_dates=True)
df.index.freq = 'W-MON'
print(df.shape)
df.head()

Tu devrais voir un tuple de la forme (N, 3) suivi des premieres lignes avec les colonnes Maintenance, Croissance et Resorption.

FileNotFoundError or wrong shape?

  • Verifie que le chemin vers ton fichier CSV est correct par rapport a l'emplacement de ton notebook.
  • La colonne d'index doit s'appeler WEEK dans ton CSV. Verifie avec un editeur de texte si necessaire.
  • Si tu vois 1 ou 2 colonnes au lieu de 3, verifie que ton CSV contient bien les colonnes Maintenance, Croissance et Resorption.
  • help me load my Linky CSV data into a pandas DataFrame
Define global parameters and prepare the sub-series

Chaque flux Linky commence a une date differente. Les flux Maintenance et Croissance demarrent en aout 2021, tandis que Resorption ne commence qu'en janvier 2023.

  • Definis les parametres globaux et la configuration par flux en collant ce code dans une nouvelle cellule :
# =============================================================================
# PARAMETRES GLOBAUX
# =============================================================================
last_forecast = '2028-12-31'
first_date_reso = '2023-01-01'
first_date = '2021-08-02'

# Parametres par flux
flux_config = {
    'Maintenance': {'start': first_date},
    'Croissance': {'start': first_date},
    'Resorption': {'start': first_date_reso},
}

What do these parameters define?

  • last_forecast = '2028-12-31' fixe la date de fin des previsions que tu genereras dans les etapes suivantes.
  • first_date et first_date_reso definissent le debut des donnees exploitables pour chaque flux. Resorption demarre plus tard car ce programme n'existait pas avant 2023.
  • flux_config centralise la configuration dans un dictionnaire pour pouvoir iterer sur les 3 flux avec une seule boucle.
  • Dans une nouvelle cellule juste apres, ajoute la boucle de filtrage et execute-la :
# Filtrer les series par date de debut et supprimer les zeros
series_dict = {}
for flux, config in flux_config.items():
    s = df[flux].loc[config['start']:].copy()
    s = s[s > 0]
    series_dict[flux] = s

print("Flux disponibles :")
for flux, s in series_dict.items():
    print(f"  {flux}: {len(s)} observations, de {s.index[0].date()} a {s.index[-1].date()}")

What does this code do?

  • La boucle filtre chaque flux a partir de sa date de debut configuree dans flux_config.
  • s = s[s > 0] supprime les semaines a zero poses qui correspondent a des periodes d'inactivite. Garder ces zeros fausserait les modeles de prevision.
  • series_dict stocke les 3 series nettoyees dans un dictionnaire pour iterer facilement dessus dans les etapes suivantes.

Avant de lire la sortie, a ton avis, combien d'observations devrait avoir le flux Resorption par rapport aux deux autres ?

Tu devrais voir un affichage listant les 3 flux avec leur nombre d'observations et la plage de dates. Resorption a nettement moins d'observations car il demarre 18 mois plus tard.

KeyError on a column name?

  • Verifie que ton DataFrame df contient bien les colonnes Maintenance, Croissance et Resorption avec ces noms exacts (majuscule en premiere lettre).
  • Lance df.columns.tolist() pour voir les noms reels de tes colonnes.
  • help me fix a KeyError when accessing columns in my Linky DataFrame

✔️ Awesome, I've got everything!

Parfait. Verifie que tu as bien sauvegarde ton notebook avant de passer a l'etape suivante.

⨯ I'd like to double check the full code

Voici les sections que tu devrais avoir dans ton notebook a ce stade. Confirme que chaque cellule est presente dans cet ordre :

  • Cellule 1 : le bloc d'imports complet (de import warnings jusqu'a import datetime).
  • Cellule 2 : le chargement CSV avec pd.read_csv et df.index.freq = 'W-MON'.
  • Cellule 3 : les parametres globaux (last_forecast, first_date_reso, first_date) et flux_config.
  • Cellule 4 : la boucle de filtrage qui cree series_dict et affiche le recapitulatif des 3 flux.

Ton environnement est pret et tes donnees sont chargees. Dans l'etape suivante, tu vas plonger dans l'analyse exploratoire pour comprendre la structure de chaque flux avant de choisir un modele.

Exploratory Data Analysis

Your environment is set up and your three Linky flux series are loaded. Before fitting any forecasting model, you need to understand what structure lives inside the data.

Time series have three key components to uncover: trend, seasonality, and stationarity. Without this diagnostic step, you risk fitting a model that ignores real patterns or hallucinates patterns that do not exist.

In this step, get ready to:
  • Visualize the raw series and decompose their seasonal patterns.
  • Run statistical tests (ADF, KPSS) to determine stationarity for each flux.
  • Examine ACF and PACF plots to identify autocorrelation structures.
Visualize the series and decompose seasonality

The first thing to do with any time series is simply look at it. A raw plot reveals obvious trends, seasonal oscillations, and anomalies that no test can replace.

  • In forecasting_linky_stats.py, add the following code below the existing Step 1 code to plot all three flux on a shared time axis:
# =============================================================================
# STEP 2 : ANALYSE EXPLORATOIRE (EDA)
# =============================================================================

# --- 2.1 Visualisation et decomposition saisonniere ---
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)
for i, (flux, s) in enumerate(series_dict.items()):
    axes[i].plot(s.index, s.values, linewidth=0.8)
    axes[i].set_title(f'Flux {flux}')
    axes[i].set_ylabel('Poses/semaine')
plt.tight_layout()
plt.show()

What does this code do?

It creates a 3-row subplot figure where each row displays one flux (Maintenance, Croissance, Resorption) over time.

The sharex=True parameter aligns all three time axes so you can compare seasonal peaks and troughs across flux at the same dates.

  • Run the cell. You should see three line charts stacked vertically, each showing weekly pose volumes over time.

Now you need to decompose each series into its structural components. Seasonal decomposition separates a time series into trend, seasonal pattern, and residual noise.

  • Add the following code below the previous block to run an additive seasonal decomposition with a 52-week period:
for flux, s in series_dict.items():
    if len(s) >= 104:
        decomp = seasonal_decompose(s, model='additive', period=52)
        fig = decomp.plot()
        fig.suptitle(f'Decomposition saisonniere - {flux}', y=1.02)
        plt.tight_layout()
        plt.show()

What does this code do?

  • seasonal_decompose splits each series into observed, trend, seasonal, and residual components.
  • period=52 tells statsmodels that one full seasonal cycle is 52 weeks (one year of weekly data).
  • The if len(s) >= 104 guard ensures at least 2 full seasonal cycles exist before attempting decomposition. Resorption started later, so this protects against insufficient data.
  • Run the cell. For each flux, you should see a 4-panel plot showing the observed series, extracted trend, repeating seasonal pattern, and leftover residuals.

Seeing a period error?

  • Make sure your DataFrame index has frequency set to W-MON. If you loaded from CSV without setting df.index.freq = 'W-MON', the decomposition cannot infer the period.
  • Check that series_dict contains data. Re-run the Step 1 cells if the dictionary is empty.
  • Why is my seasonal decompose failing?
Test stationarity with ADF and KPSS

A stationary series has constant mean and variance over time. Most forecasting models (ARIMA family) require stationarity, so you need to know whether differencing is needed.

Two complementary tests confirm this. The ADF test checks if a unit root exists (H0 = non-stationary). The KPSS test checks the opposite (H0 = stationary). Running both avoids ambiguous conclusions.

  • Add the following function definition below your decomposition code to create a reusable stationarity test:
# --- 2.2 Tests de stationnarite ---
def test_stationarity(series, name):
    adf_stat, adf_p, _, _, _, _ = adfuller(series.dropna(), autolag='AIC')
    kpss_stat, kpss_p, _, _ = kpss(series.dropna(), regression='c', nlags='auto')
    return {
        'Flux': name,
        'ADF Stat': round(adf_stat, 4),
        'ADF p-value': round(adf_p, 4),
        'ADF Conclusion': 'Stationnaire' if adf_p < 0.05 else 'Non stationnaire',
        'KPSS Stat': round(kpss_stat, 4),
        'KPSS p-value': round(kpss_p, 4),
        'KPSS Conclusion': 'Stationnaire' if kpss_p > 0.05 else 'Non stationnaire',
    }

What does this function do?

  • adfuller runs the Augmented Dickey-Fuller test. A p-value below 0.05 rejects the null hypothesis of a unit root, meaning the series IS stationary.
  • kpss runs the KPSS test with the opposite null. A p-value above 0.05 means you cannot reject stationarity.
  • When both tests agree, the conclusion is clear. When they disagree, the series may be trend-stationary (needs differencing).
  • Now add the code that runs this function on each flux and prints a summary table:
stationarity_results = []
for flux, s in series_dict.items():
    stationarity_results.append(test_stationarity(s, flux))

df_stationarity = pd.DataFrame(stationarity_results)
print("\n=== Tests de stationnarite ===")
print(df_stationarity.to_string(index=False))

What does this code do?

It iterates over all three flux, runs both stationarity tests on each, then assembles the results into a clean DataFrame with automatic conclusions.

  • Run both cells. You should see a printed table with columns for each test's statistic, p-value, and a conclusion label for Maintenance, Croissance, and Resorption.

Getting an interpolation warning from KPSS?

  • This is normal. The KPSS test sometimes reports that the p-value is outside its interpolation range. The warnings.filterwarnings('ignore') at the top of your file suppresses these. If you still see warnings, confirm that line is present and was executed.
  • My stationarity tests are returning unexpected results.
Plot ACF and PACF for each flux

The ACF (Autocorrelation Function) shows how correlated a series is with its own past values at different lags. The PACF (Partial Autocorrelation Function) isolates the direct effect of each lag after removing intermediate correlations.

These plots guide the choice of AR and MA orders in SARIMA models. Seasonal spikes at lag 52 confirm the annual cycle you saw in the decomposition.

  • Before you run this next block, predict: do you expect to see significant spikes at lag 52 in the ACF, given what the seasonal decomposition showed you?
  • Add the following code to plot ACF and PACF with 52 lags for each flux:
# --- 2.3 ACF et PACF ---
for flux, s in series_dict.items():
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
    plot_acf(s.dropna(), lags=52, ax=ax1, title=f'ACF - {flux}')
    plot_pacf(s.dropna(), lags=52, ax=ax2, title=f'PACF - {flux}')
    plt.tight_layout()
    plt.show()

What does this code do?

  • plot_acf draws the autocorrelation at each lag from 0 to 52. Bars extending beyond the blue confidence band are statistically significant.
  • plot_pacf draws the partial autocorrelation. A sharp cutoff in the PACF suggests an AR process. A slow decay suggests MA terms are needed.
  • Using 52 lags lets you see the full annual cycle and spot seasonal autocorrelation spikes.
  • Run the cell. You should see paired ACF/PACF plots for each flux. Look for slow decay in the ACF (confirming non-stationarity or strong trend) and spikes near lag 52 (confirming annual seasonality).

Plots look empty or flat?

  • Verify that series_dict is not empty. Run print(len(series_dict['Maintenance'])) to confirm data is present.
  • If the PACF plot raises a LinAlgError, your series may contain exact zeros or NaN values. The s.dropna() call should handle NaN. Check that the Step 1 filter s = s[s > 0] removed zeros.
  • My ACF/PACF plots are not displaying correctly.

✔️ Awesome, I've got everything!

Confirm you have saved forecasting_linky_stats.py and that all plots rendered correctly.

ⓧ I'd like to double check the full code

Your file should contain all the code from Step 1 (imports, parameters, data loading, series_dict), followed by the Step 2 EDA section. Scroll through and confirm these chunks appear in order:

  • The section header # STEP 2 : ANALYSE EXPLORATOIRE (EDA) followed by the 3-subplot raw series visualization with plt.show().
  • The seasonal_decompose loop with the if len(s) >= 104 guard and period=52.
  • The test_stationarity() function definition with ADF and KPSS logic.
  • The loop calling test_stationarity() on each flux, building df_stationarity, and printing the results.
  • The ACF/PACF loop using plot_acf and plot_pacf with lags=52.

You now have a complete diagnostic picture of your Linky data. The decomposition revealed the seasonal shape, the ADF/KPSS tests quantified stationarity, and the ACF/PACF plots exposed the autocorrelation structure. Next, you will use these insights to fit your first forecasting model: ETS (Exponential Smoothing).

Build the ETS Model

L'analyse exploratoire a revele la structure de chaque flux : tendance, saisonnalite annuelle, niveau de stationnarite. Il est temps de passer a la modelisation en commencant par le modele ETS (Error-Trend-Seasonal).

ETS est le modele de lissage exponentiel le plus complet. On va tester systematiquement toutes les combinaisons de composantes (error, trend, seasonal) pour trouver la meilleure configuration par flux. Avec seasonal_periods=52, certaines configurations vont echouer ou produire des residus imparfaits.

In this step, get ready to:
  • Definir une fonction de grid search ETS testant toutes les combinaisons valides.
  • Evaluer les residus du meilleur modele avec des tests statistiques (Ljung-Box, Jarque-Bera, ARCH).
  • Observer les limites du modele ETS sur les residus.
Definir la fonction de grid search ETS

Le grid search ETS teste systematiquement chaque combinaison de composantes error, trend et seasonal. Contrairement a auto_arima pour SARIMA, il n'existe pas de recherche automatique integree pour ETS.

  • Dans votre notebook, creez une nouvelle cellule sous la section EDA et ajoutez la premiere partie de la fonction grid_search_ets :
def grid_search_ets(series, seasonal_periods=52):
    """Grid search sur les composantes ETS."""
    results = []
    # Definir l'espace de recherche
    error_options = ['add', 'mul']
    trend_options = [None, 'add']
    seasonal_options = ['add', 'mul']
    damped_options = [True, False]

    for error in error_options:
        for trend in trend_options:
            for seasonal in seasonal_options:
                for damped in damped_options:
                    # Damped n'a de sens qu'avec un trend
                    if damped and trend is None:
                        continue
                    # Multiplicatif requiert des valeurs positives
                    if (error == 'mul' or seasonal == 'mul') and (series <= 0).any():
                        continue
  • Directement en dessous dans la meme cellule, ajoutez la logique de fitting et le retour du DataFrame :
                    try:
                        model = ETSModel(
                            series,
                            error=error,
                            trend=trend,
                            seasonal=seasonal,
                            seasonal_periods=seasonal_periods,
                            damped_trend=damped,
                        )
                        fit = model.fit(disp=False, maxiter=1000)
                        results.append({
                            'error': error,
                            'trend': trend,
                            'seasonal': seasonal,
                            'damped': damped,
                            'AIC': fit.aic,
                            'BIC': fit.bic,
                        })
                    except Exception:
                        continue

    return pd.DataFrame(results).sort_values('AIC').reset_index(drop=True)

Que fait cette fonction ?

  • Elle itere sur toutes les combinaisons possibles de composantes ETS : error (additif/multiplicatif), trend (aucun/additif), seasonal (additif/multiplicatif), damped (oui/non).
  • Les combinaisons invalides sont filtrees avant le fit : pas de damped sans trend, pas de multiplicatif sur des valeurs negatives ou nulles.
  • Chaque combinaison qui converge stocke son AIC et BIC. Les erreurs de convergence sont capturees silencieusement avec except Exception.
  • Le resultat est un DataFrame trie par AIC croissant. Le meilleur modele est en premiere ligne.

Pourquoi certaines combinaisons echouent ?

Avec seasonal_periods=52, le modele doit estimer 51 etats saisonniers initiaux. Cela rend l'optimisation difficile et certaines configurations ne convergent tout simplement pas. Le try/except permet de continuer la recherche sans interrompre le processus.

Maintenant, il faut une fonction pour evaluer la qualite des residus du meilleur modele.

  • Dans une nouvelle cellule en dessous, ajoutez la fonction evaluate_residuals :
def evaluate_residuals(residuals, model_name, flux_name):
    """Evaluation complete des residus d'un modele."""
    results = {}
    # Test Ljung-Box : autocorrelation des residus
    lb = acorr_ljungbox(residuals.dropna(), lags=[10, 20], return_df=True)
    results['Ljung-Box p (lag 10)'] = lb['lb_pvalue'].iloc[0]
    results['Ljung-Box p (lag 20)'] = lb['lb_pvalue'].iloc[1]
    # Test Jarque-Bera : normalite des residus
    jb_stat, jb_p, skew, kurtosis = jarque_bera(residuals.dropna())
    results['Jarque-Bera p'] = jb_p
    results['Skewness'] = skew
    results['Kurtosis'] = kurtosis
    # Test ARCH : heteroscedasticite
    arch_lm, arch_p, _, _ = het_arch(residuals.dropna(), nlags=5)
    results['ARCH LM p'] = arch_p
    results['Modele'] = model_name
    results['Flux'] = flux_name
    return results

Que testent ces 3 diagnostics ?

  • Ljung-Box : teste si les residus sont autocorreles. Si p < 0.05, le modele n'a pas capture toute la structure temporelle.
  • Jarque-Bera : teste la normalite des residus. Si p < 0.05, les residus ne suivent pas une loi normale.
  • ARCH (het_arch) : teste l'heteroscedasticite (variance non constante). Si p < 0.05, la variance des residus change dans le temps.
  • Executez ces deux cellules pour definir les fonctions (pas de sortie visible, c'est normal).
Lancer le grid search et evaluer les residus

Les fonctions sont pretes. On va maintenant lancer le grid search sur chaque flux, selectionner le meilleur modele par AIC, et evaluer ses residus.

  • Dans une nouvelle cellule, lancez le grid search ETS en ajoutant ce code :
ets_results = {}
ets_best_models = {}

for flux, s in series_dict.items():
    print(f"\n=== Grid Search ETS - {flux} ===")
    df_ets = grid_search_ets(s, seasonal_periods=52)

    if len(df_ets) > 0:
        print(df_ets.head(5).to_string(index=False))
        best = df_ets.iloc[0]
        # Re-fit le meilleur modele pour stocker l'objet complet
        best_model = ETSModel(
            s,
            error=best['error'],
            trend=best['trend'],
            seasonal=best['seasonal'],
            seasonal_periods=52,
            damped_trend=best['damped'],
        )
        best_fit = best_model.fit(disp=False, maxiter=1000)
        ets_best_models[flux] = best_fit
        ets_results[flux] = evaluate_residuals(best_fit.resid, 'ETS', flux)
        print(f"\nMeilleur ETS: error={best['error']}, trend={best['trend']}, "
              f"seasonal={best['seasonal']}, damped={best['damped']}")
        print(f"AIC: {best['AIC']:.2f}")
    else:
        print("Aucun modele ETS n'a converge.")

Que fait cette boucle ?

  • Pour chaque flux (Maintenance, Croissance, Resorption), elle lance le grid search et affiche les 5 meilleures combinaisons par AIC.
  • Le meilleur modele est re-fit et stocke dans ets_best_models pour pouvoir acceder a ses residus ensuite.
  • Les resultats des tests de residus sont stockes dans ets_results pour la comparaison finale entre modeles.
  • Executez la cellule. Le grid search peut prendre 1 a 2 minutes par flux.

Vous devriez voir un tableau des 5 meilleures combinaisons pour chaque flux, avec le meilleur modele identifie.

Le grid search prend trop longtemps ?

Avec m=52, certaines combinaisons mettent du temps a converger. Si la cellule tourne depuis plus de 5 minutes, interrompez le kernel et relancez. Le maxiter=1000 limite deja le nombre d'iterations.

help me debug my ETS grid search that takes too long

Le grid search identifie le meilleur modele. Mais un bon AIC ne garantit pas des residus propres. Generons le plot diagnostique pour inspecter visuellement les residus.

  • Dans une nouvelle cellule, ajoutez le code du plot diagnostique :
first_flux = list(ets_best_models.keys())[0]
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
residuals = ets_best_models[first_flux].resid

# Plot des residus dans le temps
axes[0, 0].plot(residuals)
axes[0, 0].set_title('Residus')
axes[0, 0].axhline(y=0, color='r', linestyle='--')

# Distribution des residus
axes[0, 1].hist(residuals.dropna(), bins=30, density=True, alpha=0.7)
axes[0, 1].set_title('Distribution des residus')

# QQ-Plot pour la normalite
from scipy import stats
stats.probplot(residuals.dropna(), plot=axes[1, 0])
axes[1, 0].set_title('QQ-Plot')

# Correlogramme des residus
plot_acf(residuals.dropna(), lags=40, ax=axes[1, 1], title='ACF Residus')
plt.suptitle(f'Diagnostique ETS - {first_flux}')
plt.tight_layout()
plt.show()

Avant d'executer cette cellule, faites une prediction : pensez-vous que les residus seront propres (pas d'autocorrelation) ou qu'il restera des patterns non captures ?

  • Executez la cellule pour generer le plot diagnostique.

Observez le graphique ACF des residus (en bas a droite). Plusieurs barres depassent l'intervalle de confiance. Le test Ljung-Box confirme : la p-value est inferieure a 0.05, ce qui signifie que les residus sont autocorreles.

Pourquoi les residus ETS montrent de l'autocorrelation ?

ETS capture la tendance et la saisonnalite, mais il ne modelise pas les structures AR/MA (autoregressive / moving average) presentes dans les donnees. Les pics dans l'ACF des residus revelent des dependances temporelles que le lissage exponentiel ne peut pas capturer.

C'est exactement ce qui motive le passage au modele SARIMA dans l'etape suivante. SARIMA ajoute les composantes AR et MA pour capturer ces structures residuelles.

Pas de plot ou erreur KeyError ?

Verifiez que la cellule du grid search a bien termine son execution et qu'au moins un modele a converge (ets_best_models ne doit pas etre vide). Si le flux Resorption n'a pas converge (moins de 2 ans de donnees), le premier flux sera Maintenance.

help me fix my ETS diagnostic plot

✔️ Awesome, I've got everything!

Verifiez que vous avez bien sauvegarde votre notebook.

ⓧ I'd like to double check the full code

Remontez dans votre notebook et confirmez que vous avez, dans l'ordre :

  • Cellule 1 : la fonction grid_search_ets complete (definition + boucle + try/except + return).
  • Cellule 2 : la fonction evaluate_residuals complete (Ljung-Box + Jarque-Bera + het_arch).
  • Cellule 3 : la boucle d'execution (ets_results et ets_best_models remplis).
  • Cellule 4 : le plot diagnostique 2x2 (residus, histogramme, QQ-plot, ACF).

Le modele ETS a trouve la meilleure combinaison de composantes, mais ses residus trahissent des dependances temporelles non capturees. Dans l'etape suivante, SARIMA va ajouter les composantes AR et MA pour corriger cette autocorrelation.

Build the SARIMA Model

Your ETS grid search revealed a clear problem. The Ljung-Box test flagged significant autocorrelation in the residuals, meaning ETS left predictable patterns in the errors it could not capture.

SARIMA (Seasonal ARIMA) captures exactly the AR/MA structures that ETS misses. Instead of manually testing hundreds of (p,d,q)(P,D,Q,52) combinations, auto_arima from pmdarima automates the search using a stepwise algorithm optimised for AIC.

In this step, get ready to:
  • Run auto_arima with stepwise search on each flux.
  • Evaluate SARIMA residuals and compare with ETS.
Run auto_arima with stepwise search

With 52-week seasonality, an exhaustive grid search over all ARIMA orders would take hours. The stepwise algorithm narrows the search intelligently, and StepwiseContext(max_dur=120) caps the total search time at 120 seconds per flux to prevent runaway computations.

  • In your notebook, create a new cell below the ETS diagnostic plot section and paste the following code:
sarima_best_models = {}
sarima_results = {}

for flux, s in series_dict.items():
    print(f"\n=== auto_arima SARIMA - {flux} ===")

    with StepwiseContext(max_dur=120):
        model = pm.auto_arima(
            s,
            start_p=0, start_q=0,
            max_p=3, max_q=3,
            start_P=0, start_Q=0,
            max_P=2, max_Q=2,
            m=52,
            seasonal=True,
            d=None, D=None,
            information_criterion='aic',
            trace=True,
            error_action='ignore',
            suppress_warnings=True,
            stepwise=True,
            n_jobs=-1,
        )

    sarima_best_models[flux] = model
    print(f"\nMeilleur modele: {model.order} x {model.seasonal_order}")
    print(f"AIC: {model.aic():.2f}")

    resid = pd.Series(model.resid(), index=s.index)
    sarima_results[flux] = evaluate_residuals(resid, 'SARIMA', flux)

What does this code do?

  • StepwiseContext(max_dur=120) limits the search to 120 seconds per flux. Without it, the algorithm could run indefinitely on highly seasonal data with m=52.
  • max_p=3, max_q=3 bounds the non-seasonal AR/MA orders, while max_P=2, max_Q=2 bounds the seasonal orders. These are reasonable upper limits for weekly data.
  • trace=True prints each candidate model's AIC as the algorithm searches, so you can follow the selection process live.
  • The loop reuses your evaluate_residuals() function from the ETS step to apply the same Ljung-Box, Jarque-Bera, and ARCH tests on the SARIMA residuals.
  • Run the cell.

This cell takes time

With m=52, each candidate SARIMA model is expensive to fit. Expect 1-3 minutes per flux depending on your machine. The trace=True output confirms the search is progressing.

You should see a stream of candidate models with their AIC values, followed by the best order printed for each flux (e.g. Meilleur modele: (1, 1, 1) x (1, 0, 1, 52)).

Search seems stuck or timed out?

  • The StepwiseContext caps at 120 seconds. If you see a timeout message, the best model found within that window is still stored.
  • If the kernel crashes (memory issue), restart it and re-run all cells from the top. Consider closing other applications to free RAM.

Help me debug my auto_arima run with m=52 seasonal data.

Evaluate SARIMA residuals and compare with ETS

The residual evaluation is already computed inside the loop. Now you need the visual diagnostic to confirm what the numbers tell you.

  • In a new cell below, add the diagnostic plot for the first flux by pasting this code:
first_flux = list(sarima_best_models.keys())[0]
sarima_best_models[first_flux].plot_diagnostics(figsize=(12, 8))
plt.suptitle(f'Diagnostique SARIMA - {first_flux}')
plt.tight_layout()
plt.show()

What does this code do?

plot_diagnostics() produces a 2x2 grid: standardized residuals over time, a histogram with KDE, a QQ-plot for normality, and a correlogram (ACF of residuals). This is the same layout you built manually for ETS, but built into pmdarima.

Before you run this cell, think about what you expect to see. Do you think the SARIMA correlogram will show fewer significant lags than the ETS correlogram?

  • Run the cell.

The correlogram should show a clear improvement over ETS. Most lags fall within the confidence bands. However, you may notice some persistent spikes that align with holiday or vacation periods.

Why do some spikes persist?

SARIMA models the regular seasonal cycle (52-week pattern) but cannot account for irregular calendar events like public holidays or school vacation periods. These events shift pose volumes in ways that fall outside the fixed seasonal structure.

This is exactly what motivates the next step: adding calendar effects as exogenous variables in a SARIMAX model.

  • Compare the AIC values printed by the SARIMA loop against the ETS AIC values from the previous step. SARIMA should show a lower (better) AIC for each flux.
  • Check the Ljung-Box p-values in sarima_results. At the first lags, they should be above 0.05, indicating the autocorrelation issue from ETS is resolved.

Ljung-Box p-values still below 0.05?

  • Check which lags are significant. If only higher lags (20+) are borderline, the model is capturing most structure. The remaining signal is likely calendar-driven.
  • If early lags (1-5) are still significant, the stepwise search may have been constrained by time. Try increasing max_dur to 180 or 240.

Help me interpret my SARIMA residual diagnostics.

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with both the auto_arima loop and the diagnostic plot cells.

ⓧ I'd like to double check the full code

Your SARIMA section should contain these two blocks in order, right after the ETS diagnostic plot section:

  • The auto_arima loop cell: initialises sarima_best_models and sarima_results dicts, iterates over series_dict, calls pm.auto_arima() inside a StepwiseContext, stores the model, prints order and AIC, then evaluates residuals.
  • The diagnostic plot cell: picks the first flux from sarima_best_models, calls plot_diagnostics(figsize=(12, 8)), adds a suptitle, and shows the figure.

Scroll up and confirm both cells are present and match. Make sure n_jobs=-1 is set in the auto_arima call and max_dur=120 is set in the StepwiseContext.

SARIMA captures the regular seasonal AR/MA structures that ETS missed. But those persistent spikes around holidays and vacations tell you the model still has a blind spot. In the next step, you will feed calendar effects directly into the model as exogenous variables with SARIMAX.

Build the SARIMAX Model with Calendar Features

Your SARIMA model improved the residual autocorrelation that ETS left behind. But if you look closely at the diagnostic plots, patterns still cluster around specific periods of the year.

Those residual spikes often coincide with French public holidays and school vacations, when field teams deploy fewer concentrators. SARIMAX lets you feed these calendar effects as exogenous variables so the model can capture what pure ARIMA structure cannot.

In this step, get ready to:
  • Build per-week calendar features from French holidays and school vacations.
  • Fit a SARIMAX model on each flux with exogenous variables and evaluate residuals.
  • Generate forecasts to 2028-12-31 with confidence intervals and compare all three models.
Build the calendar feature function

You need a function that takes a DatetimeIndex of weekly dates and returns a DataFrame with binary indicators for public holidays and school vacations in each week. The holidays package covers French public holidays. The vacances-scolaires-france package covers school vacation periods by zone.

  • In forecasting_linky_stats.py, add a new section header and the build_calendar_features function below the SARIMA diagnostic plot code:
# =============================================================================
# STEP 5 : MODELE SARIMAX AVEC EFFETS CALENDAIRES
# =============================================================================

def build_calendar_features(date_index):
    # Jours feries France
    fr_holidays = holidays.France(years=range(
        date_index.min().year, date_index.max().year + 1
    ))
    # Vacances scolaires (zone C = Paris/Ile-de-France)
    school_hol = SchoolHolidayDates()
    records = []
    for week_start in date_index:
        week_dates = pd.date_range(week_start, periods=7, freq='D')
        # Compter les jours feries dans la semaine
        n_feries = sum(1 for d in week_dates if d.date() in fr_holidays)
        # Verifier si la semaine tombe en vacances scolaires
        is_vacances = any(
            school_hol.is_holiday_for_zone(d.date(), 'C')
            for d in week_dates
        )
        records.append({
            'week': week_start,
            'n_jours_feries': n_feries,
            'is_jour_ferie': int(n_feries > 0),
            'is_vacances_scolaires': int(is_vacances),
        })
    exog = pd.DataFrame(records).set_index('week')
    return exog

What does this function do?

For each Monday in the date index, the function scans the 7 days of that week. It counts how many public holidays fall in the week (n_jours_feries), flags whether at least one exists (is_jour_ferie), and checks whether any day falls during a Zone C school vacation (is_vacances_scolaires).

Zone C covers Paris and Ile-de-France, chosen as a default reference zone. The learner can adjust this to Zone A or B for other regions.

Now call the function to build the exogenous DataFrame covering the full period from data start to the forecast horizon.

  • Add this code directly below the function definition:
# Construire exog pour toute la periode (historique + prevision)
full_index = pd.date_range(
    start=df.index.min(),
    end=last_forecast,
    freq='W-MON'
)
exog_full = build_calendar_features(full_index)
print(f"\nVariables exogenes construites: {exog_full.shape}")
print(exog_full.head(10))

Why cover the full period?

The exogenous variables must exist for both training and forecasting. Building them from the earliest data point to last_forecast = '2028-12-31' ensures one consistent DataFrame that you can slice for training and forecasting later.

  • Run the cell.
  • You should see a print statement showing the shape (several hundred rows by 3 columns) and the first 10 rows of the exogenous DataFrame with columns n_jours_feries, is_jour_ferie, and is_vacances_scolaires.

Seeing an error from SchoolHolidayDates?

If you get a date range error, check that vacances-scolaires-france is installed correctly. The package only covers dates within its internal calendar range. If your last_forecast extends beyond the package's data, the function may return partial results for far-future weeks.

help me debug my SchoolHolidayDates calendar feature construction

Fit SARIMAX and evaluate residuals

With the exogenous features ready, you can now fit a SARIMAX model on each flux. The key insight: you reuse the same ARIMA order that auto_arima found in Step 4, but add the calendar variables as exog.

  • Add the SARIMAX fitting loop below the exog_full print statements:
sarimax_best_models = {}
sarimax_results = {}

for flux, s in series_dict.items():
    print(f"\n=== SARIMAX - {flux} ===")
    # Aligner exog avec la serie d'entrainement
    exog_train = exog_full.loc[s.index]
    order = sarima_best_models[flux].order
    seasonal_order = sarima_best_models[flux].seasonal_order
    try:
        model = SARIMAX(
            s,
            exog=exog_train,
            order=order,
            seasonal_order=seasonal_order,
            enforce_stationarity=False,
            enforce_invertibility=False,
        )
        fit = model.fit(disp=False, maxiter=200)
        sarimax_best_models[flux] = fit
        sarimax_results[flux] = evaluate_residuals(fit.resid, 'SARIMAX', flux)
        print(f"Ordre: {order} x {seasonal_order}")
        print(f"AIC: {fit.aic:.2f}")
        # Afficher les coefficients exogenes avec significativite
        print("\nCoefficients exogenes:")
        for var in exog_train.columns:
            if var in fit.params.index:
                coef = fit.params[var]
                pval = fit.pvalues[var]
                sig = '***' if pval < 0.01 else '**' if pval < 0.05 else '*' if pval < 0.1 else ''
                print(f"  {var}: {coef:.4f} (p={pval:.4f}) {sig}")
    except Exception as e:
        print(f"Erreur SARIMAX pour {flux}: {e}")

What does this loop do?

  • For each flux, it slices exog_full to match the training dates and fits a SARIMAX with the same (p,d,q)(P,D,Q,52) order from SARIMA.
  • It evaluates residuals with the same Ljung-Box, Jarque-Bera, and ARCH tests as before for consistent comparison.
  • The significance stars (***) tell you which calendar features matter. A p < 0.05 means the variable has a statistically significant effect on that flux.
  • Run the cell.
  • For each flux you should see the ARIMA order, the AIC value, and the exogenous coefficients with their p-values and significance stars.

SARIMAX fails to converge?

If you see a convergence warning, the model still produces results but may not be optimal. Try increasing maxiter to 500. For the Resorption flux with fewer observations, convergence can be harder.

help me fix SARIMAX convergence issues with my Linky data

Now build the comparison tables that bring together results from all three models.

  • Add the comparison table code below the SARIMAX loop:
# Rassembler les diagnostics de residus des 3 modeles
all_results = []
for flux in series_dict.keys():
    if flux in ets_results:
        all_results.append(ets_results[flux])
    if flux in sarima_results:
        all_results.append(sarima_results[flux])
    if flux in sarimax_results:
        all_results.append(sarimax_results[flux])

df_comparison = pd.DataFrame(all_results)
cols_display = ['Flux', 'Modele', 'Ljung-Box p (lag 10)', 'Jarque-Bera p', 'ARCH LM p']
print("\n=== TABLEAU COMPARATIF DES RESIDUS ===")
print(df_comparison[cols_display].to_string(index=False))

# Comparaison AIC entre les 3 modeles
aic_comparison = []
for flux in series_dict.keys():
    row = {'Flux': flux}
    if flux in ets_best_models:
        row['AIC_ETS'] = ets_best_models[flux].aic
    if flux in sarima_best_models:
        row['AIC_SARIMA'] = sarima_best_models[flux].aic()
    if flux in sarimax_best_models:
        row['AIC_SARIMAX'] = sarimax_best_models[flux].aic
    aic_comparison.append(row)

df_aic = pd.DataFrame(aic_comparison)
print("\n=== COMPARAISON AIC ===")
print(df_aic.to_string(index=False))

What do these tables tell you?

The residual diagnostic table shows whether each model's residuals are well-behaved. A Ljung-Box p-value above 0.05 means no significant autocorrelation remains. A Jarque-Bera p-value above 0.05 suggests normality. An ARCH LM p-value above 0.05 means no heteroscedasticity.

The AIC table directly compares model fit. Lower AIC is better. For each flux, the model with the lowest AIC is your best statistical fit.

Before you run this cell, which model do you expect to have the lowest AIC for each flux?

  • Run the cell.
  • You should see two printed tables: one comparing Ljung-Box, Jarque-Bera, and ARCH p-values across all 9 model-flux combinations, and one showing AIC values per model per flux.

Missing flux in the tables?

If a flux does not appear in sarimax_best_models, the SARIMAX fit failed for that flux. Check the error message printed earlier. The comparison tables only include models that fitted successfully.

help me understand why my comparison table has missing models

Generate final forecasts with confidence intervals

With the best model identified per flux, you can now generate forecasts from the last training observation all the way to 2028-12-31. The SARIMAX model uses the future calendar features you already built in exog_full.

  • Add the forecast generation loop below the AIC comparison print:
# Previsions finales avec intervalles de confiance
for flux, s in series_dict.items():
    if flux not in sarimax_best_models:
        continue
    fit = sarimax_best_models[flux]
    forecast_start = s.index[-1] + pd.Timedelta(weeks=1)
    forecast_end = pd.Timestamp(last_forecast)
    n_forecast = len(pd.date_range(forecast_start, forecast_end, freq='W-MON'))
    forecast_index = pd.date_range(forecast_start, periods=n_forecast, freq='W-MON')
    # Extraire les variables exogenes pour la periode de prevision
    exog_forecast = exog_full.loc[exog_full.index.isin(forecast_index)]
    if len(exog_forecast) < n_forecast:
        exog_forecast = exog_full.loc[forecast_index[0]:forecast_index[-1]]
    forecast = fit.get_forecast(steps=len(exog_forecast), exog=exog_forecast)
    pred_mean = forecast.predicted_mean
    pred_ci = forecast.conf_int(alpha=0.05)
    # Visualisation historique + prevision
    fig, ax = plt.subplots(figsize=(14, 5))
    ax.plot(s.index, s.values, label='Historique', linewidth=0.8)
    ax.plot(pred_mean.index, pred_mean.values, label='Prevision SARIMAX', color='red')
    ax.fill_between(
        pred_ci.index,
        pred_ci.iloc[:, 0],
        pred_ci.iloc[:, 1],
        alpha=0.2, color='red', label='IC 95%'
    )
    ax.set_title(f"Prevision {flux} - SARIMAX jusqu'a {last_forecast}")
    ax.legend()
    ax.set_ylabel('Poses/semaine')
    plt.tight_layout()
    plt.show()

print("\n=== PREVISIONS GENEREES AVEC SUCCES ===")

What does this produce?

For each flux, get_forecast predicts weekly values from the end of training data to December 2028. The conf_int(alpha=0.05) method returns the 95% confidence interval bounds.

The plot shows historical data in blue and predictions in red, with a shaded band for uncertainty. The band widens further into the future, reflecting growing forecast uncertainty.

Before you run this, do you expect the confidence intervals to stay narrow or widen significantly by 2028?

  • Run the cell.
  • You should see one plot per flux showing the historical series in blue and the red forecast line extending to 2028, with a pink shaded 95% confidence interval that fans out over time. The final print statement confirms === PREVISIONS GENEREES AVEC SUCCES ===.

Forecast plot looks flat or explodes?

A flat forecast usually means the seasonal component dominates and the model converged to a simple average. An exploding forecast (values going to infinity) suggests the model is non-stationary. Check that enforce_stationarity=False and enforce_invertibility=False are set in your SARIMAX call.

help me interpret my SARIMAX forecast that looks unrealistic

✔️ Awesome, I've got everything!

Great. Double check you have saved forecasting_linky_stats.py and confirm all forecast plots rendered correctly.

ⓧ I'd like to double check the full code

Your file is long at this point. Scroll through and confirm these sections exist in order:

  • Section header: STEP 5 : MODELE SARIMAX AVEC EFFETS CALENDAIRES.
  • The build_calendar_features function definition.
  • The full_index and exog_full construction with print statements.
  • The sarimax_best_models fitting loop with coefficient display.
  • The df_comparison and df_aic comparison tables.
  • The forecast loop with get_forecast and plot generation.
  • Final print: === PREVISIONS GENEREES AVEC SUCCES ===.

You now have a complete statistical forecasting pipeline. Three models tested, residuals evaluated, and the best model identified per flux with forecasts extending to 2028. Next up: a Secret Mission to stress-test these results with temporal cross-validation.

Secret mission

Validation Croisee Temporelle (Expanding Window)

Les metriques calculees sur un seul jeu de donnees d'entrainement peuvent etre trompeuses. Implementez une validation croisee temporelle avec fenetre croissante pour obtenir une estimation plus fiable de la performance de chaque modele.

Clean Up Your Resources

Clean Up Your Resources

Ce projet tourne entierement en local. Il n'y a aucun cout recurrent ni ressource cloud a surveiller. Decidez si vous souhaitez garder vos fichiers, liberer la memoire, ou tout supprimer.

Resources you used:

  • Fichiers Python locaux (forecasting_linky_stats.py, cross_validation_temporal.py).
  • Notebook Jupyter avec les cellules executees.
  • Packages Python installes (pmdarima, holidays, vacances-scolaires-france).

Keep everything running

Aucune action necessaire. Choisissez cette option si vous comptez continuer a travailler sur le Projet 2 (Prophet, XGBoost) prochainement.

  • Vos fichiers restent sur votre machine, prets a etre reutilises.
  • Les packages installes restent disponibles pour vos futurs projets de time series.
  • Le kernel Jupyter consomme de la RAM tant qu'il est actif. Fermez-le si vous n'en avez plus besoin immediatement.

Pause - I'll come back to this later

Liberez la memoire en arretant le kernel Jupyter. Vos fichiers restent intacts sur le disque.

  • Dans Jupyter Notebook, cliquez sur Kernel puis Shutdown pour liberer la RAM.
  • Fermez l'onglet du navigateur et arretez le serveur Jupyter dans le terminal avec Ctrl+C.
  • Vos fichiers forecasting_linky_stats.py et cross_validation_temporal.py restent en place pour la prochaine session.

Delete - I don't want to use this again

Supprimez tous les fichiers du projet et desinstallez les packages ajoutes.

  • Arretez le kernel Jupyter et fermez le serveur comme decrit dans l'onglet Pause.
  • Supprimez les fichiers du projet en executant cette commande depuis le dossier parent:
rm forecasting_linky_stats.py cross_validation_temporal.py
  • Desinstallez les packages ajoutes pendant ce projet:
pip uninstall pmdarima holidays vacances-scolaires-france -y
  • Si vous avez un notebook .ipynb associe, supprimez-le egalement.

Nice Work!

Nice Work!

You just built a complete time series forecasting notebook that systematically compares three statistical models on real-world Linky deployment data. Your notebook identifies the best model for each flux with rigorous statistical evidence.

What you learned:

  • Conducted a full exploratory analysis with seasonal decomposition, ADF/KPSS stationarity tests, and ACF/PACF correlograms to understand the structure of 3 Linky flux series before modelling.
  • Built a systematic GridSearch ETS pipeline, observed its residual autocorrelation limits, then improved with SARIMA via auto_arima, and finally captured calendar effects with SARIMAX using French public holidays and school vacation variables.
  • Evaluated residuals rigorously with Ljung-Box, Jarque-Bera, and ARCH tests, then produced a comparative table identifying the best model per flux with forecasts through 2028.
  • Secret Mission: implemented expanding window temporal cross-validation to get a robust, fold-by-fold estimate of each model's predictive performance.

Ready to quiz yourself?