Forecast Linky Maintenance: 5 Models

Benchmark 5 forecasting models on weekly maintenance data through 2028.

Introduction

30 Second Summary

Imagine knowing exactly how many maintenance operations your team will face next week, next month, or next year. That kind of foresight turns reactive scrambling into calm, scheduled planning.

In this project, you will build a complete forecasting pipeline that benchmarks 5 different models on your weekly Linky concentrator maintenance data, ranks them by accuracy, and generates operational predictions through December 2028.

What You'll Build

You will run your notebook, point to any future week, and confidently say "we expect X maintenance operations, plus or minus Y" with full statistical backing ready to present to planning teams.

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

  • A seasonal decomposition revealing yearly patterns in your maintenance data, plus side-by-side model fits showing how ETS, SARIMA, SARIMAX, Prophet, and XGBoost each capture the signal.
  • A metrics dashboard ranking all 5 models by MAE, RMSE, and MAPE, with residual diagnostic plots proving whether each model adequately captured the structure.
  • A forecast chart through 2028 with 90% confidence intervals from all models, plus a business interpretation example translating a forecast week into actionable planning recommendations.
  • Secret Mission: Build an ensemble model that combines the best-performing models with optimized weights to outperform any individual approach.

Are there any prerequisites?

You need Python 3.12+ installed, Jupyter Notebook installed, and basic familiarity with pandas and time series concepts.

Before We Start

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

Set Up Your Environment and Load Data

Forecasting maintenance volumes requires several specialized Python libraries working together. You need statistical models, machine learning tools, and French calendar data packages to build a meaningful benchmark.

In this step, you will install all required packages, create your notebook, load your weekly maintenance data, and run an initial decomposition to reveal the underlying signal structure.

In this step, get ready to:
  • Install all forecasting and calendar packages.
  • Create your notebook file and import all libraries.
  • Load your data and visualize the seasonal decomposition.
Install the forecasting packages

This project uses five distinct forecasting approaches, each from a different library. You also need two French calendar packages to inject holiday and school-break effects later.

  • Install all required packages by running this command in your terminal:
pip install statsmodels==0.14.6 prophet==1.3.0 xgboost==3.3.0 jours-feries-france==0.7.0 vacances-scolaires-france==0.11.0 scikit-learn matplotlib pandas numpy

What does this command do?

This installs the exact pinned versions of each library your project needs. statsmodels provides ETS and SARIMA models. prophet is Facebook's forecasting tool with built-in holiday support. xgboost is a gradient-boosted tree library for non-linear modeling.

The two calendar packages (jours-feries-france and vacances-scolaires-france) provide French public holiday and school holiday data you will use as exogenous regressors in later steps.

  • Verify the core libraries installed successfully by running:
python -c "import statsmodels; import prophet; import xgboost; print('OK')"
  • You should see OK printed with no errors.

Seeing an import error?

If prophet fails to import, it may need pystan as a dependency. Run pip install pystan then retry. If xgboost fails, confirm you are running Python 3.12+ since xgboost==3.3.0 requires it.

Still stuck? Help me fix my package installation errors.

Create your notebook and import libraries

Your forecasting code lives in a single Python file that uses # %% cell markers. This format lets Jupyter Notebook and VS Code treat each section as an executable cell.

  • Create a new file called forecast_maintenance.py in a working folder on your Desktop.
  • Paste the following imports and configuration as the first cell:
# %% [markdown]
# # Prevision Hebdomadaire - Maintenance Concentrateurs Linky
# Benchmark: ETS, SARIMA, SARIMAX, Prophet, XGBoost

# %% Imports and Configuration
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')

from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungbox

from prophet import Prophet

from xgboost import XGBRegressor

from sklearn.metrics import mean_absolute_error, mean_squared_error

from jours_feries_france import JoursFeries
from vacances_scolaires_france import SchoolHolidayDates

# Configuration
first_date = '2023-01-01'
last_forecast = '2028-12-31'

What does this code do?

  • The # %% [markdown] line starts a documentation cell with the notebook title.
  • The # %% Imports and Configuration line starts an executable cell that loads every library you will use across all five models.
  • The two configuration variables first_date and last_forecast define your training window start and the forecast horizon end.
  • Run this cell in your notebook environment.
  • You should see no errors. If a library fails to import, revisit the installation step above.
Load your data and visualize the signal

With all libraries loaded, you can bring in your weekly maintenance data and immediately see what patterns exist. The seasonal decomposition will reveal whether your data has a yearly cycle, a trend, or unexplained residual noise.

  • Add a new cell below your imports and paste the data loading code:
# %% Data Loading
# Adapt this section to your data source (CSV, Excel, etc.)
# df = pd.read_csv('your_file.csv', index_col='WEEK', parse_dates=True)
# df = df.loc[first_date:]
# df = df.asfreq('W-MON')  # Adjust to your week start day

# Placeholder - replace with your actual data loading
# df must have a DatetimeIndex (weekly) and a column named 'Maintenance'
print(f"Data shape: {df.shape}")
print(f"Date range: {df.index[0]} to {df.index[-1]}")
print(df.describe())
df.plot(figsize=(14, 4), title='Weekly Maintenance Operations')
plt.tight_layout()
plt.show()

How do I load my own data?

Uncomment the pd.read_csv line (or use pd.read_excel) and point it at your actual file. Your DataFrame must have a DatetimeIndex with weekly frequency and a column named Maintenance.

The .loc[first_date:] filter trims data to start from 2023-01-01. The .asfreq('W-MON') call ensures pandas recognizes the weekly frequency. Adjust W-MON to match your week start day.

  • Run this cell after adapting the loading to your file.
  • You should see the data shape, date range starting at 2023-01-01, descriptive statistics, and a time series line plot of your weekly maintenance operations.

Seeing a NameError for df?

If you get NameError: name 'df' is not defined, it means the commented-out loading lines have not been replaced with your actual data source. Uncomment and edit the pd.read_csv line to point to your real file path.

Need help loading your data? Help me load my CSV into a DataFrame with a weekly DatetimeIndex.

  • Add a new cell for the seasonal decomposition:
# %% Seasonal Decomposition
decomposition = seasonal_decompose(df['Maintenance'], model='additive', period=52)
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.tight_layout()
plt.show()

What does the decomposition show?

The seasonal_decompose function splits your time series into three components: trend (long-term direction), seasonal (repeating 52-week pattern), and residual (leftover noise).

Setting period=52 tells the algorithm that one full cycle is 52 weeks (one year). If your maintenance data drops during holidays, that pattern will show up clearly in the seasonal component.

  • Run this cell.
  • You should see a four-panel chart with observed data, trend, seasonal pattern, and residuals.
  • Add a final cell that splits the data for evaluation and defines a reusable metrics helper:
# %% Train/Test Split
test_size = 12  # last 12 weeks
train = df.iloc[:-test_size]
test = df.iloc[-test_size:]
print(f"Train: {train.index[0]} to {train.index[-1]} ({len(train)} weeks)")
print(f"Test:  {test.index[0]} to {test.index[-1]} ({len(test)} weeks)")


