Weekly Forecast Pipeline in Python

Build a time series forecasting pipeline comparing SARIMA, ETS, and XGBoost.

Introduction

30 Second Summary

When you manage a fleet of equipment that needs replacing on a schedule, knowing how many replacements to expect each week is the difference between smooth operations and scrambling for resources. Forecasting weekly demand lets you plan staff, materials, and logistics ahead of time.

In this project, you will build a complete weekly forecasting pipeline in a Jupyter Notebook that uses real French weather data from Open-Meteo, official public holidays, and school calendar data to predict Linky meter replacements across four model families.

What You'll Build

You run your notebook and get side-by-side forecasts from four competing models on one chart, a performance table ranking them by three metrics, and a visual analysis proving which validation strategy holds up best.

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

A statistical diagnostics dashboard showing STL decomposition, ACF/PACF plots, and stationarity test results that prove your data's seasonal structure and non-stationarity.
Four trained forecasting models (SARIMA, SARIMAX, ETS, XGBoost) with overlaid predictions you can compare visually on a single chart.
A walk-forward cross-validation comparison table (MAE, RMSE, MAPE) with sensitivity analysis showing how the train/test split date affects each model's robustness.
Secret Mission: implement a Diebold-Mariano statistical test and a weighted ensemble that outperforms any individual model.

Are there any prerequisites?

You should be comfortable with Python (pandas, numpy, matplotlib). You also need Python 3.12+ already installed on your machine. No prior time series experience is required.

Before We Start

Pour construire un pipeline de forecasting fiable, tu as besoin d'un environnement Python isole. Sans ca, les conflits de dependances entre librairies cassent la reproductibilite des resultats.

Dans cette etape, tu vas creer un dossier projet, un environnement virtuel, et installer toutes les librairies necessaires. Tu prepareras aussi ton fichier notebook dans VS Code pour pouvoir executer du code cellule par cellule.

In this step, get ready to:

Creer un environnement virtuel et installer les dependances.
Creer le notebook Jupyter dans VS Code et verifier les imports.

Creer l'environnement virtuel et installer les dependances

Un environnement virtuel isole les librairies de ce projet du reste de ton systeme. Chaque projet a ses propres versions de packages, sans interferences.

Appuie sur Cmd+Space (macOS) ou la touche Windows (Windows) pour ouvrir la recherche.
Tape Terminal et appuie sur Entree pour l'ouvrir.

Deplace-toi sur ton Bureau en executant cette commande :

cd ~/Desktop

Cree le dossier du projet et entre dedans en executant ces commandes :

mkdir linky-forecast
cd linky-forecast

Confirme que tu es dans le bon dossier en executant :

pwd

Tu devrais voir un chemin qui finit par linky-forecast.

Cree l'environnement virtuel en executant cette commande :

macOS

python3 -m venv .venv

Windows

python -m venv .venv

Active l'environnement virtuel :

macOS

source .venv/bin/activate

Windows

.venv\Scripts\activate

Tu devrais voir (.venv) apparaitre au debut de ta ligne de commande. Cela confirme que l'environnement est actif.

Tu ne vois pas (.venv) dans le terminal ?

Verifie que tu es bien dans le dossier linky-forecast avec pwd. Le dossier .venv doit etre visible quand tu lances ls -a.

Sur Windows, si tu as une erreur de permissions, execute Set-ExecutionPolicy -Scope CurrentUser RemoteSigned dans PowerShell puis reessaie l'activation.

Aide-moi a activer mon environnement virtuel Python.

Installe toutes les librairies du projet en executant cette commande :

pip install numpy pandas matplotlib statsmodels pmdarima xgboost scikit-learn ipykernel

Que fait cette commande ?

statsmodels fournit SARIMAX, ETS, STL, les tests ADF et KPSS.
pmdarima automatise la recherche des ordres SARIMA via auto_arima.
xgboost entraine le modele gradient boosting pour la comparaison.
scikit-learn fournit les metriques et le walk-forward split.
ipykernel permet a VS Code d'utiliser ton environnement virtuel comme kernel Jupyter.

L'installation prend 1 a 2 minutes. Tu verras defiler les packages telecharges puis un message Successfully installed a la fin.

L'installation echoue ?

Verifie que ton environnement est actif (tu dois voir (.venv) dans le prompt). Si tu as une erreur sur xgboost, verifie ta version Python avec python --version. La version 3.12 ou superieure est requise.

Aide-moi a resoudre un probleme d'installation pip.

A ton avis, que va afficher la commande de verification si toutes les librairies sont bien installees ?

Verifie que tout est en place en executant cette commande :

python -c "import statsmodels; import pmdarima; import xgboost; print('OK')"

Le terminal affiche OK sans erreur. Les trois librairies principales sont fonctionnelles.

Creer le notebook Jupyter dans VS Code et verifier les imports

Tu vas travailler dans un fichier .py avec des marqueurs # %% qui delimitent les cellules. VS Code les reconnait comme un notebook interactif grace a l'extension Jupyter.

Appuie sur Cmd+Space (macOS) ou la touche Windows (Windows) pour ouvrir la recherche.
Tape Visual Studio Code et appuie sur Entree pour l'ouvrir.
Dans VS Code, clique sur l'icone Extensions dans la barre laterale gauche (ou appuie sur Cmd+Shift+X sur macOS, Ctrl+Shift+X sur Windows).
Tape Jupyter dans la barre de recherche des extensions.
Clique sur Install sur l'extension Jupyter de Microsoft.
Installe aussi l'extension Python de Microsoft si elle n'est pas deja presente.

Pourquoi ces extensions ?

L'extension Jupyter detecte les marqueurs # %% dans un fichier .py et ajoute des boutons Run Cell au-dessus de chaque cellule. Tu executes le code section par section, comme dans un notebook classique.

Ouvre le dossier linky-forecast dans VS Code en cliquant sur File dans la barre de menu.
Clique sur Open Folder.
Selectionne le dossier linky-forecast sur ton Bureau, puis clique sur Open.