# %% Helper: Metrics
def compute_metrics(actual, predicted, model_name):
    mae = mean_absolute_error(actual, predicted)
    rmse = np.sqrt(mean_squared_error(actual, predicted))
    mape = np.mean(np.abs((actual - predicted) / actual)) * 100
    return {'Model': model_name, 'MAE': round(mae, 2), 'RMSE': round(rmse, 2), 'MAPE': round(mape, 2)}


results = []

What does this code do?

  • The train/test split holds back the last 12 weeks for evaluation. Every model will be tested against these same 12 weeks so comparisons are fair.
  • The compute_metrics function calculates three metrics (MAE, RMSE, MAPE) for any model's predictions. You will call it for each model to build the comparison table.
  • The empty results list will accumulate each model's scores as you train them in later steps.
  • Run this cell.

Before you check the output, what do you expect the printed train and test ranges to show?

  • You should see two lines: the training set covering from your start date through 12 weeks before the end, and the test set covering the final 12 weeks.

Train set shows fewer weeks than expected?

If the numbers look wrong, confirm your data loading step applied .asfreq('W-MON') correctly. Missing frequency metadata causes iloc to count rows but dates may not align to weekly boundaries.

Help me fix my train/test split.

✔️ Awesome, I've got everything!

Great. Double check you have saved forecast_maintenance.py with all four cells (imports, data loading, decomposition, train/test split).

ⓧ I'd like to double check the full code

Your file should contain these four cell blocks in order. Scroll up and confirm each one is present:

  • Cell 1: Imports and Configuration (starting with # %% [markdown] and ending with last_forecast = '2028-12-31').
  • Cell 2: Data Loading (starting with # %% Data Loading and ending with plt.show()).
  • Cell 3: Seasonal Decomposition (starting with # %% Seasonal Decomposition and ending with plt.show()).
  • Cell 4: Train/Test Split and Helper (starting with # %% Train/Test Split and ending with results = []).

Your environment is ready and your data is split for evaluation. Next up, you will fit the first two baseline models and see how pure statistical approaches handle your weekly maintenance series.

Build ETS and SARIMA Baseline Models

Your seasonal decomposition revealed the underlying patterns in the maintenance data. Now the question is: can a model capture those patterns well enough to forecast future weeks?

In this step, you'll fit two classical statistical models as baselines. If they perform well, complex models may be unnecessary. If they struggle, especially around holidays, it proves the need for calendar-aware approaches in the next step.

In this step, get ready to:
  • Fit an ETS model with additive seasonality and evaluate its 12-week forecast against the test set.
  • Fit a SARIMA model and compare its metrics to ETS.
  • Run residual diagnostics to assess what the baseline models missed.
Fit the ETS baseline model

ETS stands for Error, Trend, Seasonality. It decomposes your time series into these three components and uses them to generate forecasts. You already have your train and test DataFrames ready from Step 1.

  • Create a new cell in your notebook and fit the ETS model by running this code:
# %% ETS Model
ets_model = ETSModel(
    train['Maintenance'],
    error='add',
    trend='add',
    seasonal='add',
    seasonal_periods=52
)
ets_fit = ets_model.fit(disp=False)
print(ets_fit.summary())

ets_pred = ets_fit.forecast(steps=test_size)
results.append(compute_metrics(test['Maintenance'].values, ets_pred.values, 'ETS'))

plt.figure(figsize=(14, 4))
plt.plot(train.index, train['Maintenance'], label='Train')
plt.plot(test.index, test['Maintenance'], label='Test', color='orange')
plt.plot(test.index, ets_pred.values, label='ETS Forecast', color='green', linestyle='--')
plt.legend()
plt.title('ETS - Fit vs Actual')
plt.tight_layout()
plt.show()

What does this code do?

  • ETSModel with all-additive components tells the model to look for a linear trend and repeating seasonal pattern every 52 weeks (one year).
  • .forecast(steps=test_size) generates predictions for the 12 held-out weeks, so you can measure accuracy on data the model never saw.
  • compute_metrics calculates MAE, RMSE, and MAPE, then appends the results to your results list for comparison later.
  • Run the cell. You should see a printed summary table of ETS parameters, followed by a plot showing the green dashed forecast line against the orange test data.

Seeing an error about seasonal_periods?

Make sure your train DataFrame has at least 52 rows. ETS needs at least one full seasonal cycle to estimate the pattern. If your data starts later than 2023-01-01 or your total data is too short, adjust first_date accordingly.

Help me debug my ETS model fitting error

Fit the SARIMA model

SARIMA extends ARIMA with a seasonal component. It models both short-term autocorrelation and repeating seasonal patterns through differencing and lagged terms.

  • Create a new cell and fit the SARIMA model by running this code:
# %% SARIMA Model
sarima_model = SARIMAX(
    train['Maintenance'],
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 52),
    enforce_stationarity=False,
    enforce_invertibility=False
)
sarima_fit = sarima_model.fit(disp=False)
print(sarima_fit.summary())

sarima_pred = sarima_fit.forecast(steps=test_size)
results.append(compute_metrics(test['Maintenance'].values, sarima_pred.values, 'SARIMA'))

plt.figure(figsize=(14, 4))
plt.plot(train.index, train['Maintenance'], label='Train')
plt.plot(test.index, test['Maintenance'], label='Test', color='orange')
plt.plot(test.index, sarima_pred.values, label='SARIMA Forecast', color='red', linestyle='--')
plt.legend()
plt.title('SARIMA - Fit vs Actual')
plt.tight_layout()
plt.show()

What does this code do?

  • order=(1, 1, 1) means one autoregressive lag, one difference, and one moving average term for short-term patterns.
  • seasonal_order=(1, 1, 1, 52) adds the same structure at the 52-week seasonal lag. This captures year-over-year repeating patterns.
  • enforce_stationarity=False and enforce_invertibility=False relax numeric constraints that can prevent convergence on real-world data.

SARIMA with s=52 is slow

Fitting SARIMA with a seasonal period of 52 is computationally intensive. Expect this cell to run for 1 to 5 minutes depending on your machine. The progress bar may not appear, so be patient.

  • Run the cell. After it finishes, you should see a printed summary table and a plot with the red dashed SARIMA forecast line against the orange test data.

Cell running for more than 5 minutes?

SARIMA with seasonal_order=(1, 1, 1, 52) on weekly data can take a while. If it exceeds 10 minutes, try restarting your kernel and running all cells from the top. Make sure no other heavy processes are consuming CPU.

Help me troubleshoot my slow SARIMA fit

Plot residual diagnostics for both models

Model metrics tell you how far off the predictions are. Residual diagnostics tell you why. If residuals show remaining patterns (autocorrelation, skew), the model is leaving signal on the table.

  • Create a new cell and plot the residual diagnostics for both models by running this code:
# %% Residual Diagnostics - ETS and SARIMA
fig, axes = plt.subplots(2, 3, figsize=(15, 8))

# ETS residuals
ets_resid = ets_fit.resid
axes[0, 0].hist(ets_resid, bins=20, edgecolor='black')
axes[0, 0].set_title('ETS - Residual Distribution')
from scipy import stats
stats.probplot(ets_resid, plot=axes[0, 1])
axes[0, 1].set_title('ETS - Q-Q Plot')
plot_acf(ets_resid, ax=axes[0, 2], lags=40)
axes[0, 2].set_title('ETS - Residual ACF')

# SARIMA residuals
sarima_resid = sarima_fit.resid
axes[1, 0].hist(sarima_resid, bins=20, edgecolor='black')
axes[1, 0].set_title('SARIMA - Residual Distribution')
stats.probplot(sarima_resid, plot=axes[1, 1])
axes[1, 1].set_title('SARIMA - Q-Q Plot')
plot_acf(sarima_resid, ax=axes[1, 2], lags=40)
axes[1, 2].set_title('SARIMA - Residual ACF')

plt.tight_layout()
plt.show()

What do these three diagnostic plots show?

  • Histogram shows whether residuals are roughly normally distributed. Skew or heavy tails suggest the model is systematically over- or under-predicting certain periods.
  • Q-Q plot compares residual quantiles to a theoretical normal distribution. Points deviating from the diagonal indicate non-normality.
  • ACF plot shows autocorrelation in the residuals. Significant spikes mean the model left predictable structure uncaptured. This is the most important one for deciding if you need a better model.

Before you run this cell, what do you expect the ACF plots to show? Do you think these baseline models fully captured the seasonal signal, or is there remaining autocorrelation?

  • Run the cell. You should see a 2x3 grid of diagnostic plots. Check the ACF columns (rightmost) for significant spikes beyond the blue confidence band.

If you see spikes in the ACF at holiday-adjacent lags, that is evidence these models cannot account for calendar effects. This is the baseline quality that the calendar-enhanced models in the next step need to beat.

Seeing a scipy import error?

The from scipy import stats line imports scipy which ships with most scientific Python distributions. If you see an import error, install it by running pip install scipy in your terminal, then restart your kernel.

Help me fix my scipy import error in Jupyter

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all three cells (ETS Model, SARIMA Model, Residual Diagnostics) present and run successfully.

ⓧ I'd like to double check the full code

Scroll up through your notebook and confirm these three cells are present in order, each one run with output visible:

  • The # %% ETS Model cell with ETSModel fitting, forecast, metrics append, and plot.
  • The # %% SARIMA Model cell with SARIMAX fitting, forecast, metrics append, and plot.
  • The # %% Residual Diagnostics - ETS and SARIMA cell with the 2x3 grid of histogram, Q-Q, and ACF plots for both models.

You now have two baseline models with concrete metrics. Both models treat every week the same. They have no awareness that holiday weeks reduce maintenance activity or that school break periods shift workload patterns. In the next step, you'll add French calendar effects and see if that closes the gap in residual autocorrelation.

Add Calendar Effects with SARIMAX and Prophet

Your ETS and SARIMA baselines treat every week identically. But maintenance activity on Linky concentrators drops predictably during public holidays and school vacation periods.

In this step, you'll build calendar features from two French reference libraries, feed them into a SARIMAX model as exogenous regressors, and then fit a Prophet model that handles holidays natively. If calendar effects matter, these models should beat the baselines.

In this step, get ready to:
  • Build calendar features (public holidays and school vacation proportions) for each week.
  • Fit a SARIMAX model with calendar exogenous variables and evaluate it.
  • Fit a Prophet model with built-in French holiday handling and evaluate it.
Build calendar features for each week

The baseline models have no idea that Christmas week or August vacation behaves differently from a normal week. You need to encode that knowledge as numeric features the model can learn from.

  • In your notebook, create a new cell below the residual diagnostics section and define the build_calendar_features function by pasting this code:
# %% Calendar Features Construction
def build_calendar_features(date_index):
    """Build weekly calendar features: nb jours feries + pct vacances scolaires."""
    school_holidays = SchoolHolidayDates()
    records = []
    for week_start in date_index:
        week_dates = pd.date_range(week_start, periods=7, freq='D')
        # Count public holidays in the week
        nb_feries = sum(
            1 for d in week_dates
            if JoursFeries.is_bank_holiday(d.date(), zone="Metropole")
        )
        # Proportion of school holiday days
        nb_vacances = sum(
            1 for d in week_dates
            if school_holidays.is_holiday(d.date())
        )
        pct_vacances = nb_vacances / 7.0
        records.append({'nb_jours_feries': nb_feries, 'pct_vacances': pct_vacances})
    return pd.DataFrame(records, index=date_index)

What does this code do?

For each week in the date index, the function iterates over the 7 days and counts two things.

  • nb_jours_feries: how many public holidays (jours fériés) fall in that week, using the jours-feries-france library.
  • pct_vacances: what proportion of the week's days fall during school holidays (vacances scolaires), using vacances-scolaires-france.

These two numbers give the model a way to learn that weeks with holidays have lower maintenance activity.

  • Now call this function for your train, test, full dataset, and future forecast horizon by adding another cell with this code:
# Build for training period
calendar_train = build_calendar_features(train.index)
calendar_test = build_calendar_features(test.index)
calendar_all = build_calendar_features(df.index)

# Build for forecast horizon
future_index = pd.date_range(
    start=df.index[-1] + pd.Timedelta(weeks=1),
    end=last_forecast,
    freq=df.index.freq
)
calendar_future = build_calendar_features(future_index)

print("Calendar features sample:")
print(calendar_all.head(10))

What does this code do?

  • calendar_train, calendar_test, and calendar_all align calendar features to your existing weekly indexes.
  • future_index creates a weekly DatetimeIndex from the end of your data through 2028-12-31 for generating forecasts later.
  • calendar_future builds the same features for the forecast horizon so SARIMAX can use them when predicting ahead.
  • Run both cells in order.

You should see a table showing 10 rows with columns nb_jours_feries and pct_vacances. Weeks containing a public holiday will show values of 1 or higher in the first column.

Seeing an error about the zone parameter?

Some versions of jours-feries-france expect the accented spelling "Métropole" instead of "Metropole". If you see an error related to the zone, try switching to the accented version in the function.

help me fix the zone parameter error in build_calendar_features

Fit SARIMAX with calendar regressors

SARIMAX extends SARIMA by accepting exogenous variables. By passing calendar_train as exog, the model learns how holiday counts and vacation proportions shift the maintenance level up or down.

  • Create a new cell and fit the SARIMAX model with calendar effects by pasting this code:
# %% SARIMAX with Calendar Effects
sarimax_model = SARIMAX(
    train['Maintenance'],
    exog=calendar_train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 52),
    enforce_stationarity=False,
    enforce_invertibility=False
)
sarimax_fit = sarimax_model.fit(disp=False)
print(sarimax_fit.summary())

sarimax_pred = sarimax_fit.forecast(steps=test_size, exog=calendar_test)
results.append(compute_metrics(test['Maintenance'].values, sarimax_pred.values, 'SARIMAX'))

What does this code do?

  • The SARIMAX constructor takes the same order=(1,1,1) and seasonal_order=(1,1,1,52) as the pure SARIMA, but adds exog=calendar_train to incorporate calendar effects.
  • The summary will show coefficients and p-values for nb_jours_feries and pct_vacances. Low p-values mean the calendar features are statistically significant.
  • When forecasting, you pass exog=calendar_test so the model knows which future weeks contain holidays.
  • Run the cell. This may take 1-3 minutes due to the seasonal period of 52 weeks.

You should see the SARIMAX summary table printed. Look at the coefficients section for your exogenous regressors. Significant p-values (below 0.05) confirm that holidays affect maintenance volumes.

  • Now plot the SARIMAX forecast against actuals and run residual diagnostics by adding this cell:
plt.figure(figsize=(14, 4))
plt.plot(train.index, train['Maintenance'], label='Train')
plt.plot(test.index, test['Maintenance'], label='Test', color='orange')
plt.plot(test.index, sarimax_pred.values, label='SARIMAX Forecast', color='purple', linestyle='--')
plt.legend()
plt.title('SARIMAX (Calendar Effects) - Fit vs Actual')
plt.tight_layout()
plt.show()

# SARIMAX Residuals
sarimax_resid = sarimax_fit.resid
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].hist(sarimax_resid, bins=20, edgecolor='black')
axes[0].set_title('SARIMAX - Residual Distribution')
stats.probplot(sarimax_resid, plot=axes[1])
axes[1].set_title('SARIMAX - Q-Q Plot')
plot_acf(sarimax_resid, ax=axes[2], lags=40)
axes[2].set_title('SARIMAX - Residual ACF')
plt.tight_layout()
plt.show()

What do these plots show?

The first chart overlays the SARIMAX forecast (purple dashed) on your test data (orange). Compare visually to the SARIMA chart from the previous step to see if the calendar-enhanced model tracks better during holiday periods.

The residual diagnostics (histogram, Q-Q, ACF) tell you whether adding calendar effects reduced the leftover structure. If the ACF shows fewer significant spikes than the pure SARIMA residuals, the calendar features captured real signal.

  • Run the cell.

You should see the SARIMAX forecast line tracking the test data, and three residual diagnostic plots below it.

SARIMAX taking too long to fit?

With seasonal_order=(1,1,1,52), the model estimates many parameters. On slower machines this can take up to 5 minutes. If it stalls beyond that, check your kernel is not stuck and restart it.

my SARIMAX model is taking too long to fit, what can I do?

Fit Prophet with French holidays

Prophet takes a different approach. Instead of manually engineering calendar columns, it has built-in holiday handling. You add French holidays with a single method call, and Prophet learns their effect during fitting.

Why use both SARIMAX and Prophet for calendar effects?

SARIMAX gives you explicit coefficients you can interpret ("each public holiday reduces maintenance by X units"). Prophet gives you a decomposable component plot showing holiday effects visually.

Comparing both tells you whether the signal is real and robust, or just an artifact of one model's assumptions.

  • Create a new cell and fit the Prophet model by pasting this code:
# %% Prophet Model
df_prophet = df.reset_index().rename(columns={df.index.name or 'index': 'ds', 'Maintenance': 'y'})
train_prophet = df_prophet.iloc[:-test_size]
test_prophet = df_prophet.iloc[-test_size:]

prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False
)
prophet_model.add_country_holidays(country_name='FR')
prophet_model.fit(train_prophet)