Cree le fichier du pipeline en executant cette commande dans le terminal integre de VS Code (Ctrl+` pour l'ouvrir) :

macOS

touch linky_forecast_pipeline.py

Windows

New-Item linky_forecast_pipeline.py

Confirme que le fichier existe en executant :

ls

Tu devrais voir linky_forecast_pipeline.py dans la liste.

Clique sur linky_forecast_pipeline.py dans l'explorateur de fichiers de VS Code pour l'ouvrir.
Selectionne ton interpreteur Python en appuyant sur Cmd+Shift+P (macOS) ou Ctrl+Shift+P (Windows).
Tape Python: Select Interpreter et choisis l'interpreteur qui contient .venv dans son chemin.

Tape cette premiere cellule de test dans le fichier pour verifier que les imports fonctionnent :

# %% Test des imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import STL
from pmdarima import auto_arima
from xgboost import XGBRegressor
from sklearn.metrics import mean_absolute_error

print("Tous les imports OK")

Que fait ce code ?

Cette cellule importe les librairies principales du pipeline. Si un import echoue, tu verras immediatement une erreur ModuleNotFoundError avec le nom du module manquant.

Sauvegarde le fichier avec Cmd+S (macOS) ou Ctrl+S (Windows).
Clique sur Run Cell qui apparait au-dessus de la ligne # %% Test des imports.

VS Code ouvre un panneau Python Interactive sur le cote. Tu devrais voir Tous les imports OK s'afficher dans la sortie.

Run Cell n'apparait pas ?

Recharge la fenetre VS Code avec Cmd+Shift+P (macOS) ou Ctrl+Shift+P (Windows) puis tape Developer: Reload Window. Verifie que l'extension Jupyter est bien installee et active.

Si VS Code te demande de selectionner un kernel, choisis l'interpreteur .venv qui apparait dans la liste.

Aide-moi a faire fonctionner les cellules Jupyter dans VS Code.

✔️ Awesome, I've got everything!

Parfait. Ton environnement est pret et les cellules Jupyter fonctionnent. Supprime le contenu de test du fichier pour repartir d'une page blanche au prochain step.

ⓧ I'd like to double check the full code

Ton fichier linky_forecast_pipeline.py devrait contenir uniquement la cellule de test des imports. Verifie que tu as le marqueur # %% en premiere ligne, suivi des imports et du print. Si tout fonctionne, tu peux vider le fichier pour le prochain step.

Ton environnement est pret. Tu as Python isole, toutes les librairies installees, et un fichier notebook fonctionnel dans VS Code. Dans le prochain step, tu vas generer les donnees synthetiques qui miment le remplacement Linky et les visualiser pour la premiere fois.

Set Up the Notebook and Imports

Your environment is ready with all the libraries installed. Now you need a workspace where you can run code in small, testable chunks and see results inline.

A Jupyter notebook gives you exactly that. In this step, you'll create the notebook file and load every library the pipeline needs, confirming they all import cleanly before writing any analysis code.

In this step, get ready to:

Create the Jupyter notebook file in VS Code.
Write and run all the imports to confirm the libraries load without errors.

Create the notebook file

VS Code's Jupyter extension supports a "percent format" for notebooks. Instead of a .ipynb file, you write a standard .py file and mark cells with # %% comments. VS Code detects these markers and lets you run each cell interactively.

In VS Code, press Cmd+N (macOS) or Ctrl+N (Windows) to create a new file.
Save it immediately with Cmd+S (macOS) or Ctrl+S (Windows).
Name the file linky_forecast_pipeline.py and save it inside your project folder.

Why a .py file instead of .ipynb?

The percent format (.py with # %% markers) gives you the interactivity of a notebook with the simplicity of a plain Python file.

It's easier to version-control, works with any linter, and runs as a normal script if needed.

Write and run the imports

The pipeline uses multiple libraries for different stages: data manipulation, plotting, statistical tests, and machine learning models. Loading them all upfront ensures nothing is missing before you start the analysis.

Copy the following code into linky_forecast_pipeline.py:

# %% [markdown]
# # Pipeline de Forecast Hebdo - Remplacement Linky Post-Deploiement

# %% Imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests
import holidays
from statsmodels.tsa.seasonal import STL, seasonal_decompose
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
from pmdarima import auto_arima
from xgboost import XGBRegressor
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error
import warnings
warnings.filterwarnings("ignore")

What does this code do?

# %% [markdown] marks a documentation cell. The title appears as a heading when VS Code renders the notebook.
# %% Imports starts an executable code cell. Everything below it until the next # %% marker runs together.
requests handles HTTP calls to the Open-Meteo weather API and the French school calendar API.
holidays provides the official list of French public holidays for computing real working days per week.
statsmodels provides the statistical tests (ADF, KPSS), decomposition (STL), and models (SARIMAX, ETS).
pmdarima wraps auto_arima for automatic SARIMA order selection.
xgboost provides the gradient boosting model for non-linear feature-based forecasting.

Save the file with Cmd+S (macOS) or Ctrl+S (Windows).
Click Run Cell above the # %% Imports marker, or press Ctrl+Enter with your cursor inside the cell.

What should I see?

An Interactive Window opens on the right side of VS Code. The cell runs silently with no output (imports don't print anything). If you see no errors, all libraries loaded successfully.

🙋‍♀️ Seeing a ModuleNotFoundError?

Make sure the Python interpreter in VS Code is set to your virtual environment. Click the interpreter name in the bottom-right status bar and select the one inside your .venv folder.
If the error names a specific package, switch to your terminal and run pip install <package-name> with your virtual environment activated.

Why can't VS Code find my installed packages?

Before you run the next command, do you think it will print OK or raise an error?

Confirm the libraries are available outside the notebook too by running this in your terminal:

python -c "import statsmodels; import pmdarima; import xgboost; import holidays; import requests; print('OK')"

The terminal prints OK with no errors. Your notebook and all its dependencies are confirmed working.

Command prints an error instead of OK?

Check that your virtual environment is activated in the terminal. You should see (.venv) at the start of your prompt.
If you see No module named followed by a package name, reinstall that specific package with pip install and try again.

Help me debug why my Python imports are failing in the terminal.

✔️ Awesome, I've got everything!

Great. Double check you've saved linky_forecast_pipeline.py and move on to the next step.

ⓧ I'd like to double check the full code

Your linky_forecast_pipeline.py should contain exactly this:

# %% [markdown]
# # Pipeline de Forecast Hebdo - Remplacement Linky Post-Deploiement

# %% Imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests
import holidays
from statsmodels.tsa.seasonal import STL, seasonal_decompose
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
from pmdarima import auto_arima
from xgboost import XGBRegressor
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error
import warnings
warnings.filterwarnings("ignore")

Your notebook is live and every library is confirmed working. Next up, you'll fetch real French weather data from Open-Meteo and build the exogenous variables from official sources.

Generate Synthetic Data and Visualize

Fetch Real Data and Build the Series

Your environment is locked and loaded. Now you need real-world data to feed your forecasting pipeline.

Instead of synthetic placeholders, you will pull actual historical temperature data from the Open-Meteo API (free, no key required), compute real French working days using the holidays library, and fetch official school holiday periods from the government open data portal.

In this step, get ready to:

Fetch 4 years of real daily temperature data from Open-Meteo and compute a weekly France-wide average.
Compute real weekly working days and school holiday indicators from official French sources.
Build the target series, visualize the data, and split into train/test sets.

Fetch real temperature and exogenous data

Open-Meteo provides free ERA5 reanalysis data from 1940 onward. You will query 6 cities spread across France (Lille, Paris, Strasbourg, Lyon, Bordeaux, Marseille) and average their daily temperatures to approximate a national mean, then resample to weekly.

In your notebook linky_forecast_pipeline.py, create a new cell below your imports and paste the following code to fetch temperature data and build the exogenous variables:

# %% Fetch real exogenous data from French sources

# 6 cities covering France north-to-south, east-to-west
cities = {
    "Lille": (50.63, 3.06),
    "Paris": (48.86, 2.35),
    "Strasbourg": (48.58, 7.75),
    "Lyon": (45.76, 4.84),
    "Bordeaux": (44.84, -0.58),
    "Marseille": (43.30, 5.37),
}

# Fetch daily mean temperature from Open-Meteo Historical API
start_date = "2021-01-04"
end_date = "2024-12-29"
temp_frames = []

for city, (lat, lon) in cities.items():
    url = "https://archive-api.open-meteo.com/v1/archive"
    params = {
        "latitude": lat,
        "longitude": lon,
        "start_date": start_date,
        "end_date": end_date,
        "daily": "temperature_2m_mean",
        "timezone": "Europe/Paris",
    }
    resp = requests.get(url, params=params)
    data = resp.json()
    temp_frames.append(
        pd.Series(data["daily"]["temperature_2m_mean"],
                  index=pd.to_datetime(data["daily"]["time"]),
                  name=city)
    )
    print(f"  {city}: {len(data['daily']['time'])} days fetched")

# Average across cities and resample to weekly (Monday)
temp_daily = pd.concat(temp_frames, axis=1).mean(axis=1)
temp_weekly = temp_daily.resample("W-MON").mean()
print(f"\nTemperature: {len(temp_weekly)} weeks")

What does this code do?

The 6 cities span the country from Lille (north) to Marseille (south), giving a representative cross-section of French climate.
archive-api.open-meteo.com/v1/archive returns ERA5 reanalysis daily data. No API key needed.
temperature_2m_mean is the 24-hour average air temperature at 2 meters above ground, in Celsius.
Resampling with resample("W-MON").mean() converts daily values to weekly averages aligned to Mondays.

Still in the same cell, add the following code below to compute working days, school holidays, and build the final DataFrame:

# Compute real working days per week using French public holidays
fr_holidays = holidays.FR(years=range(2021, 2025))

def working_days_in_week(week_start):
    """Count weekdays that are not public holidays in a Mon-Sun week."""
    count = 0
    for i in range(7):
        day = week_start + pd.Timedelta(days=i)
        if day.weekday() < 5 and day not in fr_holidays:
            count += 1
    return count

jours_ouvres = pd.Series(
    [working_days_in_week(d) for d in temp_weekly.index],
    index=temp_weekly.index
)

# Fetch school holidays from data.education.gouv.fr
cal_url = ("https://data.education.gouv.fr/api/explore/v2.1/catalog/datasets/"
           "fr-en-calendrier-scolaire/records?limit=100"
           "&refine=location%3AM%C3%A9tropole")
cal_resp = requests.get(cal_url)
cal_records = cal_resp.json().get("results", [])

# Build a set of dates that fall inside any vacation period
vacation_dates = set()
for rec in cal_records:
    start = pd.to_datetime(rec["start_date"]).normalize()
    end = pd.to_datetime(rec["end_date"]).normalize()
    vacation_dates.update(pd.date_range(start, end))

# For each week: 1 if any day falls in school holidays, else 0
vacances = pd.Series(
    [int(any((d + pd.Timedelta(days=i)) in vacation_dates
             for i in range(7)))
     for d in temp_weekly.index],
    index=temp_weekly.index
)
print(f"Working days range: {jours_ouvres.min()}-{jours_ouvres.max()}")
print(f"Weeks with school holidays: {vacances.sum()} / {len(vacances)}")

What does this code do?

holidays.FR() generates the official list of French public holidays (11 per year including Easter Monday, Ascension, etc.).
working_days_in_week() counts weekdays (Mon-Fri) that are not public holidays. Most weeks have 5, some have 4.
The school calendar API returns vacation periods for metropolitan France. Each record has a start_date and end_date that define one holiday block.
vacances is binary: 1 if any day in the week overlaps with school holidays, 0 otherwise.

Add a final block below to build the synthetic target series and assemble the DataFrame:

# Build the synthetic replacement target (real exog, synthetic target)
np.random.seed(42)
n_weeks = len(temp_weekly)
dates = temp_weekly.index

# Trend: aging fleet drives more replacements over time
trend = np.linspace(200, 350, n_weeks)

# Seasonal: summer peaks (weeks 18-39)
week_of_year = dates.isocalendar().week.values.astype(int)
seasonal = np.where(
    (week_of_year >= 18) & (week_of_year <= 39),
    50 + 30 * np.sin(np.pi * (week_of_year - 18) / 21),
    -20 + 10 * np.sin(2 * np.pi * week_of_year / 52)
)

# Exogenous effect from real variables
exog_effect = (2.0 * temp_weekly.values
               + 15 * jours_ouvres.values
               - 30 * vacances.values)

# Add noise and enforce realistic floor
noise = np.random.normal(0, 15, n_weeks)
y_values = trend + seasonal + exog_effect + noise
y_values = np.maximum(y_values, 50)

# Assemble final DataFrame
df = pd.DataFrame({
    "remplacements": y_values,
    "temperature": temp_weekly.values,
    "jours_ouvres": jours_ouvres.values,
    "vacances": vacances.values,
}, index=dates)
df.index.freq = "W-MON"

print(f"Dataset: {len(df)} weeks")
print(df.head())

Why is the target still synthetic?

Real Linky replacement data is confidential (Enedis internal). The exogenous drivers (temperature, working days, holidays) are all real, which gives the model realistic correlations to learn from.

The synthetic target combines these real drivers with a linear trend and noise. Your SARIMAX model will find genuine weather-replacement correlations, not fabricated ones.

Run the cell by clicking the play button (or pressing Shift+Enter).

You should see output confirming 208 weeks of data from 2021-01-04 to 2024-12-30, followed by the first 5 rows of your DataFrame.

Seeing a connection or JSON error?

Make sure you have internet access. The Open-Meteo API and data.education.gouv.fr both require an active connection.
If you see KeyError: 'daily', check that the dates in start_date and end_date are in the correct YYYY-MM-DD format.
If the school calendar API returns an empty list, verify the URL has the correct refine=location%3AMétropole parameter. The dataset may paginate — 100 records covers most vacation periods.

Help me debug my data fetching cell.

Visualize the raw series and exogenous variables

Numbers alone do not confirm that your data pipeline worked correctly. A multi-panel chart lets you visually verify that the real temperature cycle, school holidays, and the resulting replacement pattern all align as expected.

Create a new cell below and paste the following visualization code:

# %% Visualisation de la serie brute
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)

axes[0].plot(df.index, df["remplacements"], color="steelblue", linewidth=0.8)
axes[0].set_title("Remplacements Linky hebdomadaires")
axes[0].set_ylabel("Nombre de remplacements")

axes[1].plot(df.index, df["temperature"], color="darkorange", linewidth=0.8)
axes[1].set_title("Temperature moyenne hebdomadaire")
axes[1].set_ylabel("Temperature (C)")

axes[2].bar(df.index, df["vacances"], color="gray", alpha=0.5, width=5)
axes[2].set_title("Periodes de vacances scolaires")
axes[2].set_ylabel("Vacances (0/1)")

plt.tight_layout()
plt.show()

What does this chart show?

Panel 1 (replacements): the target series. You should spot an upward drift AND recurring summer peaks.
Panel 2 (temperature): a smooth annual cycle that correlates with summer peaks in replacements.
Panel 3 (vacances): binary spikes during school holidays. These cause dips in replacements because fewer technicians are available.

Before you run this cell, what do you expect to see in the top panel? Think about how a trend plus seasonality would look over 4 years.

Run the cell.

You should see a three-panel chart. The top panel displays 4 years of weekly data with a clearly rising baseline and peaks recurring every summer (May-September). The middle panel shows a smooth temperature sine wave. The bottom panel shows sparse gray bars during holiday periods.

Chart not displaying?

If you see a blank output, make sure you are running the notebook with the Jupyter extension in VS Code. The interactive window must be active for plots to render inline.
Check that you ran the data generation cell above first. The visualization depends on the df DataFrame existing in memory.

Help me get my matplotlib chart to display in VS Code.

Split train and test sets

Time series data cannot be shuffled randomly. You must respect chronological order so the model never sees future data during training. An 80/20 temporal split gives 166 weeks for training and 42 weeks for evaluation.

Create a new cell below and paste this code to perform the temporal split:

# %% Split train/test
split_idx = int(len(df) * 0.8)
train = df.iloc[:split_idx]
test = df.iloc[split_idx:]

y_train = train["remplacements"]
y_test = test["remplacements"]
exog_train = train[["temperature", "jours_ouvres", "vacances"]]
exog_test = test[["temperature", "jours_ouvres", "vacances"]]

print(f"Train: {len(train)} semaines | Test: {len(test)} semaines")

What does this code do?

split_idx takes the first 80% of rows chronologically (166 weeks) as training data.
y_train and y_test isolate just the target column (number of replacements).
exog_train and exog_test hold the exogenous variables separately. SARIMAX will need these later.

Run the cell.

You should see: Train: 166 semaines | Test: 42 semaines.

Why 80/20 and not 70/30?

With only 208 data points, a 70/30 split would leave 62 test points but reduce training to just 145 weeks. That is barely 2.8 seasonal cycles for a model with m=52 seasonality. 80/20 gives 3.2 full cycles for training while retaining enough test data (42 weeks, nearly a full year) to evaluate forecast quality across seasons.

🙋‍♀️ Seeing different numbers?

Confirm that n_weeks is set to 208 in your data generation cell. A different value changes the split sizes.
Make sure you did not accidentally filter or drop rows from df between the generation and split cells.

Help me check why my train/test split sizes are wrong.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook file.

ⓧ I'd like to double check the full code

Your notebook should now contain these cells in order. Scroll up and confirm each is present:

Cell 1: The markdown title cell (# %% [markdown]).
Cell 2: Imports cell (# %% Imports) with all library imports including requests and holidays.
Cell 3: Data fetching cell (# %% Fetch real exogenous data) with Open-Meteo API calls, working days computation, school holidays fetch, and DataFrame assembly.
Cell 4: Visualization cell (# %% Visualisation de la serie brute) with the 3-panel figure.
Cell 5: Train/test split cell (# %% Split train/test) producing y_train, y_test, exog_train, exog_test.

Your data is generated, visualized, and split. Next up, you will run statistical tests to formally prove that the series is non-stationary and identify the seasonal structure that your models will need to handle.

Run Statistical Tests and Decomposition

Your synthetic data is ready and you can see the trend and seasonality in the plot. But if you tried to fit a SARIMA model directly on this raw series, it would diverge or produce auto-correlated residuals.

The reason is simple: the series is non-stationary, and you haven't identified the correct differencing order. In this step, you'll run formal statistical tests to prove this, then decompose the signal to understand its components before any modeling.

In this step, get ready to:

Run ADF and KPSS stationarity tests with cross-interpretation.
Decompose the series into trend, seasonality, and residual using STL.
Plot ACF and PACF to identify seasonal periodicity and candidate orders.

Run stationarity tests: ADF and KPSS

Before choosing a model, you need to answer a critical question: is this series stationary? Stationarity means the statistical properties (mean, variance) don't change over time. SARIMA requires stationarity after differencing, so you need to know how many times to difference.

You'll use two complementary tests. The ADF test's null hypothesis is "the series is non-stationary." The KPSS test flips it: its null hypothesis is "the series is stationary." Running both gives you a cross-check that avoids false conclusions.

In your notebook, create a new cell below the train/test split cell and add the ADF test by pasting this code:

# %% Test ADF
adf_result = adfuller(y_train, autolag="AIC")
print("=" * 50)
print("TEST ADF (Augmented Dickey-Fuller)")
print(f"  Statistique ADF : {adf_result[0]:.4f}")
print(f"  P-value         : {adf_result[1]:.4f}")
print(f"  Lags utilises   : {adf_result[2]}")
print(f"  Conclusion      : {'Stationnaire' if adf_result[1] < 0.05 else 'NON stationnaire'}")

What does this code do?

adfuller(y_train, autolag="AIC") runs the Augmented Dickey-Fuller test on the training data, automatically selecting the number of lags using AIC.
If the p-value is above 0.05, you cannot reject the null hypothesis of non-stationarity. The series needs differencing before SARIMA can work.

Before you run this cell, what do you expect the conclusion to be? Think about what you saw in the visualization: a clear upward trend means the mean is changing over time.

Run the cell.

You should see a p-value greater than 0.05 and the conclusion "NON stationnaire". The ADF test confirms your series has a unit root.

Now create another new cell directly below and add the KPSS test to cross-check:

# %% Test KPSS
kpss_result = kpss(y_train, regression="c", nlags="auto")
print("=" * 50)
print("TEST KPSS")
print(f"  Statistique KPSS : {kpss_result[0]:.4f}")
print(f"  P-value          : {kpss_result[1]:.4f}")
print(f"  Conclusion       : {'Stationnaire' if kpss_result[1] > 0.05 else 'NON stationnaire'}")

Why use two tests instead of one?

ADF and KPSS have opposite null hypotheses. When both agree (ADF says non-stationary AND KPSS says non-stationary), you have strong evidence. If they disagree, the series might be trend-stationary or borderline.

The regression="c" parameter tests stationarity around a constant (level stationarity). The nlags="auto" lets statsmodels pick the optimal lag length automatically.

Run the KPSS cell.

KPSS also concludes "NON stationnaire". Both tests agree: this series needs at least one order of differencing (d=1) before fitting a SARIMA model.

Seeing a warning about p-value interpolation?

KPSS may print a warning that the p-value is outside the table range. This is normal when the test statistic is extreme. The conclusion is still valid. The warnings.filterwarnings("ignore") line in your imports cell should suppress it, but if you removed that line, the warning is harmless.

Help me interpret my ADF and KPSS results

Decompose the series with STL

You've confirmed non-stationarity. Now you need to see exactly what's inside the signal. STL (Seasonal-Trend decomposition using LOESS) splits your series into three additive components: trend, seasonal pattern, and residual noise.

Create a new cell below the KPSS cell and add the STL decomposition by pasting this code:

# %% Decomposition STL
stl = STL(y_train, period=52, robust=True)
stl_result = stl.fit()

fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)
axes[0].plot(stl_result.observed, color="steelblue")
axes[0].set_title("Observed")
axes[1].plot(stl_result.trend, color="darkorange")
axes[1].set_title("Trend")
axes[2].plot(stl_result.seasonal, color="green")
axes[2].set_title("Seasonal")
axes[3].plot(stl_result.resid, color="red", alpha=0.7)
axes[3].set_title("Residual")
plt.tight_layout()
plt.show()

What does this decomposition reveal?

period=52 tells STL the seasonality repeats every 52 weeks (annual cycle). This matches the summer peaks you saw in the raw data.
robust=True makes the decomposition resistant to outliers by using iterative reweighting.
The trend panel confirms the upward drift (why ADF detected non-stationarity). The seasonal panel shows the repeating annual pattern. The residual panel should look like random noise if the decomposition captured everything.

Run the cell.

You should see a 4-panel plot. The top panel shows the original series, followed by a smooth upward trend, a repeating green seasonal wave, and red residuals that look like random noise.

Plot not showing or cells running out of order?

Make sure you ran all previous cells in order (imports, data generation, split). If y_train is not defined, re-run the split cell first. If the plot appears but is blank, check that plt.show() is at the end of the cell.

Help me fix my STL decomposition plot

Plot ACF and PACF to identify seasonal periodicity

The final diagnostic is the autocorrelation function (ACF) and partial autocorrelation function (PACF). These plots show how correlated the series is with its own past values at different lags. Significant spikes at lag 52 confirm the annual seasonality and tell you the seasonal period (m=52) for SARIMA.

Create a new cell below the STL cell and add the ACF/PACF plots by pasting this code:

# %% ACF et PACF
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(y_train, lags=60, ax=axes[0], title="ACF")
plot_pacf(y_train, lags=60, ax=axes[1], title="PACF")
plt.tight_layout()
plt.show()

How to read ACF and PACF plots?

ACF shows total correlation at each lag. A slow decay suggests non-stationarity (trend). Spikes at lag 52 confirm annual seasonality.
PACF shows correlation at lag k after removing the effect of intermediate lags. Significant early spikes help identify the AR order (p) for SARIMA.
lags=60 shows correlations up to 60 weeks, enough to reveal the annual (52-week) spike clearly.

Before running this cell, do you expect to see significant spikes at any particular lag number? Think about the seasonality period you set in the STL decomposition.

Run the cell.

You should see significant spikes at lag 52 in both plots, confirming the annual seasonal periodicity. The ACF also shows a slow decay pattern typical of non-stationary data.

Spikes not clearly visible at lag 52?

If the spikes at lag 52 look small, check that your data generation cell used np.random.seed(42). A different seed produces different noise levels. The spike should still be above the blue confidence band. If nothing appears at all, make sure y_train has 166 values by running len(y_train) in a separate cell.

Help me understand my ACF and PACF plots

💡 What do these diagnostics tell you for modeling?

You now have three confirmed facts: the series is non-stationary (needs d=1 differencing), it has annual seasonality (m=52, needs D=1 seasonal differencing), and the ACF/PACF suggest low AR and MA orders. In the next step, auto_arima will use these constraints to find the optimal SARIMA configuration automatically.

✔️ Awesome, I've got everything!

Great. Make sure you've saved your notebook file.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm you have these 4 new cells in order after the train/test split cell:

# %% Test ADF cell with adfuller(y_train, autolag="AIC") and print statements.
# %% Test KPSS cell with kpss(y_train, regression="c", nlags="auto") and print statements.
# %% Decomposition STL cell with STL(y_train, period=52, robust=True) and 4-panel plot.
# %% ACF et PACF cell with plot_acf and plot_pacf with lags=60.

Your diagnostic toolkit is complete. You've proven the series is non-stationary, confirmed the 52-week seasonal cycle, and identified the differencing orders. Next up, you'll use these insights to train SARIMA, SARIMAX, and ETS models that properly handle these characteristics.

Train SARIMA, SARIMAX, and ETS Models

Your diagnostics confirmed what you suspected: the series is non-stationary (both ADF and KPSS agree) and has strong annual seasonality at lag 52.

Now you need models that handle these properties natively. In this step, you'll train three statistical forecasting models that each tackle trend and seasonality in a different way, then generate predictions on your test set.

In this step, get ready to:

Use auto_arima to automatically identify the best SARIMA order.
Train SARIMAX with exogenous variables (temperature, working days, holidays).
Train ETS (Holt-Winters additive) as a complementary baseline.

Fit SARIMA via auto_arima

A SARIMA model extends ARIMA with seasonal parameters (P, D, Q, m). Instead of guessing those orders manually from the ACF/PACF plots, auto_arima tests combinations using AIC to find the best fit.

Create a new cell in your notebook below the ACF/PACF cell, and paste the following code:

# %% SARIMA via auto_arima
print("Recherche du meilleur SARIMA via auto_arima...")
sarima_model = auto_arima(
    y_train,
    seasonal=True,
    m=52,
    start_p=0, start_q=0,
    max_p=3, max_q=3,
    start_P=0, start_Q=0,
    max_P=2, max_Q=2,
    d=1, D=1,
    trace=True,
    error_action="ignore",
    suppress_warnings=True,
    stepwise=True,
    information_criterion="aic"
)
print(f"\nMeilleur modele SARIMA: {sarima_model.order} x {sarima_model.seasonal_order}")

# Forecast SARIMA
sarima_forecast = sarima_model.predict(n_periods=len(test))
sarima_forecast = pd.Series(sarima_forecast, index=test.index)

What does this code do?

auto_arima performs a stepwise search over ARIMA orders (p, d, q) and seasonal orders (P, D, Q) with m=52 (annual seasonality on weekly data). You fix d=1 and D=1 because the tests confirmed non-stationarity requires one regular and one seasonal differencing.

trace=True prints each candidate model and its AIC score so you can watch the search progress. The final .predict() generates forecasts for the full test period (42 weeks).

Run the cell.

You should see a list of ARIMA combinations with their AIC scores scrolling by, followed by a summary line like Meilleur modele SARIMA: (p, 1, q) x (P, 1, Q, 52). The exact orders depend on the stepwise search.

Cell runs for a very long time?

With m=52, the seasonal ARIMA search is computationally heavy. The stepwise algorithm keeps it manageable (typically 1-3 minutes). If it exceeds 5 minutes, check that stepwise=True is set. Without stepwise, the search space explodes.

Help me debug why auto_arima is running too long on my weekly data with m=52

Train SARIMAX with exogenous variables

SARIMA only uses past values of the series itself. SARIMAX adds external regressors (temperature, working days, holidays) that you know influence replacement volumes. It reuses the same (p,d,q)(P,D,Q,m) orders that auto_arima found.

Create a new cell below the SARIMA cell, and paste the following code:

# %% SARIMAX avec variables exogenes
print("\nEntrainement SARIMAX avec variables exogenes...")
sarimax_model = SARIMAX(
    y_train,
    exog=exog_train,
    order=sarima_model.order,
    seasonal_order=sarima_model.seasonal_order,
    enforce_stationarity=False,
    enforce_invertibility=False
)
sarimax_fit = sarimax_model.fit(disp=False)
sarimax_forecast = sarimax_fit.forecast(steps=len(test), exog=exog_test)
sarimax_forecast.index = test.index
print("SARIMAX entraine.")

What does this code do?

SARIMAX takes the same ARIMA structure but feeds exog_train (temperature, jours_ouvres, vacances) during training, and exog_test during forecasting. This lets the model learn how external factors shift replacement volumes above or below the seasonal baseline.

enforce_stationarity=False and enforce_invertibility=False prevent numerical errors when the optimizer explores parameter regions near the boundary of stationarity.

Run the cell.

You should see SARIMAX entraine. printed. This confirms the model fitted and generated 42 weeks of forecasts.

Seeing a convergence warning?

Convergence warnings are common with high-order seasonal models. The warnings.filterwarnings("ignore") line in your imports cell suppresses them. If you still see one, confirm that line is present and the imports cell was executed.

Help me understand why SARIMAX shows convergence warnings

Train ETS (Holt-Winters additive)

ETS (Error, Trend, Seasonality) decomposes the series into additive components and forecasts each one with exponential smoothing. It provides a fundamentally different approach from the ARIMA family, which makes it a valuable comparison baseline.

Create a new cell below the SARIMAX cell, and paste the following code:

# %% ETS (Holt-Winters)
print("\nEntrainement ETS (Holt-Winters additif)...")
ets_model = ExponentialSmoothing(
    y_train,
    trend="add",
    seasonal="add",
    seasonal_periods=52,
    initialization_method="estimated"
)
ets_fit = ets_model.fit(optimized=True)
ets_forecast = ets_fit.forecast(steps=len(test))
ets_forecast.index = test.index
print("ETS entraine.")

What does this code do?

ExponentialSmoothing with trend="add" and seasonal="add" implements the additive Holt-Winters method. It separates the series into level + trend + season and updates each with a smoothing parameter.

initialization_method="estimated" lets statsmodels choose starting values via maximum likelihood, which is more robust than heuristic initialization for high-frequency seasonal data.

Run the cell.

You should see ETS entraine. printed. All three statistical models are now fitted and their forecasts stored.

Seeing a ValueError about seasonal_periods?

If you see an error mentioning the seasonal period, confirm that seasonal_periods=52 matches what you used in the STL decomposition earlier. Also verify your training set has at least 2 full seasonal cycles (104+ observations). Your y_train should have 166 rows.

Help me fix the seasonal_periods error in ExponentialSmoothing

Before you move on, take a moment to predict: which of the three models do you think will track the test data most closely? SARIMA (autoregressive structure only), SARIMAX (same structure plus external variables), or ETS (smoothing-based decomposition)?

Why SARIMAX should outperform here

Your synthetic data was generated with an explicit exogenous effect (2.0 * temperature + 15 * jours_ouvres - 30 * vacances). SARIMAX is the only model here that receives those variables. SARIMA and ETS must approximate those effects through autoregressive memory or smoothing alone.

✔️ Awesome, I've got everything!

Great. Confirm you have saved linky_forecast_pipeline.py and all three cells ran without errors.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm these three new cells appear in order after your ACF/PACF cell:

# %% SARIMA via auto_arima cell with auto_arima(...) and sarima_model.predict().
# %% SARIMAX avec variables exogenes cell with SARIMAX(...) and sarimax_fit.forecast().
# %% ETS (Holt-Winters) cell with ExponentialSmoothing(...) and ets_fit.forecast().

Each cell should have printed its confirmation message when you ran it.

You now have three trained statistical models with forecasts stored in sarima_forecast, sarimax_forecast, and ets_forecast. Next up, you'll add XGBoost to the mix with feature engineering and pit all four models against each other with walk-forward cross-validation.

Train XGBoost and Compare All Models

Les modeles statistiques (SARIMA, SARIMAX, ETS) capturent bien la structure lineaire de la serie. Mais ils peuvent rater des interactions non-lineaires entre les features calendaires, les lags et les variables exogenes.

Dans ce step, tu entraines un XGBoost avec du feature engineering temporel, puis tu compares les 4 familles de modeles avec une validation walk-forward rigoureuse. Tu analyseras aussi l'impact du choix de la strategie de validation (expanding vs sliding window) sur les metriques.

In this step, get ready to:

Construire la fonction de feature engineering (lags, rolling stats, Fourier, calendaire) et entrainer XGBoost.
Implementer la walk-forward cross-validation avec expanding et sliding window pour les 4 modeles.
Comparer les modeles avec un tableau de metriques (MAE, RMSE, MAPE) et analyser la sensibilite a la strategie de validation.

Build the feature engineering function for XGBoost

Les modeles statistiques utilisent la structure temporelle directement. XGBoost a besoin qu'on lui fournisse explicitement des features derivees de la serie : lags, moyennes glissantes, termes de Fourier et variables calendaires.

Cree une nouvelle cellule dans ton notebook et definis la fonction de feature engineering en copiant ce code :

# %% XGBoost - Feature Engineering
def create_features(df, target_col="remplacements"):
    """Cree les features temporelles pour XGBoost."""
    features = pd.DataFrame(index=df.index)
    # Features calendaires
    features["week_of_year"] = df.index.isocalendar().week.astype(int)
    features["month"] = df.index.month
    features["quarter"] = df.index.quarter
    # Fourier terms (saisonnalite annuelle)
    features["sin_52"] = np.sin(2 * np.pi * features["week_of_year"] / 52)
    features["cos_52"] = np.cos(2 * np.pi * features["week_of_year"] / 52)
    features["sin_26"] = np.sin(2 * np.pi * features["week_of_year"] / 26)
    features["cos_26"] = np.cos(2 * np.pi * features["week_of_year"] / 26)
    # Lags
    y = df[target_col]
    features["lag_1"] = y.shift(1)
    features["lag_2"] = y.shift(2)
    features["lag_4"] = y.shift(4)
    features["lag_52"] = y.shift(52)
    # Rolling stats
    features["rolling_mean_4"] = y.shift(1).rolling(window=4).mean()
    features["rolling_std_4"] = y.shift(1).rolling(window=4).std()
    features["rolling_mean_12"] = y.shift(1).rolling(window=12).mean()
    # Variables exogenes
    features["temperature"] = df["temperature"]
    features["jours_ouvres"] = df["jours_ouvres"]
    features["vacances"] = df["vacances"]
    return features

What does this function do?

La fonction create_features construit 4 familles de features a partir de la serie brute :

Features calendaires (week_of_year, month, quarter) pour capter la position dans l'annee.
Termes de Fourier (sin/cos a periodicites 52 et 26 semaines) pour modeliser la saisonnalite de facon continue.
Lags (1, 2, 4, 52 semaines) et rolling stats (moyenne, ecart-type sur 4 et 12 semaines) pour capturer l'auto-correlation.
Les variables exogenes (temperature, jours ouvres, vacances) deja presentes dans le DataFrame.

Dans la meme cellule, juste apres la fonction, ajoute le code de preparation des features et d'entrainement du modele :

# Preparer les features sur tout le dataset
all_features = create_features(df)
all_features = all_features.dropna()
df_aligned = df.loc[all_features.index]

# Re-split apres feature engineering
split_date = train.index[-1]
X_train_xgb = all_features.loc[all_features.index <= split_date]
X_test_xgb = all_features.loc[all_features.index > split_date]
y_train_xgb = df_aligned.loc[X_train_xgb.index, "remplacements"]
y_test_xgb = df_aligned.loc[X_test_xgb.index, "remplacements"]

# Entrainement XGBoost
print("\nEntrainement XGBoost...")
xgb_model = XGBRegressor(
    n_estimators=200,
    max_depth=5,
    learning_rate=0.05,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42
)
xgb_model.fit(X_train_xgb, y_train_xgb, verbose=False)
xgb_forecast = pd.Series(xgb_model.predict(X_test_xgb), index=X_test_xgb.index)
print("XGBoost entraine.")

Why dropna and re-split?

Les lags et rolling stats introduisent des valeurs manquantes en debut de serie (les 52 premieres semaines n'ont pas de lag_52). Le dropna les supprime, puis on re-split train/test en utilisant split_date pour garder la coherence temporelle.

Execute cette cellule.

Tu devrais voir s'afficher "XGBoost entraine." dans la sortie.

Seeing a NameError on create_features?

Assure-toi que la definition de la fonction et le code d'entrainement sont dans la meme cellule. Si tu les as separes, execute d'abord la cellule contenant la fonction.

help me debug a NameError on create_features in my XGBoost cell

Compare all forecasts visually and numerically

Tu as maintenant 4 forecasts. Le graphique superpose te permet de voir immediatement quel modele suit le mieux les variations reelles.

Cree une nouvelle cellule et ajoute le graphique de comparaison des forecasts :

# %% Visualisation des forecasts superposes
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(y_train.index, y_train, color="steelblue", label="Train", linewidth=0.8)
ax.plot(y_test.index, y_test, color="black", label="Reel (test)", linewidth=1.5)
ax.plot(sarima_forecast.index, sarima_forecast, "--", color="red", label="SARIMA", alpha=0.8)
ax.plot(sarimax_forecast.index, sarimax_forecast, "--", color="green", label="SARIMAX", alpha=0.8)
ax.plot(ets_forecast.index, ets_forecast, "--", color="purple", label="ETS", alpha=0.8)
ax.plot(xgb_forecast.index, xgb_forecast, "--", color="darkorange", label="XGBoost", alpha=0.8)
ax.axvline(x=split_date, color="gray", linestyle=":", alpha=0.7)
ax.legend()
ax.set_title("Comparaison des Forecasts - Remplacement Linky")
ax.set_ylabel("Nombre de remplacements")
plt.tight_layout()
plt.show()

Execute la cellule.

Tu devrais voir un graphique avec la serie d'entrainement en bleu, le reel en noir, et les 4 forecasts en couleurs pointillees. SARIMAX (vert) et XGBoost (orange) suivent en general mieux les variations que SARIMA seul.

Pour quantifier cette observation, tu as besoin de metriques numeriques. La cellule suivante calcule MAE, RMSE et MAPE pour chaque modele.

Cree une nouvelle cellule et ajoute la fonction de metriques et le tableau comparatif :

# %% Metriques de comparaison
def compute_metrics(y_true, y_pred, model_name):
    """Calcule MAE, RMSE et MAPE."""
    mae = mean_absolute_error(y_true, y_pred)
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
    return {"Modele": model_name, "MAE": round(mae, 2), "RMSE": round(rmse, 2), "MAPE (%)": round(mape, 2)}

# Aligner les index pour la comparaison
common_idx = y_test.index.intersection(xgb_forecast.index)
metrics = []
metrics.append(compute_metrics(y_test.loc[common_idx], sarima_forecast.loc[common_idx], "SARIMA"))
metrics.append(compute_metrics(y_test.loc[common_idx], sarimax_forecast.loc[common_idx], "SARIMAX"))
metrics.append(compute_metrics(y_test.loc[common_idx], ets_forecast.loc[common_idx], "ETS"))
metrics.append(compute_metrics(y_test.loc[common_idx], xgb_forecast.loc[common_idx], "XGBoost"))
metrics_df = pd.DataFrame(metrics).set_index("Modele")
print("\n" + "=" * 50)
print("COMPARAISON DES MODELES")
print("=" * 50)
print(metrics_df.to_string())

What does this code do?

La fonction compute_metrics calcule trois metriques complementaires pour chaque modele :

MAE (Mean Absolute Error) : erreur moyenne en valeur absolue, facile a interpreter.
RMSE (Root Mean Squared Error) : penalise davantage les grosses erreurs.
MAPE (Mean Absolute Percentage Error) : erreur en pourcentage, comparable entre series d'echelles differentes.

Le common_idx aligne les index car XGBoost perd quelques semaines a cause du dropna sur les lags.

Execute la cellule.

Tu devrais voir un tableau affichant MAE, RMSE et MAPE pour les 4 modeles. SARIMAX et XGBoost devraient avoir les meilleurs scores.

Seeing a KeyError on common_idx?

Verifie que tu as bien execute la cellule XGBoost precedente et que xgb_forecast existe. Si tu as relance le kernel, re-execute toutes les cellules dans l'ordre.

help me debug a KeyError with common_idx in my metrics cell

Implement walk-forward cross-validation

Le tableau precedent evalue les modeles sur un seul split train/test. Ce n'est pas suffisant pour conclure. La walk-forward cross-validation re-entraine chaque modele sur des fenetres croissantes et mesure la stabilite des performances.

Why walk-forward instead of random CV?

En series temporelles, le futur ne doit jamais fuiter dans l'entrainement. Le walk-forward respecte l'ordre chronologique : on entraine sur le passe, on teste sur le futur immediat, puis on avance la fenetre.

L'expanding window utilise tout l'historique disponible. Le sliding window garde une taille fixe (ici 104 semaines). Comparer les deux revele si les donnees anciennes aident ou polluent les predictions.

Cree une nouvelle cellule et commence la fonction walk_forward_cv avec les branches SARIMA et SARIMAX :

# %% Walk-Forward Cross-Validation
def walk_forward_cv(df, model_type, n_splits=4, window_type="expanding"):
    """Walk-forward CV avec expanding ou sliding window."""
    tscv = TimeSeriesSplit(n_splits=n_splits)
    y = df["remplacements"]
    exog = df[["temperature", "jours_ouvres", "vacances"]]
    scores = []
    for fold, (train_idx, test_idx) in enumerate(tscv.split(y)):
        y_tr = y.iloc[train_idx]
        y_te = y.iloc[test_idx]
        exog_tr = exog.iloc[train_idx]
        exog_te = exog.iloc[test_idx]
        if window_type == "sliding":
            min_size = min(104, len(train_idx))
            y_tr = y_tr.iloc[-min_size:]
            exog_tr = exog_tr.iloc[-min_size:]
        try:
            if model_type == "SARIMA":
                model = auto_arima(y_tr, seasonal=True, m=52, d=1, D=1,
                                   stepwise=True, suppress_warnings=True,
                                   error_action="ignore", max_p=2, max_q=2,
                                   max_P=1, max_Q=1)
                pred = model.predict(n_periods=len(y_te))
            elif model_type == "SARIMAX":
                model = SARIMAX(y_tr, exog=exog_tr, order=(1,1,1),
                               seasonal_order=(1,1,0,52),
                               enforce_stationarity=False,
                               enforce_invertibility=False)
                fit = model.fit(disp=False)
                pred = fit.forecast(steps=len(y_te), exog=exog_te)

Continue dans la meme cellule avec les branches ETS et XGBoost, puis la fin de la fonction :

            elif model_type == "ETS":
                model = ExponentialSmoothing(y_tr, trend="add", seasonal="add",
                                           seasonal_periods=52,
                                           initialization_method="estimated")
                fit = model.fit(optimized=True)
                pred = fit.forecast(steps=len(y_te))
            elif model_type == "XGBoost":
                features_all = create_features(df)
                features_all = features_all.dropna()
                valid_train = features_all.index.isin(y_tr.index)
                valid_test = features_all.index.isin(y_te.index)
                X_tr = features_all.loc[valid_train]
                X_te = features_all.loc[valid_test]
                y_tr_xgb = df.loc[X_tr.index, "remplacements"]
                y_te_xgb = df.loc[X_te.index, "remplacements"]
                if len(X_tr) == 0 or len(X_te) == 0:
                    continue
                xgb = XGBRegressor(n_estimators=200, max_depth=5,
                                   learning_rate=0.05, random_state=42)
                xgb.fit(X_tr, y_tr_xgb, verbose=False)
                pred = xgb.predict(X_te)
                y_te = y_te_xgb
            rmse = np.sqrt(mean_squared_error(y_te, pred))
            scores.append(rmse)
        except Exception:
            continue
    return scores

What does this function do?

La fonction walk_forward_cv utilise TimeSeriesSplit pour creer 4 folds chronologiques. A chaque fold, elle entraine le modele demande sur le passe et mesure le RMSE sur le futur.

En mode sliding, elle tronque l'historique a 104 semaines (2 ans). Cela teste si les donnees tres anciennes aident ou degradent la prediction.

Juste apres la fonction (toujours dans la meme cellule), ajoute les boucles qui executent la CV pour les 4 modeles :

print("\nWalk-Forward CV (expanding window)...")
results_expanding = {}
for model_name in ["SARIMA", "SARIMAX", "ETS", "XGBoost"]:
    scores = walk_forward_cv(df, model_name, n_splits=4, window_type="expanding")
    results_expanding[model_name] = scores
    print(f"  {model_name}: RMSE moyen = {np.mean(scores):.2f} (+/- {np.std(scores):.2f})")

print("\nWalk-Forward CV (sliding window)...")
results_sliding = {}
for model_name in ["SARIMA", "SARIMAX", "ETS", "XGBoost"]:
    scores = walk_forward_cv(df, model_name, n_splits=4, window_type="sliding")
    results_sliding[model_name] = scores
    print(f"  {model_name}: RMSE moyen = {np.mean(scores):.2f} (+/- {np.std(scores):.2f})")

Why does this cell take a while?

La cellule re-entraine chaque modele 4 fois (4 folds) x 2 strategies (expanding + sliding) = 8 entrainements par modele. SARIMA avec auto_arima est le plus lent car il cherche les meilleurs ordres a chaque fold. Attends 1 a 3 minutes selon ta machine.

Execute la cellule. Attends que les 8 lignes de resultats s'affichent (4 modeles x 2 strategies).

Tu devrais voir le RMSE moyen et l'ecart-type pour chaque modele, d'abord en expanding window puis en sliding window.

Cell seems stuck or very slow?

C'est normal que SARIMA prenne du temps avec m=52. Si ca depasse 5 minutes, interromps la cellule et reduis n_splits a 3 temporairement pour valider que le code fonctionne.

help me speed up my walk-forward CV execution

Visualize the sensitivity to validation strategy

Les scores numeriques sont la, mais un graphique revele immediatement quels modeles sont stables et lesquels sont sensibles au choix de la fenetre.

Cree une nouvelle cellule et ajoute le bar chart comparant expanding vs sliding window :

# %% Visualisation expanding vs sliding
fig, ax = plt.subplots(figsize=(10, 5))
models = list(results_expanding.keys())
x = np.arange(len(models))
width = 0.35
expanding_means = [np.mean(results_expanding[m]) for m in models]
sliding_means = [np.mean(results_sliding[m]) for m in models]
bars1 = ax.bar(x - width/2, expanding_means, width, label="Expanding Window", color="steelblue")
bars2 = ax.bar(x + width/2, sliding_means, width, label="Sliding Window", color="darkorange")
ax.set_xticks(x)
ax.set_xticklabels(models)
ax.set_ylabel("RMSE moyen (walk-forward CV)")
ax.set_title("Impact de la strategie de validation sur la performance")
ax.legend()
plt.tight_layout()
plt.show()

What does this chart reveal?

Chaque paire de barres compare le RMSE moyen d'un modele en expanding (bleu) vs sliding (orange) window. Si un modele performe mieux en sliding, c'est que les donnees anciennes le perturbent. Si expanding est meilleur, l'historique complet aide la prediction.

A ton avis, les modeles performent-ils mieux avec plus d'historique (expanding) ou avec une fenetre fixe (sliding)?

Execute la cellule.

Tu devrais voir un bar chart avec 4 paires de barres. Les modeles statistiques (SARIMA, ETS) montrent generalement peu de difference entre les deux strategies. XGBoost peut etre plus sensible au volume de donnees d'entrainement.

Chart shows empty or zero-height bars?

Verifie que la cellule walk-forward CV precedente a termine sans erreur et que results_expanding et results_sliding contiennent des scores non-vides. Affiche-les avec un print pour diagnostiquer.

help me debug empty walk-forward CV results

✔️ Awesome, I've got everything!

Parfait. Verifie que tu as bien sauvegarde ton notebook.

ⓧ I'd like to double check the full code

Scroll up through your notebook and confirm these cells are present in order :

Cellule XGBoost - Feature Engineering : contient create_features, all_features, split, XGBRegressor training, et xgb_forecast.
Cellule Visualisation des forecasts superposes : contient le plot avec 4 modeles en pointilles + la serie reelle en noir.
Cellule Metriques de comparaison : contient compute_metrics, common_idx, et le print de metrics_df.
Cellule Walk-Forward Cross-Validation : contient walk_forward_cv (la fonction entiere), les boucles expanding et sliding, et les prints de resultats.
Cellule Visualisation expanding vs sliding : contient le bar chart comparatif.

Tu as maintenant un pipeline de forecasting complet avec 4 modeles compares rigoureusement. Tu sais quel modele performe le mieux et si le choix de la strategie de validation change le classement. Dans la Secret Mission, tu iras plus loin en testant si cette difference est statistiquement significative.

Secret mission

Test de Diebold-Mariano et Ensembling

Ton tableau de metriques montre un classement des modeles, mais peut-on prouver statistiquement qu'un modele bat l'autre ? Le test de Diebold-Mariano repond a cette question. Combine-le avec un ensemble pondere pour potentiellement surpasser tout modele individuel.

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs or cloud resources to worry about. Decide whether to keep your files and virtual environment in place, or delete them to free up disk space.

Resources you used:

Virtual environment (.venv folder) with installed Python packages (statsmodels, pmdarima, xgboost, scikit-learn).
Notebook file linky_forecast_pipeline.py on your local disk.

Keep everything running

No action needed. Choose this if you plan to keep iterating on your forecasting pipeline or want to experiment with different models and parameters.

Your virtual environment stays activated and ready to use.
You can reopen the notebook in VS Code at any time and re-run cells.

Pause - I'll come back to this later

Deactivate the virtual environment to free up your terminal session. Your files and packages remain on disk for when you return.

In your terminal, deactivate the virtual environment by running:

deactivate

When you come back, navigate to your project folder and reactivate with source .venv/bin/activate (macOS/Linux) or .venv\Scripts\activate (Windows).

Delete - I don't want to use this again

Remove the virtual environment and project files entirely. This frees up disk space and removes all installed packages.

Deactivate the virtual environment first by running:

deactivate

Navigate to your Desktop by running:

cd ~/Desktop

Delete the project folder (macOS/Linux) by running:

rm -rf linky-forecast

On Windows, delete the project folder by running:

rmdir /s /q linky-forecast

This removes the .venv folder (with all installed packages), the notebook linky_forecast_pipeline.py, and any generated output.

Nice Work!

You just built a complete weekly forecasting pipeline for Linky meter replacements. From synthetic data generation to model comparison with rigorous validation, you now have a reproducible notebook you can plug into real operational data by swapping a single CSV.

You've learned how to:

Build a reproducible forecasting pipeline that generates realistic synthetic data, runs statistical diagnostics, trains 4 model families (SARIMA, SARIMAX, ETS, XGBoost), and outputs a ranked comparison table.
Prove non-stationarity and seasonal structure using ADF, KPSS, STL decomposition, and ACF/PACF analysis before fitting any model.
Stress-test model robustness with walk-forward cross-validation (expanding and sliding window) and quantify how sensitive each model is to the train/test cutoff date.
Secret Mission: implement a Diebold-Mariano test from scratch and build a weighted ensemble that outperforms both individual best models.

Ready to quiz yourself?