future_prophet = prophet_model.make_future_dataframe(periods=test_size, freq='W')
forecast_prophet = prophet_model.predict(future_prophet)

prophet_pred = forecast_prophet.iloc[-test_size:]['yhat'].values
results.append(compute_metrics(test['Maintenance'].values, prophet_pred, 'Prophet'))

What does this code do?

  • Prophet requires a DataFrame with columns ds (dates) and y (values). The rename call converts your data to this format.
  • weekly_seasonality=False is disabled because your data is already at weekly granularity. Yearly seasonality captures the annual cycle.
  • add_country_holidays(country_name='FR') loads all French public holidays automatically. No manual feature engineering needed.
  • make_future_dataframe generates the dates Prophet needs for prediction, and predict returns the forecast with confidence intervals.
  • Run the cell. Prophet will print some fitting messages.

You should see Prophet's fitting progress complete without errors. The results list now contains 4 entries.

  • Now visualize Prophet's components and residuals by adding this final cell:
# Prophet components
fig = prophet_model.plot_components(forecast_prophet)
plt.tight_layout()
plt.show()

plt.figure(figsize=(14, 4))
plt.plot(train.index, train['Maintenance'], label='Train')
plt.plot(test.index, test['Maintenance'], label='Test', color='orange')
plt.plot(test.index, prophet_pred, label='Prophet Forecast', color='brown', linestyle='--')
plt.legend()
plt.title('Prophet - Fit vs Actual')
plt.tight_layout()
plt.show()

# Prophet Residuals
prophet_fitted = forecast_prophet.iloc[:len(train_prophet)]['yhat'].values
prophet_resid = train_prophet['y'].values - prophet_fitted
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].hist(prophet_resid, bins=20, edgecolor='black')
axes[0].set_title('Prophet - Residual Distribution')
stats.probplot(prophet_resid, plot=axes[1])
axes[1].set_title('Prophet - Q-Q Plot')
plot_acf(prophet_resid, ax=axes[2], lags=40)
axes[2].set_title('Prophet - Residual ACF')
plt.tight_layout()
plt.show()

What do these plots show?

The component plot breaks Prophet's forecast into its pieces: trend, yearly seasonality, and holiday effects. Look at the holidays panel to see whether Prophet detected meaningful dips on French public holidays.

The fit-vs-actual chart and residual diagnostics let you compare Prophet's quality to the SARIMAX results you saw earlier. If both models show holiday effects matter, you have strong evidence that calendar-awareness improves forecasts.

  • Run the cell.

Before you look at the output, what do you expect the Prophet component plot to show for holiday effects? Do you think holidays will increase or decrease maintenance?

You should see the Prophet component plot with separate panels for trend, yearly seasonality, and holidays. The holidays panel should show downward spikes on major holiday dates. Below that, the fit-vs-actual chart shows Prophet's brown dashed line against your orange test data, followed by three residual diagnostic plots.

Prophet component plot missing the holidays panel?

If the holidays panel is absent, check that add_country_holidays(country_name='FR') was called before .fit(). Prophet only shows components it was configured with before fitting.

my Prophet model is not showing the holidays component, how do I fix it?

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all four new cells: Calendar Features Construction, SARIMAX with Calendar Effects, Prophet Model, and Prophet components/residuals.

ⓧ I'd like to double check the full code

Your notebook should now contain all the cells from the previous steps, plus the four new cells you added in this step. Scroll through and confirm each section heading is present, in order:

  • Calendar Features Construction (the build_calendar_features function definition).
  • The cell calling build_calendar_features for train, test, all, and future indexes.
  • SARIMAX with Calendar Effects (model fit, forecast, and metrics).
  • SARIMAX plot and residual diagnostics.
  • Prophet Model (data prep, fit, forecast, and metrics).
  • Prophet component plot, fit-vs-actual chart, and residual diagnostics.

You now have 4 models evaluated: ETS, SARIMA, SARIMAX with calendar effects, and Prophet with French holidays. The calendar-aware models should show improved metrics, especially around holiday periods. Next up, you'll add XGBoost to see whether a tree-based model can capture non-linear patterns the statistical models miss.

Build the XGBoost Forecasting Model

Your four statistical models all share a fundamental assumption: the relationship between time, seasonality, and maintenance volume is linear. But real-world maintenance data often contains non-linear interactions between features that linear models cannot capture.

In this step, you will train an XGBoost model that uses gradient-boosted decision trees to learn complex patterns. Unlike ARIMA-family models, XGBoost requires you to manually engineer features from the time series. In return, it can discover interactions that statistical models miss entirely.

In this step, get ready to:
  • Engineer time-based, lag, rolling, and calendar features for XGBoost.
  • Train an XGBRegressor and evaluate it on the test set.
  • Analyze residual diagnostics to check for remaining structure.
Create the XGBoost feature engineering function

XGBoost does not understand time natively. You need to encode temporal information as numeric columns: which week of the year it is, what the value was 1 week ago, what the rolling average looks like, and whether holidays are present.

  • In your notebook, create a new cell below the Prophet section and define the feature engineering function by pasting this code:
# %% XGBoost Feature Engineering
def create_xgb_features(dataframe, calendar_df):
    """Create time-based + calendar features for XGBoost."""
    feat = pd.DataFrame(index=dataframe.index)
    feat['week_of_year'] = dataframe.index.isocalendar().week.astype(int)
    feat['month'] = dataframe.index.month
    feat['quarter'] = dataframe.index.quarter
    feat['year'] = dataframe.index.year
    # Lag features
    feat['lag_1'] = dataframe['Maintenance'].shift(1)
    feat['lag_4'] = dataframe['Maintenance'].shift(4)
    feat['lag_52'] = dataframe['Maintenance'].shift(52)
    # Rolling features
    feat['rolling_4'] = dataframe['Maintenance'].rolling(4).mean()
    feat['rolling_12'] = dataframe['Maintenance'].rolling(12).mean()
    # Calendar features
    feat['nb_jours_feries'] = calendar_df['nb_jours_feries']
    feat['pct_vacances'] = calendar_df['pct_vacances']
    return feat

What does this function do?

This function converts raw time series data into a feature matrix XGBoost can learn from. It builds four types of features:

  • Time encodings (week_of_year, month, quarter, year) tell the model where it is in the calendar cycle.
  • Lag features (lag_1, lag_4, lag_52) give the model recent and year-ago values as predictive context.
  • Rolling means (rolling_4, rolling_12) smooth out noise and capture short-term and medium-term trends.
  • Calendar features (nb_jours_feries, pct_vacances) reuse the holiday and school break data from Step 3.
  • Now build the feature matrix and split it into train/test sets by adding this code in the same cell:
xgb_features = create_xgb_features(df, calendar_all)
xgb_features['target'] = df['Maintenance']
xgb_features = xgb_features.dropna()

# Split
xgb_train = xgb_features.iloc[:-test_size]
xgb_test = xgb_features.iloc[-test_size:]

feature_cols = [c for c in xgb_features.columns if c != 'target']
X_train = xgb_train[feature_cols]
y_train = xgb_train['target']
X_test = xgb_test[feature_cols]
y_test = xgb_test['target']

Why drop NaN rows?

The lag and rolling features create NaN values at the beginning of the dataset. For example, lag_52 needs 52 weeks of history before it can produce a value. Calling .dropna() removes those incomplete rows so XGBoost receives only valid data.

  • Run the cell.
  • After the cell finishes, verify the feature matrix shape by running this in a new cell:
print(f"Features: {X_train.shape[1]} columns, Train: {X_train.shape[0]} rows, Test: {X_test.shape[0]} rows")

You should see 11 feature columns, with the train set having fewer rows than your original train DataFrame (because of the dropped NaN rows from lag_52).

Seeing 0 rows or an error?

If X_train has 0 rows, your dataset may be too short for lag_52 to produce any valid values. Check that your data starts from 2023-01-01 and has at least 64 weeks of history (52 for the lag + 12 for the test set).

If you see a KeyError on Maintenance, confirm that your DataFrame column is named exactly Maintenance (case-sensitive).

Help me debug my XGBoost feature matrix

Train XGBoost and inspect feature importance
  • Create a new cell and train the XGBoost model by pasting this code:
# %% XGBoost Training
xgb_model = XGBRegressor(
    n_estimators=200,
    max_depth=5,
    learning_rate=0.1,
    random_state=42
)
xgb_model.fit(X_train, y_train)

xgb_pred = xgb_model.predict(X_test)
results.append(compute_metrics(y_test.values, xgb_pred, 'XGBoost'))

What do these hyperparameters control?

  • n_estimators=200 builds 200 sequential trees, each correcting the errors of the previous ones.
  • max_depth=5 limits each tree to 5 levels deep, preventing overfitting to noise.
  • learning_rate=0.1 shrinks each tree's contribution so the ensemble learns gradually.
  • random_state=42 ensures reproducible results across runs.
  • Now add the feature importance visualization and test predictions plot in the same cell:
# Feature importance
importances = pd.Series(xgb_model.feature_importances_, index=feature_cols).sort_values()
plt.figure(figsize=(10, 6))
importances.plot(kind='barh')
plt.title('XGBoost - Feature Importance')
plt.tight_layout()
plt.show()

plt.figure(figsize=(14, 4))
plt.plot(train.index, train['Maintenance'], label='Train')
plt.plot(test.index, test['Maintenance'], label='Test', color='orange')
plt.plot(xgb_test.index, xgb_pred, label='XGBoost Forecast', color='darkgreen', linestyle='--')
plt.legend()
plt.title('XGBoost - Fit vs Actual')
plt.tight_layout()
plt.show()

What do these plots show?

The feature importance bar chart shows which inputs the model relies on most. Lag features and rolling means typically dominate because recent history is the strongest predictor of near-future maintenance volume.

The test predictions plot overlays the XGBoost forecast on the actual test data. Look for how well it tracks the true values compared to the SARIMA/SARIMAX plots from earlier steps.

  • Run the cell.
  • You should see a horizontal bar chart of feature importances and a time series plot with the XGBoost forecast (dark green dashed line) overlaid on the orange test data.

Feature importance plot is empty or flat?

If all bars are the same height or the plot appears empty, check that X_train has more than a handful of rows. XGBoost needs sufficient training data to differentiate between features.

Help me troubleshoot my XGBoost feature importance

Plot XGBoost residual diagnostics

Residual analysis reveals whether the model left structure in the errors. If the residuals still show patterns, the model is missing signal that a different approach might capture.

  • Add the residual diagnostic plots in a new cell by pasting this code:
# XGBoost Residuals
xgb_resid = y_test.values - xgb_pred
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].hist(xgb_resid, bins=20, edgecolor='black')
axes[0].set_title('XGBoost - Residual Distribution')
stats.probplot(xgb_resid, plot=axes[1])
axes[1].set_title('XGBoost - Q-Q Plot')
plot_acf(pd.Series(xgb_resid), ax=axes[2], lags=min(10, len(xgb_resid) - 1))
axes[2].set_title('XGBoost - Residual ACF')
plt.tight_layout()
plt.show()

How to read the three residual plots

  • Histogram: a roughly bell-shaped distribution centered near zero means the model is not systematically over- or under-predicting.
  • Q-Q plot: points following the diagonal line confirm residuals are approximately normal. Deviations at the tails indicate extreme prediction errors.
  • ACF: spikes above the blue confidence band mean the model missed recurring patterns at that lag. No significant spikes means the model captured the temporal structure well.
  • Run the cell.

Before you look at the output, predict: do you expect XGBoost residuals to show less autocorrelation (fewer ACF spikes) than the ETS and SARIMA residuals from earlier?

You should see a 3-panel figure with the histogram, Q-Q plot, and ACF. Note that the ACF uses fewer lags than earlier models because the test set has only 12 data points. Compare the residual distribution to what you saw for ETS and SARIMA.

Seeing a ValueError about lags?

The lags=min(10, len(xgb_resid) - 1) expression handles this by limiting lags to one fewer than the number of residuals. If you still see an error, make sure you have exactly 12 test observations (matching test_size = 12).

Help me fix the ACF lag error in my XGBoost residuals

Why XGBoost uses fewer ACF lags

The earlier ETS and SARIMA residual ACF used 40 lags because those residuals spanned the entire training set (100+ points). XGBoost residuals here are computed only on the 12-week test set. With just 12 observations, computing 40 lags would be statistically meaningless.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook with all XGBoost cells.

ⓧ I'd like to double check the full code

Scroll up and confirm you have these code blocks in order, each in their own cell or combined into one XGBoost cell:

  • The create_xgb_features function definition with all 11 feature columns.
  • The feature matrix build, .dropna(), and train/test split into X_train, y_train, X_test, y_test.
  • XGBRegressor training with n_estimators=200, max_depth=5, learning_rate=0.1, random_state=42, predictions, and metrics append.
  • Feature importance bar chart and test predictions vs actuals plot.
  • Residual diagnostics: histogram, Q-Q plot, and ACF with lags=min(10, len(xgb_resid) - 1).

You now have all 5 models trained and evaluated. The results list contains MAE, RMSE, and MAPE for ETS, SARIMA, SARIMAX, Prophet, and XGBoost. Next up, you will compare them side-by-side, pick the winners, and generate final forecasts through 2028.

Compare Models and Generate Operational Forecasts

You now have test-set metrics for all 5 models, but individual numbers in isolation tell you nothing about relative performance.

This step puts all models side-by-side, declares a winner, retrains on all available data, and produces the forecast that planning teams actually need: a specific week, a predicted volume, and a confidence interval translated into operational language.

In this step, get ready to:
  • Rank all 5 models by MAE, RMSE, and MAPE in a comparison table.
  • Retrain models on full data and generate forecasts through December 2028 with 90% confidence intervals.
  • Write a business interpretation translating a forecast week into actionable planning.
Build the model comparison table and bar chart

Your results list already holds metrics dictionaries for all 5 models. Turning that into a ranked DataFrame makes the winner immediately visible.

  • Create a new cell in your notebook and paste the following code to build the comparison table and plot metrics side-by-side:
# %% Model Comparison
results_df = pd.DataFrame(results).sort_values('MAPE')
print("\n" + "=" * 60)
print("MODEL COMPARISON - RANKED BY MAPE")
print("=" * 60)
print(results_df.to_string(index=False))
print("=" * 60)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
results_df_plot = results_df.set_index('Model')
results_df_plot['MAE'].plot(kind='bar', ax=axes[0], color='steelblue')
axes[0].set_title('MAE by Model')
axes[0].tick_params(axis='x', rotation=45)
results_df_plot['RMSE'].plot(kind='bar', ax=axes[1], color='coral')
axes[1].set_title('RMSE by Model')
axes[1].tick_params(axis='x', rotation=45)
results_df_plot['MAPE'].plot(kind='bar', ax=axes[2], color='seagreen')
axes[2].set_title('MAPE (%) by Model')
axes[2].tick_params(axis='x', rotation=45)
plt.tight_layout()
plt.show()

What does this code do?

It converts the results list into a DataFrame sorted by MAPE (Mean Absolute Percentage Error). The model at the top of the table produces forecasts closest to actuals in percentage terms.

The three bar charts let you visually compare whether one model dominates across all metrics or whether the ranking shifts depending on which error measure you prioritize.

  • Run the cell.

You should see a printed table with all 5 models ranked by MAPE, followed by a 3-panel bar chart showing MAE, RMSE, and MAPE comparisons.

Why MAPE as the primary ranking?

MAPE expresses error as a percentage of the actual value, making it interpretable regardless of scale. A planning team understands "our forecast is off by 8%" better than "MAE is 47.3 operations."

Seeing an empty table or KeyError?

Make sure all 5 model cells ran successfully before this cell. The results list should contain exactly 5 dictionaries. Run len(results) to confirm.

Help me debug my model comparison table.

Retrain all models on full data and forecast through 2028

During evaluation you held back the last 12 weeks for testing. Now that you know which models perform well, retrain every model on all available data so the forecasts benefit from the most recent observations.

  • Retrain the three statistical models (ETS, SARIMA, SARIMAX) on the full dataset by pasting this into a new cell:
# %% Final Forecasts - Retrain on All Data
print("\nRetraining on all data for final forecasts through", last_forecast)

# ETS Final
ets_final = ETSModel(
    df['Maintenance'], error='add', trend='add', seasonal='add', seasonal_periods=52
).fit(disp=False)
ets_forecast = ets_final.get_forecast(steps=len(future_index))
ets_fc_mean = ets_forecast.predicted_mean
ets_fc_ci = ets_forecast.conf_int(alpha=0.1)

# SARIMA Final
sarima_final = SARIMAX(
    df['Maintenance'], order=(1, 1, 1), seasonal_order=(1, 1, 1, 52),
    enforce_stationarity=False, enforce_invertibility=False
).fit(disp=False)
sarima_forecast = sarima_final.get_forecast(steps=len(future_index))
sarima_fc_mean = sarima_forecast.predicted_mean
sarima_fc_ci = sarima_forecast.conf_int(alpha=0.1)

# SARIMAX Final
sarimax_final = SARIMAX(
    df['Maintenance'], exog=calendar_all, order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 52),
    enforce_stationarity=False, enforce_invertibility=False
).fit(disp=False)
sarimax_forecast = sarimax_final.get_forecast(steps=len(future_index), exog=calendar_future)
sarimax_fc_mean = sarimax_forecast.predicted_mean
sarimax_fc_ci = sarimax_forecast.conf_int(alpha=0.1)

What does this code do?

Each model is re-fit on the entire df (no train/test split). The .get_forecast() method with alpha=0.1 produces 90% confidence intervals. SARIMAX receives calendar_future as exogenous variables so it knows about future holidays.

SARIMA fitting can be slow

Fitting SARIMA and SARIMAX with seasonal_order=(1, 1, 1, 52) on weekly data can take 1-5 minutes. Let the cell finish before proceeding.

  • Run the cell. You should see the "Retraining on all data" confirmation print while the models fit.
  • Now retrain Prophet on the full dataset by pasting this into a new cell:
# Prophet Final
df_prophet_full = df.reset_index().rename(columns={df.index.name or 'index': 'ds', 'Maintenance': 'y'})
prophet_final = Prophet(yearly_seasonality=True, weekly_seasonality=False, daily_seasonality=False)
prophet_final.add_country_holidays(country_name='FR')
prophet_final.fit(df_prophet_full)
future_prophet_full = prophet_final.make_future_dataframe(periods=len(future_index), freq='W')
forecast_prophet_full = prophet_final.predict(future_prophet_full)
prophet_fc_mean = forecast_prophet_full.iloc[-len(future_index):]['yhat'].values
prophet_fc_lower = forecast_prophet_full.iloc[-len(future_index):]['yhat_lower'].values
prophet_fc_upper = forecast_prophet_full.iloc[-len(future_index):]['yhat_upper'].values

What does this code do?

Prophet produces its own uncertainty intervals via yhat_lower and yhat_upper columns. These are extracted for the future portion of the prediction to align with the other models' confidence intervals.

  • Run the cell. Prophet will print fitting logs as it trains on the full dataset.
  • Now generate XGBoost's recursive forecast by pasting this into a new cell:
# XGBoost Final (recursive forecasting)
xgb_final = XGBRegressor(n_estimators=200, max_depth=5, learning_rate=0.1, random_state=42)
xgb_features_full = create_xgb_features(df, calendar_all).dropna()
xgb_final.fit(xgb_features_full[[c for c in xgb_features_full.columns if c != 'target']],
              df.loc[xgb_features_full.index, 'Maintenance'])

# Recursive forecast for XGBoost
xgb_fc_values = []
history = df['Maintenance'].copy()
for i, future_date in enumerate(future_index):
    feat_row = pd.DataFrame(index=[future_date])
    feat_row['week_of_year'] = future_date.isocalendar()[1]
    feat_row['month'] = future_date.month
    feat_row['quarter'] = future_date.quarter
    feat_row['year'] = future_date.year
    feat_row['lag_1'] = history.iloc[-1]
    feat_row['lag_4'] = history.iloc[-4] if len(history) >= 4 else history.iloc[-1]
    feat_row['lag_52'] = history.iloc[-52] if len(history) >= 52 else history.iloc[-1]
    feat_row['rolling_4'] = history.iloc[-4:].mean()
    feat_row['rolling_12'] = history.iloc[-12:].mean()
    feat_row['nb_jours_feries'] = calendar_future.loc[future_date, 'nb_jours_feries']
    feat_row['pct_vacances'] = calendar_future.loc[future_date, 'pct_vacances']
    pred = xgb_final.predict(feat_row[feature_cols])[0]
    xgb_fc_values.append(pred)
    history = pd.concat([history, pd.Series([pred], index=[future_date])])

xgb_fc_mean = np.array(xgb_fc_values)

Why recursive forecasting for XGBoost?

XGBoost has no native multi-step forecast method. It predicts one step at a time, then feeds that prediction back as a lag feature for the next step. This is called recursive (or iterated) forecasting.

The trade-off: prediction errors accumulate over time because each future step depends on the previous (possibly wrong) prediction. This is why XGBoost forecasts tend to drift more than native time series models over long horizons.

  • Run the cell. This loop iterates over every future week through 2028, so it may take 10-30 seconds.

Seeing a KeyError on calendar_future?

Make sure the calendar features cell from Step 3 ran successfully and that calendar_future covers the full forecast horizon. Run print(calendar_future.index[0], calendar_future.index[-1]) to verify the range matches future_index.

Help me fix a KeyError in my XGBoost recursive forecast loop.

  • Plot all 5 model forecasts on one chart by pasting this into a new cell:
# %% Plot All Forecasts
fig, ax = plt.subplots(figsize=(16, 6))
ax.plot(df.index, df['Maintenance'], 'k-', label='Historical Data', linewidth=1)
ax.plot(future_index, ets_fc_mean, label='ETS', alpha=0.8)
ax.fill_between(future_index, ets_fc_ci.iloc[:, 0], ets_fc_ci.iloc[:, 1], alpha=0.1)
ax.plot(future_index, sarima_fc_mean, label='SARIMA', alpha=0.8)
ax.fill_between(future_index, sarima_fc_ci.iloc[:, 0], sarima_fc_ci.iloc[:, 1], alpha=0.1)
ax.plot(future_index, sarimax_fc_mean, label='SARIMAX', alpha=0.8)
ax.fill_between(future_index, sarimax_fc_ci.iloc[:, 0], sarimax_fc_ci.iloc[:, 1], alpha=0.1)
ax.plot(future_index, prophet_fc_mean, label='Prophet', alpha=0.8)
ax.fill_between(future_index, prophet_fc_lower, prophet_fc_upper, alpha=0.1)
ax.plot(future_index, xgb_fc_mean, label='XGBoost', alpha=0.8)
ax.axvline(x=df.index[-1], color='gray', linestyle=':', label='Forecast Start')
ax.set_title('Forecasts Through 2028 - All Models with 90% Confidence Intervals')
ax.set_xlabel('Date')
ax.set_ylabel('Maintenance Operations')
ax.legend(loc='upper left')
plt.tight_layout()
plt.show()

What does this chart show?

Historical data appears as a black line on the left. After the grey vertical "Forecast Start" marker, each model's forecast extends as a coloured line with a shaded 90% confidence band around it.

Models that agree on the future trend reinforce confidence. Where they diverge, the confidence intervals show how uncertain each model is about its own prediction.

  • Run the cell.

You should see a wide chart with historical data on the left and 5 coloured forecast lines extending rightward through 2028, each with a shaded confidence band.

Translate a forecast into business language

A planning team does not read MAPE tables. They need a specific statement: "For this week, expect this volume, with this uncertainty, and here is what to do about it."

  • Compute the forecast values for a target Christmas week by pasting this into a new cell:
# %% Business Interpretation Example
# Pick a Christmas week as example
target_week = pd.Timestamp('2025-12-22')
# Find closest forecast week
closest_idx = np.argmin(np.abs(future_index - target_week))
target_date = future_index[closest_idx]

# Use best model (SARIMAX as example - replace with actual best)
best_forecast_value = sarimax_fc_mean.iloc[closest_idx]
best_ci_low = sarimax_fc_ci.iloc[closest_idx, 0]
best_ci_high = sarimax_fc_ci.iloc[closest_idx, 1]
annual_mean = df['Maintenance'].mean()
pct_vs_mean = ((best_forecast_value - annual_mean) / annual_mean) * 100

What does this code do?

It picks the week of 2025-12-22 (Christmas) as a concrete example. The code finds the closest forecast week in future_index, extracts the SARIMAX point forecast and 90% confidence bounds, and computes how far that forecast deviates from the historical annual mean.

  • Now add the interpretation output by pasting this directly below in the same cell (or in a new cell):
print("\n" + "=" * 70)
print("INTERPRETATION METIER - EXEMPLE DE PREVISION")
print("=" * 70)
print(f"\nSemaine cible : {target_date.strftime('%Y-%m-%d')} (semaine de Noel)")
print(f"Modele utilise : SARIMAX (meilleur MAPE sur validation)")
print(f"\nPrevision : {best_forecast_value:.0f} operations de maintenance")
print(f"Intervalle de confiance 90% : [{best_ci_low:.0f}, {best_ci_high:.0f}]")
print(f"Moyenne annuelle historique : {annual_mean:.0f}")
print(f"Ecart vs moyenne : {pct_vs_mean:+.1f}%")
print(f"\nRecommandation operationnelle :")
if pct_vs_mean < -10:
    print(f"  -> Semaine a activite reduite ({pct_vs_mean:+.1f}% vs moyenne).")
    print(f"  -> Reduire les equipes de techniciens planifiees.")
    print(f"  -> Reporter les maintenances non-urgentes a la semaine suivante.")
elif pct_vs_mean > 10:
    print(f"  -> Semaine a forte activite ({pct_vs_mean:+.1f}% vs moyenne).")
    print(f"  -> Renforcer les equipes de techniciens.")
    print(f"  -> Anticiper les approvisionnements en materiel.")
else:
    print(f"  -> Activite dans la norme ({pct_vs_mean:+.1f}% vs moyenne).")
    print(f"  -> Maintenir le dimensionnement standard des equipes.")
print("=" * 70)

What does this output?

A formatted business interpretation paragraph. It names the target week, the predicted maintenance volume, the 90% confidence interval, how that compares to the annual average, and a concrete operational recommendation (reduce teams, reinforce teams, or maintain standard staffing).

This is the template you would hand to a planning team. They see one number, one range, one recommendation.

  • Before you run this cell, predict: do you expect the Christmas week forecast to be above or below the annual average? Why?
  • Run the cell.

You should see a formatted block starting with "INTERPRETATION METIER" showing the predicted maintenance volume for the Christmas week, its 90% confidence interval, the percentage deviation from the annual mean, and an operational recommendation. Christmas is a holiday period, so expect a negative deviation with a "reduce teams" recommendation.

Seeing a negative forecast or NaN values?

If best_forecast_value is negative or NaN, the SARIMAX model may have diverged during the long-horizon forecast. Check that sarimax_fc_mean contains reasonable values by running sarimax_fc_mean.describe(). If the model is unstable, replace SARIMAX with whichever model ranked first in your comparison table.

Help me fix NaN or negative values in my forecast output.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook file.

ⓧ I'd like to double check the full code

The full notebook is over 30 lines, so confirm each section is present by scrolling through your notebook. You should see these cells in order:

  • Imports and Configuration.
  • Data Loading.
  • Seasonal Decomposition.
  • Train/Test Split + Helper: Metrics.
  • ETS Model.
  • SARIMA Model.
  • Residual Diagnostics (ETS and SARIMA).
  • Calendar Features Construction.
  • SARIMAX with Calendar Effects.
  • Prophet Model.
  • XGBoost Feature Engineering + Training.
  • Model Comparison (table + bar chart).
  • Final Forecasts - Retrain on All Data (ETS, SARIMA, SARIMAX, Prophet, XGBoost recursive).
  • Plot All Forecasts.
  • Business Interpretation Example.

You now have a complete forecasting pipeline: 5 models ranked, operational forecasts through 2028, and a template for translating any future week into a planning decision. This is the deliverable you hand to operational teams.

Secret mission

Ensemble Model with Optimized Weights

No single model dominates across all time periods. Some models nail holiday dips while others capture the long-term trend better. Can you build an ensemble that combines all 5 models with optimized weights to beat every individual model's MAPE?

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your files and packages for future use, or delete them entirely.

Resources you used:

  • Local notebook file: forecast_maintenance.py.
  • Secret mission file: forecast_maintenance_ensemble.py.
  • Python packages installed via pip: statsmodels, prophet, xgboost, jours-feries-france, vacances-scolaires-france, scikit-learn, matplotlib, pandas, numpy.

Keep everything running

No action needed. Choose this if you plan to continue experimenting with your forecasting models or want to extend the analysis with new data.

  • Your notebook files remain on your Desktop, ready to rerun anytime.
  • Your installed Python packages stay available for other projects.

Pause - I'll come back to this later

Shut down any running Jupyter kernel to free up memory. Your files and packages stay installed for when you return.

  • In Jupyter Notebook, click Kernel, then Shut Down to release memory.
  • Close the Jupyter browser tab.
  • In your terminal, press Ctrl+C to stop the Jupyter server if it is still running.

Delete - I don't want to use this again

Remove all project files and uninstall the packages you installed for this project.

  • Shut down any running Jupyter kernel and stop the Jupyter server with Ctrl+C in your terminal.
  • Delete the project files by running:
cd ~/Desktop
rm forecast_maintenance.py
rm forecast_maintenance_ensemble.py
  • Uninstall the forecasting packages by running:
pip uninstall statsmodels prophet xgboost jours-feries-france vacances-scolaires-france -y

Keep common data science packages

Packages like pandas, numpy, matplotlib, and scikit-learn are useful for many Python projects. Only uninstall them if you are certain you will not need them again.

That's a Wrap!

That's a Wrap!

Nice work. You just built a complete forecasting pipeline that benchmarks 5 different models on real weekly maintenance data and produces operational forecasts through 2028.

You've learned how to:

  • Trained and evaluated 5 forecasting models (ETS, SARIMA, SARIMAX with French calendar effects, Prophet, and XGBoost) on weekly Linky concentrator maintenance data, comparing them with MAE, RMSE, and MAPE metrics plus residual diagnostics.
  • Generated operational forecasts through December 2028 with 90% confidence intervals, retrained on all available data.
  • Produced a business interpretation framework that translates statistical forecast outputs into actionable planning recommendations for operational teams.
  • Secret Mission: Built an optimized ensemble model using scipy.optimize to find optimal weights that outperform any individual model on MAPE.

Ready to quiz yourself?