M5 Retail Demand Forecasting

Compare retail demand forecasts and turn errors into inventory decisions.

Introduction

30 Second Summary

Retail teams must decide how much stock to prepare before tomorrow's customers arrive. A poor estimate creates missed sales or inventory left sitting on shelves.

In this project, you will build a reproducible retail demand case study in Google Colab using public M5 sales data. The finished notebook compares a seasonal-naive baseline, SARIMA, and Prophet on the same unseen 28-day period to produce an inventory recommendation.

What You'll Build

When you run the finished notebook, you'll see daily retail demand become a visual dashboard before three forecasts compete for the top spot.

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

  • A visible demand story that shows the changing level, weekly pattern, recent movement, and unusual demand days in one dashboard.
  • A fair model contest where seasonal naive, SARIMA, and Prophet forecast the same unseen 28-day period.
  • An operational recommendation backed by accuracy, asymmetric inventory cost, uncertainty intervals, and production monitoring.
  • Secret Mission: Stress-test whether the champion changes as stockout costs rise from balanced to much more expensive than overstock.

Are there any prerequisites?

You need a Google account with access to Google Colab.

This hard project assumes advanced Python skills plus experience with time-indexed data and out-of-sample forecast evaluation.

Before We Start

This opening commitment anchors the project in a retail replenishment decision before any modelling begins. You are testing whether added forecasting complexity earns its operational cost across a 28-day horizon.

Set Up a Reproducible Colab Runtime

A forecasting comparison is credible when another reviewer can recreate its environment. Managed Google Colab runtimes change over time, which can make the same notebook produce different results.

This step pins Colab to the 2026.07 runtime. It also installs exact library versions before any project imports.

In this step, get ready to:
  • Create a Colab notebook with the CPU-only 2026.07 runtime.
  • Record the project dependencies in requirements.txt.
  • Import the project libraries with verified package versions.
Create and pin the Colab notebook

The runtime determines the operating system plus the Python environment behind your notebook. Pinning it gives every reviewer the same starting point.

  • Go to Google Colab in your browser.
  • Complete Google sign-in if Colab prompts you.
  • Create a blank notebook from the Colab landing page.
  • Confirm that the notebook save status shows Google Drive.

You now have a notebook that Colab saves to Google Drive. The next setting controls the environment that runs its cells.

  • Click Runtime in the top menu.
  • Select Change runtime type.
  • Open the Runtime Version dropdown.
  • Select 2026.07.
  • Confirm that no GPU or TPU accelerator is enabled.
  • Confirm the runtime change.

This runtime switch can sit connecting for two or three minutes. Colab is rebuilding the notebook session during that wait.

Use the value shown in the Runtime Version dropdown to choose your outcome below.

✔️ The pinned runtime is selected

Your notebook is attached to the runtime required by this case study. Keep the hardware accelerator on its CPU-only setting.

ⓧ A different runtime is selected

  • Return to the Runtime menu.
  • Select Change runtime type.
  • Choose 2026.07 from the Runtime Version dropdown.
  • Confirm the runtime change.

ⓧ The pinned runtime is missing

  • Confirm that your Google account is signed into Colab.
  • Refresh the browser page.
  • Return to the Runtime Version dropdown.

Still missing the runtime?

Colab must show the pinned runtime before you continue. Help me find the pinned Colab runtime.

Install the pinned packages

Colab starts each managed runtime with its own virtual machine. A package manifest gives every reviewer the same dependency set.

  • Select the folder icon in Colab's left sidebar.
  • Use the Files pane's new-file control to create requirements.txt.
  • Paste these package pins into requirements.txt:
numpy==2.5.4
pandas==3.0.6
matplotlib==3.11.2
statsmodels==0.15.0
prophet==1.5.0

What does this file capture?

  • The runtime provides Python 3.12.13.
  • The manifest pins NumPy to 2.5.4.
  • The manifest pins pandas to 3.0.6.
  • The manifest pins Matplotlib to 3.11.2.
  • The manifest pins statsmodels to 0.15.0.
  • The manifest pins Prophet to 1.5.0.
  • Save requirements.txt.
  • Return to the Files pane.

Good, the dependency contract now exists in the Colab session. You'll see requirements.txt listed in the Files pane.

The first package installation can take a few minutes while Colab downloads compatible files. Let the cell finish before adding any imports.

  • Add the notebook's first code cell.
  • Install the pinned packages by pasting this command into the setup cell:
%pip install numpy==2.5.4 pandas==3.0.6 matplotlib==3.11.2 statsmodels==0.15.0 prophet==1.5.0

What does this command do?

The %pip install command installs the five exact package versions into the active Colab runtime. Running it before project imports prevents the notebook from loading a different preinstalled version.

  • Run the setup cell with its play button.
  • Accept a runtime restart if Colab requests one after installation.

The setup cell finishes with all five pinned packages installed. That installation gives the project a repeatable library environment.

Package installation failed?

  • Confirm that the notebook still uses the 2026.07 runtime with Python 3.12.13.
  • Rerun the setup cell once if the download was interrupted.
  • Check the final installation lines for the package that stopped the process.
  • Help me debug the Colab package installation.
Import and verify the environment

A successful installation only proves that package files were downloaded. Importing each library proves that the environment can load everything the forecasting case study needs.

  • Add a code cell beneath the setup cell.
  • Import the project libraries by pasting this code:
from pathlib import Path
from urllib.request import urlretrieve
from zipfile import ZipFile
import sys

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import prophet
import statsmodels
from prophet import Prophet
from statsmodels.tsa.statespace.sarimax import SARIMAX

What do these imports prepare?

  • The standard-library imports prepare file paths plus downloads plus ZIP archive access.
  • The data imports prepare numerical arrays plus time-indexed tables.
  • The forecasting imports prepare SARIMA plus Prophet models.
  • The plotting imports prepare the demand dashboard plus forecast comparisons.
  • Run the imports cell with its play button.

You'll see the cell finish successfully without an import exception. Every project library can now load in the pinned runtime.

An import did not load?

  • Run the package installation cell again.
  • Run the imports cell after the installation finishes.
  • Help me fix the failed import.

The imports prove that each library is available. Version dictionaries now turn that successful load into reproducibility evidence.

  • Add another code cell beneath the imports.
  • Define the expected plus actual versions by pasting this code:
EXPECTED_VERSIONS = {
    "python": "3.12.13",
    "numpy": "2.5.4",
    "pandas": "3.0.6",
    "matplotlib": "3.11.2",
    "statsmodels": "0.15.0",
    "prophet": "1.5.0",
}

ACTUAL_VERSIONS = {
    "python": ".".join(map(str, sys.version_info[:3])),
    "numpy": np.__version__,
    "pandas": pd.__version__,
    "matplotlib": matplotlib.__version__,
    "statsmodels": statsmodels.__version__,
    "prophet": prophet.__version__,
}

How are the versions compared?

The EXPECTED_VERSIONS dictionary records the environment this project requires. The ACTUAL_VERSIONS dictionary reads the versions loaded by the current runtime.

  • Run the version-dictionary cell with its play button.

You'll see the cell finish successfully. Both version dictionaries are now defined for the final check.

Can't read the package versions?

  • Run the imports cell before the version-dictionary cell.
  • Confirm that every imported package name matches the version dictionary.
  • Help me repair the version dictionaries.

Before you run the final check, do you expect every loaded version to match its pinned value?

  • Add a final code cell beneath the version dictionaries.
  • Print plus assert every version by pasting this code:
print("Runtime versions:")
for package_name, version in ACTUAL_VERSIONS.items():
    print(f"  {package_name}: {version}")
    assert version == EXPECTED_VERSIONS[package_name], (
        f"Expected {package_name} {EXPECTED_VERSIONS[package_name]}, got {version}"
    )

What does the assertion prove?

The loop prints the loaded version for every runtime component. Each assertion stops execution if an actual version differs from its expected value.

  • Run the assertion cell with its play button.

You'll see Runtime versions: followed by six verified version lines. The cell reaches the end without raising an assertion.

  • Python shows 3.12.13.
  • NumPy shows 2.5.4.
  • pandas shows 3.0.6.
  • Matplotlib shows 3.11.2.
  • statsmodels shows 0.15.0.
  • Prophet shows 1.5.0.

That's the reproducibility proof in place: your notebook now records the exact environment behind every result.

A version assertion stopped the cell?

  • Compare the actual version in the final output with its expected version.
  • Confirm that the installation cell ran before every project import.
  • Restart the runtime if Colab loaded a package before its pinned version was installed.
  • Help me resolve the version mismatch.

✔️ Awesome, I've got everything!

Great, your notebook uses the pinned runtime. Confirm that requirements.txt is saved in the Colab session.

ⓧ I'd like to double check the full code

Your complete requirements.txt file contains these five package pins:

numpy==2.5.4
pandas==3.0.6
matplotlib==3.11.2
statsmodels==0.15.0
prophet==1.5.0

These exact lines preserve the dependency versions used by the notebook. Keep their spelling plus version numbers unchanged.

Your Colab runtime can now reproduce every result in this case study. Next, you'll turn the M5 sales archive into a visible daily demand story.

Build the M5 Demand Dashboard

The reproducible Google Colab runtime from the previous step is ready. This step turns the public M5 files into one daily demand series that you can inspect before forecasting.

A credible forecast must explain the changing demand level. It must also capture weekly structure and unusual days.

In this step, get ready to:
  • Download the M5 archive when needed and aggregate its daily sales columns in chunks.
  • Map the M5 day identifiers to a validated daily date index.
  • Display a four-panel dashboard and write a short interpretation of its demand signals.
Load M5 demand without expanding the hierarchy

The M5 archive stores daily values as wide d_ columns across thousands of item series. Chunked daily demand aggregation sums those columns while keeping memory use controlled.

The loader has two logical halves. Keep both blocks in one Colab code cell so they form a single load_m5_total() function.

  • Add a code cell directly below the version-check cell in the notebook from the previous step.
  • Paste the first loader section into the new cell:
DATA_URL = "https://zenodo.org/records/12636070/files/m5-forecasting-accuracy.zip"
ZIP_PATH = Path("m5-forecasting-accuracy.zip")
HORIZON = 28
WEEK_LENGTH = 7
UNDERAGE_COST = 8.0
OVERAGE_COST = 2.0


def find_archive_member(members, filename):
    matches = [member for member in members if member.endswith(filename)]
    if len(matches) != 1:
        raise RuntimeError(f"Expected one {filename} member, found {matches}")
    return matches[0]


def load_m5_total():
    if not ZIP_PATH.exists():
        print(f"Downloading {DATA_URL}")
        urlretrieve(DATA_URL, str(ZIP_PATH))

    with ZipFile(ZIP_PATH) as archive:
        members = archive.namelist()
        sales_member = find_archive_member(members, "sales_train_evaluation.csv")
        calendar_member = find_archive_member(members, "calendar.csv")

        with archive.open(sales_member) as sales_file:
            header = pd.read_csv(sales_file, nrows=0)
        day_columns = [column for column in header.columns if column.startswith("d_")]
        if not day_columns:
            raise RuntimeError("No M5 day columns beginning with d_ were found")

What does this loader prepare?

  • The constants record the archive location and the assumptions used throughout the case study.
  • The find_archive_member() helper confirms that each required filename appears exactly once inside the archive.
  • The first half of load_m5_total() downloads m5-forecasting-accuracy.zip only when the runtime does not already have it.
  • The header read identifies every daily column while excluding product metadata from the aggregation.
  • Place the cursor on a new line at the end of the same code cell.
  • Append the second loader section:

        total_units = np.zeros(len(day_columns), dtype="float64")
        with archive.open(sales_member) as sales_file:
            chunks = pd.read_csv(sales_file, usecols=day_columns, dtype=np.float32, chunksize=4000)
            for chunk in chunks:
                total_units += chunk.sum(axis=0).to_numpy(dtype="float64")

        with archive.open(calendar_member) as calendar_file:
            calendar = pd.read_csv(calendar_file, usecols=["d", "date"], parse_dates=["date"])

    date_lookup = calendar.set_index("d").loc[day_columns, "date"]
    demand = pd.Series(total_units, index=pd.DatetimeIndex(date_lookup.to_numpy()), name="units").sort_index()
    demand.index.name = "date"
    demand = demand.asfreq("D")
    assert demand.index.is_unique
    assert demand.index.is_monotonic_increasing
    assert not demand.isna().any()
    assert len(demand) == len(day_columns)
    return demand

Why chunk the sales file?

A fixed-size NumPy array stores one total for each day. Each pandas chunk adds its item-level sales into those totals before the next chunk is read.

The calendar lookup replaces M5 day identifiers with dates. The assertions confirm that the resulting series is sorted and unique.

A short loader check gives you evidence that the archive can be downloaded and aggregated before you build the chart.

  • Add a temporary code cell directly below the completed loader cell.
  • Paste this loader check into the temporary cell:
demand = load_m5_total()
print(f"Loaded {len(demand):,} daily observations")
print(f"Date range: {demand.index.min().date()} to {demand.index.max().date()}")

What does this check prove?

The call creates demand as a daily series named units. The two printed lines expose its observation count and date coverage.

Before you run this cell, consider whether aggregation will produce hundreds or thousands of daily observations.

  • Run the loader-check cell using the play button beside it.

The first run downloads the 48.0 MB archive, so the cell may stay active while the transfer finishes. You should then see Loaded 1,941 daily observations followed by the date range.

That is the memory-conscious load complete. Your notebook now has a validated daily series without creating a full long-format table.

Loader cell not completing?

  • Confirm that the package-import cell has run in the current Colab session.
  • Confirm that both loader blocks are in the same code cell with their indentation preserved.
  • Rerun the cell if the runtime disconnected during the archive download.
  • Ask for help with the loader output.
Build the four-panel dashboard

The dashboard combines the full demand history with weekly shape and recent behavior. It also introduces reusable MAE and RMSE functions for later forecast comparisons.

The plotting function is split across two blocks for readability. Append both blocks to one new code cell before running it.

  • Add a code cell directly below the temporary loader-check cell.
  • Paste the first dashboard section into the new cell:
def mae(actual, forecast):
    return float(np.mean(np.abs(np.asarray(actual, dtype="float64") - np.asarray(forecast, dtype="float64"))))


def rmse(actual, forecast):
    return float(np.sqrt(np.mean((np.asarray(actual, dtype="float64") - np.asarray(forecast, dtype="float64")) ** 2)))


def asymmetric_cost(actual, forecast, underage_cost, overage_cost):
    actual_array = np.asarray(actual, dtype="float64")
    forecast_array = np.asarray(forecast, dtype="float64")
    shortage = np.maximum(actual_array - forecast_array, 0.0)
    excess = np.maximum(forecast_array - actual_array, 0.0)
    return float(np.sum(underage_cost * shortage + overage_cost * excess))


def plot_eda(demand):
    rolling_median = demand.rolling(window=28, center=True, min_periods=14).median()
    deviation = (demand - rolling_median).abs()
    anomalies = demand[deviation >= deviation.quantile(0.99)]
    weekday_labels = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]
    weekday_average = demand.groupby(demand.index.dayofweek).mean().reindex(range(7))
    fig, axes = plt.subplots(2, 2, figsize=(15, 9))
    axes[0, 0].plot(demand.index, demand.values, color="#4C78A8", alpha=0.55, label="Daily units")
    axes[0, 0].plot(rolling_median.index, rolling_median.values, color="#F58518", linewidth=2, label="28-day rolling median")
    axes[0, 0].set_title("M5 total daily unit demand")
    axes[0, 0].legend()

What does the first dashboard section do?

  • The scoring functions convert forecast errors into accuracy and inventory-cost measures.
  • The 28-day rolling median reveals movement in the central demand level while reducing daily noise.
  • The 99th-percentile threshold isolates the largest deviations from that rolling median.
  • The first Matplotlib panel overlays daily units with the rolling median.
  • Place the cursor on a new line after axes[0, 0].legend() in the same cell.
  • Append the rest of the dashboard code:
    axes[0, 1].bar(weekday_labels, weekday_average.values, color="#72B7B2")
    axes[0, 1].set_title("Average demand by weekday")
    axes[0, 1].tick_params(axis="x", rotation=35)
    recent = demand.iloc[-180:]
    axes[1, 0].plot(recent.index, recent.values, color="#54A24B")
    axes[1, 0].set_title("Most recent 180 days")
    axes[1, 1].plot(demand.index, demand.values, color="#B8B8B8", alpha=0.55)
    axes[1, 1].scatter(anomalies.index, anomalies.values, color="#E45756", s=24, label="Top 1% deviations")
    axes[1, 1].set_title("Large deviations from rolling median")
    axes[1, 1].legend()
    fig.autofmt_xdate()
    fig.tight_layout()
    plt.show()


demand = load_m5_total()
print(f"Loaded {len(demand):,} daily observations")
print(f"Date range: {demand.index.min().date()} to {demand.index.max().date()}")
plot_eda(demand)

How do the four views work together?

  • The weekday bar chart exposes the repeating seven-day demand shape.
  • The recent-history panel isolates the final 180 observations from the longer series.
  • The anomaly panel highlights the top one percent of rolling-median deviations.
  • The final lines reload the validated series and display the complete dashboard.

Before you run the completed cell, consider which signal will stand out most clearly across the four panels.

  • Run the completed dashboard cell using the play button beside it.

You should see the observation count and date range again. A four-panel dashboard should appear beneath them with total demand, weekday averages, recent history, and highlighted deviations.

You have turned a wide retail hierarchy into a visible demand story. The trend and weekly pattern are now concrete evidence for the models you test next.

Dashboard missing a panel?

  • Confirm that both dashboard blocks are in the same code cell.
  • Check that the second block remains indented until the final plt.show() line.
  • Rerun the loader cell if ZIP_PATH or find_archive_member() is unavailable in the current runtime.
  • Ask for help with the dashboard layout.
  • Delete the earlier temporary loader-check cell that contains only demand = load_m5_total() and the two print statements.
Interpret and verify the demand story

A portfolio chart becomes useful when your interpretation connects what is visible to the forecasting problem. Focus on evidence that a model must reproduce or explain.

  • Add a markdown cell directly beneath the dashboard using the notebook toolbar.
  • Write one sentence describing the direction of the 28-day rolling median.
  • Write one sentence identifying the strongest weekday demand pattern.
  • Write one sentence connecting the latest 180 days or highlighted anomalies to forecast difficulty.

Before the final check, predict whether rerunning the completed dashboard cell will reproduce every panel without an assertion failure.

  • Run the completed dashboard cell once more using the play button beside it.

You should see Loaded 1,941 daily observations and a date range. The four-panel chart should include the total series with its rolling median, weekday averages, the most recent 180 days, and top-one-percent deviations.

✔️ Awesome, I've got everything!

Great work. Your notebook now reproduces the validated M5 demand series and its four-panel dashboard from the pinned environment.

ⓧ I'd like to double check the full code

Compare your cumulative notebook code with this reference. The cells should follow the same order shown here.

from pathlib import Path
from urllib.request import urlretrieve
from zipfile import ZipFile
import sys

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import prophet
import statsmodels
from prophet import Prophet
from statsmodels.tsa.statespace.sarimax import SARIMAX

EXPECTED_VERSIONS = {
    "python": "3.12.13",
    "numpy": "2.5.4",
    "pandas": "3.0.6",
    "matplotlib": "3.11.2",
    "statsmodels": "0.15.0",
    "prophet": "1.5.0",
}

ACTUAL_VERSIONS = {
    "python": ".".join(map(str, sys.version_info[:3])),
    "numpy": np.__version__,
    "pandas": pd.__version__,
    "matplotlib": matplotlib.__version__,
    "statsmodels": statsmodels.__version__,
    "prophet": prophet.__version__,
}

print("Runtime versions:")
for package_name, version in ACTUAL_VERSIONS.items():
    print(f"  {package_name}: {version}")
    assert version == EXPECTED_VERSIONS[package_name], (
        f"Expected {package_name} {EXPECTED_VERSIONS[package_name]}, got {version}"
    )

DATA_URL = "https://zenodo.org/records/12636070/files/m5-forecasting-accuracy.zip"
ZIP_PATH = Path("m5-forecasting-accuracy.zip")
HORIZON = 28
WEEK_LENGTH = 7
UNDERAGE_COST = 8.0
OVERAGE_COST = 2.0


def find_archive_member(members, filename):
    matches = [member for member in members if member.endswith(filename)]
    if len(matches) != 1:
        raise RuntimeError(f"Expected one {filename} member, found {matches}")
    return matches[0]


def load_m5_total():
    if not ZIP_PATH.exists():
        print(f"Downloading {DATA_URL}")
        urlretrieve(DATA_URL, str(ZIP_PATH))

    with ZipFile(ZIP_PATH) as archive:
        members = archive.namelist()
        sales_member = find_archive_member(members, "sales_train_evaluation.csv")
        calendar_member = find_archive_member(members, "calendar.csv")

        with archive.open(sales_member) as sales_file:
            header = pd.read_csv(sales_file, nrows=0)
        day_columns = [column for column in header.columns if column.startswith("d_")]
        if not day_columns:
            raise RuntimeError("No M5 day columns beginning with d_ were found")

        total_units = np.zeros(len(day_columns), dtype="float64")
        with archive.open(sales_member) as sales_file:
            chunks = pd.read_csv(sales_file, usecols=day_columns, dtype=np.float32, chunksize=4000)
            for chunk in chunks:
                total_units += chunk.sum(axis=0).to_numpy(dtype="float64")

        with archive.open(calendar_member) as calendar_file:
            calendar = pd.read_csv(calendar_file, usecols=["d", "date"], parse_dates=["date"])

    date_lookup = calendar.set_index("d").loc[day_columns, "date"]
    demand = pd.Series(total_units, index=pd.DatetimeIndex(date_lookup.to_numpy()), name="units").sort_index()
    demand.index.name = "date"
    demand = demand.asfreq("D")
    assert demand.index.is_unique
    assert demand.index.is_monotonic_increasing
    assert not demand.isna().any()
    assert len(demand) == len(day_columns)
    return demand


def mae(actual, forecast):
    return float(np.mean(np.abs(np.asarray(actual, dtype="float64") - np.asarray(forecast, dtype="float64"))))


def rmse(actual, forecast):
    return float(np.sqrt(np.mean((np.asarray(actual, dtype="float64") - np.asarray(forecast, dtype="float64")) ** 2)))


def asymmetric_cost(actual, forecast, underage_cost, overage_cost):
    actual_array = np.asarray(actual, dtype="float64")
    forecast_array = np.asarray(forecast, dtype="float64")
    shortage = np.maximum(actual_array - forecast_array, 0.0)
    excess = np.maximum(forecast_array - actual_array, 0.0)
    return float(np.sum(underage_cost * shortage + overage_cost * excess))


def plot_eda(demand):
    rolling_median = demand.rolling(window=28, center=True, min_periods=14).median()
    deviation = (demand - rolling_median).abs()
    anomalies = demand[deviation >= deviation.quantile(0.99)]
    weekday_labels = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]
    weekday_average = demand.groupby(demand.index.dayofweek).mean().reindex(range(7))
    fig, axes = plt.subplots(2, 2, figsize=(15, 9))
    axes[0, 0].plot(demand.index, demand.values, color="#4C78A8", alpha=0.55, label="Daily units")
    axes[0, 0].plot(rolling_median.index, rolling_median.values, color="#F58518", linewidth=2, label="28-day rolling median")
    axes[0, 0].set_title("M5 total daily unit demand")
    axes[0, 0].legend()
    axes[0, 1].bar(weekday_labels, weekday_average.values, color="#72B7B2")
    axes[0, 1].set_title("Average demand by weekday")
    axes[0, 1].tick_params(axis="x", rotation=35)
    recent = demand.iloc[-180:]
    axes[1, 0].plot(recent.index, recent.values, color="#54A24B")
    axes[1, 0].set_title("Most recent 180 days")
    axes[1, 1].plot(demand.index, demand.values, color="#B8B8B8", alpha=0.55)
    axes[1, 1].scatter(anomalies.index, anomalies.values, color="#E45756", s=24, label="Top 1% deviations")
    axes[1, 1].set_title("Large deviations from rolling median")
    axes[1, 1].legend()
    fig.autofmt_xdate()
    fig.tight_layout()
    plt.show()


demand = load_m5_total()
print(f"Loaded {len(demand):,} daily observations")
print(f"Date range: {demand.index.min().date()} to {demand.index.max().date()}")
plot_eda(demand)

Your demand dashboard now exposes the benchmark patterns hiding inside the M5 hierarchy. Next, you will reserve an unseen holdout and test how far a repeated week can take you.

Test the Seasonal Baseline

The demand dashboard in Google Colab revealed changing demand levels plus a weekly rhythm. The next decision is whether a simple forecasting rule can turn that rhythm into a credible forecast.

A seasonal-naive forecast repeats the most recent observed week across the next four weeks. This creates a leakage-safe benchmark that every more complex model must earn the right to replace.

In this step, get ready to:
  • Reserve the final 28 days as an untouched holdout.
  • Generate four forecast weeks from the last seven training values.
  • Evaluate the baseline through numeric errors plus a holdout chart.
Build a leakage-safe seasonal baseline

The final 28 observations simulate demand that was unknown when the forecast was created. Keeping those values out of train gives you an honest out-of-sample test.

  • Add a code cell directly beneath your written EDA interpretation.
  • Build the holdout split plus seasonal forecast by running this code:
train = demand.iloc[:-HORIZON]
test = demand.iloc[-HORIZON:]

seasonal_naive_values = np.tile(
    train.iloc[-WEEK_LENGTH:].to_numpy(),
    HORIZON // WEEK_LENGTH,
)
seasonal_naive = pd.Series(
    seasonal_naive_values,
    index=test.index,
    name="Seasonal Naive",
)

baseline_mae = mae(test, seasonal_naive)
baseline_rmse = rmse(test, seasonal_naive)
print(f"Seasonal naive MAE: {baseline_mae:,.2f}")
print(f"Seasonal naive RMSE: {baseline_rmse:,.2f}")

How does the baseline stay leakage-safe?

  • The train series stops before the final HORIZON observations.
  • The test series preserves those final observations for evaluation.
  • NumPy's np.tile() repeats only the last WEEK_LENGTH training values.
  • The seasonal_naive series uses test.index to align each forecast with its actual date.
  • The mae() result measures the average absolute miss in daily units.
  • The rmse() result gives larger misses more influence.
  • Confirm the cell completes without an exception.
  • Read the MAE line beneath the cell.
  • Read the RMSE line beneath the cell.

You should see one numeric value for seasonal-naive MAE plus one numeric value for seasonal-naive RMSE.

Baseline cell not completing?

  • Confirm the earlier setup cells still show completed runs.
  • Rerun the data-loading cell if the runtime was restarted.
  • Rerun this baseline cell after demand is available.
  • Help me debug my seasonal-naive baseline cell.
Plot the repeated week against reality

MAE plus RMSE summarize the forecast gap. A date-aligned chart reveals when the repeated weekly shape follows demand successfully.

Before you run this cell, do you think one repeated seven-day shape will track every turn in the holdout?

  • Create the holdout comparison by running this code:
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(test.index, test.values, color="black", linewidth=2, label="Actual")
ax.plot(seasonal_naive.index, seasonal_naive.values, color="#4C78A8", linestyle="--", label="Seasonal naive")
ax.set_title("Designed obstacle: can an advanced model beat the repeated week?")
ax.legend()
fig.tight_layout()
plt.show()

What does this chart show?

  • Matplotlib places the untouched actual values on the same date axis as the baseline.
  • The black line represents actual holdout demand.
  • The dashed line represents four copies of the final training week.
  • The shared axis exposes each date where the forecast diverges from reality.

You should see the dashed forecast repeat the same seven-point shape four times. The actual line changes in ways that a fixed weekly cycle cannot reproduce.

You have cleared the key benchmark: every later model must now beat a measured leakage-safe forecast.

Chart missing a line?

  • Check that the baseline cell created test plus seasonal_naive.
  • Rerun the baseline cell before rerunning the chart cell.
  • Help me fix my seasonal-naive holdout chart.
  • Add a Markdown cell beneath the chart.
  • Record how the repeated weekly cycle compares with the changing holdout in one or two sentences.

Your notebook should now end with two error values, a 28-day comparison chart, plus your written interpretation of the baseline's limitation.

Next, you will challenge this hurdle with SARIMA plus Prophet on the same holdout.

Challenge the Baseline

Your leakage-safe seasonal-naive forecast now sets a weekly benchmark on the untouched 28-day holdout. Every challenger must face the same unseen demand.

The repeated week cannot adapt when demand shifts away from its last observed level. It also offers no model-based uncertainty.

You will fit SARIMA on the existing training data. You will give Prophet the identical holdout challenge.

Residual diagnostics expose leftover structure. Uncertainty intervals show a plausible range around each point forecast.

In this step, get ready to:
  • Fit SARIMA on the training series and inspect its residual diagnostics.
  • Fit Prophet for the same forecast horizon and inspect its components.
  • Compare both challengers against the actual holdout and seasonal-naive forecast.
Fit SARIMA and inspect residuals

The selected SARIMA specification models short-term dependence after regular differencing. Its seven-day seasonal structure targets the weekly pattern visible in demand.

Residuals hold the variation left after fitting. Their diagnostic grid helps you judge whether the model has missed a systematic pattern.

  • Add a code cell beneath the baseline chart cell in the notebook from earlier.
  • Paste this SARIMA workflow into the new cell:
sarima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 7),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
sarima_result = sarima_model.fit()
sarima_forecast = sarima_result.get_forecast(steps=HORIZON)
sarima_frame = sarima_forecast.summary_frame(alpha=0.05)
sarima_prediction = sarima_frame["mean"].clip(lower=0.0)
sarima_prediction.index = test.index
sarima_lower = sarima_frame["mean_ci_lower"].clip(lower=0.0)
sarima_lower.index = test.index
sarima_upper = sarima_frame["mean_ci_upper"].clip(lower=0.0)
sarima_upper.index = test.index

sarima_result.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

How does this SARIMA workflow work?

  • The regular order uses one autoregressive term, one level of differencing, and one moving-average term.
  • The seasonal order applies the same structure across a seven-day cycle.
  • The forecast frame supplies the mean forecast plus its lower and upper confidence bounds.
  • Clipping each forecast series at zero prevents impossible negative unit demand.

Before you run this cell, do you expect the residuals to look entirely pattern-free? Keep your prediction in mind.

  • Run the SARIMA cell from the notebook toolbar.

You will see a four-panel standardized-residual diagnostic grid. This confirms that the fitted model is ready for residual inspection.

SARIMA cell not completing?

  • Return to the version-check cell from earlier if the model import is unavailable.
  • Confirm that the baseline cells have already defined train, test, and HORIZON.
  • Ask for help with the fitted model cell.
Fit Prophet and inspect components

Prophet expects dates in a ds column. It expects numeric observations in a y column.

This configuration enables weekly structure plus yearly structure. Its component chart separates those patterns from the fitted trend.

  • Add a code cell beneath the SARIMA diagnostic output.
  • Paste this Prophet workflow into the new cell:
prophet_train = train.rename_axis("ds").reset_index(name="y")
prophet_future = pd.DataFrame({"ds": test.index})
prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False,
    interval_width=0.95,
)
prophet_model.fit(prophet_train)
prophet_forecast = prophet_model.predict(prophet_future)
prophet_prediction = pd.Series(
    np.maximum(prophet_forecast["yhat"].to_numpy(), 0.0),
    index=test.index,
    name="Prophet",
)
prophet_lower = pd.Series(
    np.maximum(prophet_forecast["yhat_lower"].to_numpy(), 0.0),
    index=test.index,
)
prophet_upper = pd.Series(
    np.maximum(prophet_forecast["yhat_upper"].to_numpy(), 0.0),
    index=test.index,
)

prophet_model.plot_components(prophet_forecast, figsize=(12, 8))
plt.tight_layout()
plt.show()

What does this Prophet workflow do?

  • The training table converts the existing time series into Prophet's required date and target columns.
  • The future table limits prediction to the same 28 dates used by the untouched test set.
  • The model estimates a trend plus weekly and yearly seasonal components.
  • The three forecast series capture the point estimate plus its 95 percent uncertainty interval.

Before you run this cell, which component do you expect to show the strongest structure? Keep your prediction in mind.

  • Run the Prophet cell from the notebook toolbar.

You will see separate views for trend, weekly seasonality, and yearly seasonality. These plots reveal the structures behind Prophet's holdout forecast.

Prophet components missing?

  • Confirm that the training columns are spelled ds and y.
  • Return to the package-version cell from earlier if model fitting cannot start.
  • Ask for help checking the Prophet input tables.
Compare every holdout forecast

Point forecasts show each model's preferred path. Shaded intervals reveal how much uncertainty surrounds the two advanced forecasts.

Every line uses the untouched test index. That shared timeline makes the visual comparison fair.

  • Add a code cell beneath the Prophet component output.
  • Paste this holdout comparison into the new cell:
fig, axes = plt.subplots(2, 1, figsize=(13, 10))
axes[0].plot(test.index, test.values, color="black", linewidth=2.5, label="Actual")
axes[0].plot(seasonal_naive.index, seasonal_naive.values, linestyle="--", label="Seasonal naive")
axes[0].plot(sarima_prediction.index, sarima_prediction.values, label="SARIMA")
axes[0].plot(prophet_prediction.index, prophet_prediction.values, label="Prophet")
axes[0].fill_between(test.index, sarima_lower.values, sarima_upper.values, alpha=0.12, label="SARIMA 95% interval")
axes[0].fill_between(test.index, prophet_lower.values, prophet_upper.values, alpha=0.12, label="Prophet 95% interval")
axes[0].set_title("28-day M5 holdout forecast comparison")
axes[0].legend(ncol=2)

fig.tight_layout()
plt.show()

How does the comparison stay fair?

  • Each point forecast is aligned to the same 28 dates from test.index.
  • Each shaded region spans the lower and upper limits from one challenger.
  • The seasonal-naive line remains visible as the low-maintenance benchmark.
  • Matplotlib creates a second panel that is reserved for the model scorecard in the next step.

Before you run this cell, do you expect all three point forecasts to follow the same path? Keep your prediction in mind.

  • Run the comparison cell from the notebook toolbar.

In the upper panel, you will see actual demand plus seasonal naive, SARIMA, and Prophet. You will also see shaded 95 percent intervals around both advanced forecasts.

You have now forced both advanced models to answer the same holdout challenge. Their evidence is ready for a business decision.

Missing a forecast line or interval?

  • Confirm that the SARIMA cell completed before the comparison cell.
  • Confirm that the Prophet cell completed before the comparison cell.
  • Ask for help tracing the missing series.

Both challengers now face the same 28-day evidence test. Next, you will convert their errors into a scorecard and choose an operational champion.

Select a Champion and Production Plan

Your notebook now compares seasonal naive with SARIMA and Prophet on the same unseen 28-day period. You have accuracy evidence plus uncertainty and diagnostic evidence for each challenger.

A retail decision also needs to account for asymmetric inventory cost. This step turns the forecast evidence into a champion recommendation with a safe production fallback.

In this step, get ready to:
  • Build a scorecard using MAE, RMSE, and illustrative inventory cost.
  • Select a champion using technical and operational evidence.
  • Write a production plan with monitoring, alerting, and fallback behavior.
Build the model scorecard

A fair scorecard stores every holdout forecast beside the same actual values. This shared table keeps each metric tied to the identical 28-day test period.

  • Add a new code cell beneath the forecast comparison output.
  • Create the holdout predictions table by running this code:
predictions = pd.DataFrame(
    {
        "actual": test,
        "seasonal_naive": seasonal_naive,
        "sarima": sarima_prediction,
        "prophet": prophet_prediction,
    }
)

What does this table contain?

  • The actual column holds the unseen demand values from the holdout.
  • The remaining columns hold predictions from the three models.
  • The shared date index keeps every prediction aligned with its matching actual value.
  • Run the cell.
  • Inspect the live predictions object in the notebook variable viewer.

You should see 28 rows across the four columns named actual, seasonal_naive, sarima, and prophet.

Do the rows or columns look wrong?

  • Confirm that each forecast uses test.index from the same 28-day holdout.
  • Check that the four column names match the code exactly.
  • Ask for help with the alignment if the table contains missing values. Help me debug the index alignment in my predictions dataframe.

MAE measures the typical absolute miss. RMSE puts more weight on larger misses while the illustrative cost function prices shortages more heavily than excess inventory.

  • Add another code cell beneath the predictions table.
  • Calculate the three score types by running this code:
comparison = pd.DataFrame(
    {
        "Model": ["Seasonal Naive", "SARIMA", "Prophet"],
        "MAE": [
            mae(test, seasonal_naive),
            mae(test, sarima_prediction),
            mae(test, prophet_prediction),
        ],
        "RMSE": [
            rmse(test, seasonal_naive),
            rmse(test, sarima_prediction),
            rmse(test, prophet_prediction),
        ],
        "Illustrative Cost": [
            asymmetric_cost(test, seasonal_naive, UNDERAGE_COST, OVERAGE_COST),
            asymmetric_cost(test, sarima_prediction, UNDERAGE_COST, OVERAGE_COST),
            asymmetric_cost(test, prophet_prediction, UNDERAGE_COST, OVERAGE_COST),
        ],
    }
).sort_values("MAE").reset_index(drop=True)

How does the scorecard work?

  • The same mae() and rmse() functions score every model against test.
  • The asymmetric_cost() function charges UNDERAGE_COST for each underforecast unit.
  • The function charges OVERAGE_COST for each overforecast unit.
  • The final operations sort the models from lowest MAE to highest MAE.

The values of UNDERAGE_COST and OVERAGE_COST represent editable scenario assumptions of $8 and $2. They do not describe historical M5 business costs.

  • Run the scorecard cell.
  • Inspect the live comparison object in the notebook variable viewer.

You should see three model rows across the columns Model, MAE, RMSE, and Illustrative Cost.

Seeing missing or infinite scores?

  • Confirm that each prediction series contains 28 numeric values.
  • Check that UNDERAGE_COST is 8.0.
  • Check that OVERAGE_COST is 2.0.
  • Ask for help if one model still produces an invalid score. Help me trace the invalid score.

A portfolio chart lets a reviewer connect the metric ranking to the daily forecast behavior. It also shows whether a lower average error hides a visibly poor day.

  • Add a new code cell beneath the scorecard cell.
  • Create the forecast and MAE comparison figure by running this code:
fig, axes = plt.subplots(2, 1, figsize=(13, 10))
axes[0].plot(test.index, test.values, color="black", linewidth=2.5, label="Actual")
axes[0].plot(seasonal_naive.index, seasonal_naive.values, linestyle="--", label="Seasonal naive")
axes[0].plot(sarima_prediction.index, sarima_prediction.values, label="SARIMA")
axes[0].plot(prophet_prediction.index, prophet_prediction.values, label="Prophet")
axes[0].fill_between(test.index, sarima_lower.values, sarima_upper.values, alpha=0.12, label="SARIMA 95% interval")
axes[0].fill_between(test.index, prophet_lower.values, prophet_upper.values, alpha=0.12, label="Prophet 95% interval")
axes[0].set_title("28-day M5 holdout forecast comparison")
axes[0].legend(ncol=2)

axes[1].bar(comparison["Model"], comparison["MAE"], color=["#4C78A8", "#F58518", "#54A24B"])
axes[1].set_title("Model comparison by MAE")
axes[1].set_ylabel("MAE in daily units")

fig.tight_layout()
plt.show()

What does this figure prove?

  • The upper panel compares all point forecasts against the actual holdout.
  • The shaded bands show the uncertainty intervals from SARIMA and Prophet.
  • The lower panel turns the MAE ranking into a direct visual comparison.
  • Run the chart cell.
  • Compare the shortest MAE bar with the forecast lines in the upper panel.

You should see the 28-day comparison above a three-bar MAE chart. The shortest bar identifies the lowest-MAE model.

Does the figure fail to render?

  • Confirm that the comparison cell ran before this chart cell.
  • Check that both interval pairs use the same test.index as the forecasts.
  • Ask for help if Matplotlib reports mismatched dimensions. Help me fix the chart dimensions.
Choose the champion and write the recommendation

The lowest-MAE model wins the accuracy comparison. The lowest-cost model wins under the current shortage and excess assumptions.

Before you run the next cell, think about whether those two winners will be the same.

  • Add a new code cell beneath the comparison figure.
  • Print the scorecard and both champion names by running this code:
print("\nModel scorecard")
print(comparison.to_string(index=False, float_format=lambda value: f"{value:,.2f}"))

mae_champion = comparison.loc[comparison["MAE"].idxmin(), "Model"]
cost_champion = comparison.loc[comparison["Illustrative Cost"].idxmin(), "Model"]
print(f"\nMAE champion: {mae_champion}")
print(f"Illustrative-cost champion: {cost_champion}")
print(
    "Cost assumptions: underforecast unit = $8, overforecast unit = $2. "
    "These are editable scenario assumptions, not historical M5 costs."
)

How are the champions selected?

  • The printed scorecard preserves the MAE ordering for a quick technical ranking.
  • The idxmin() calls locate the smallest MAE and illustrative cost.
  • The final message keeps the business assumptions visible beside the recommendation.
  • Run the champion cell.
  • Compare the two printed champion names.

You should see the sorted scorecard followed by an MAE champion and an illustrative-cost champion. Their names may agree or differ because the two criteria penalize errors differently.

Do the champion names look inconsistent?

  • Compare each printed champion with the smallest value in its scorecard column.
  • Confirm that comparison contains numeric metric columns.
  • Ask for help if the selected names do not match the minima. Help me verify the champion selection.

The executive recommendation is where the modeling evidence becomes a decision. A credible recommendation acknowledges operational tradeoffs instead of treating a small metric lead as automatically decisive.

  • Add a new markdown cell beneath the champion output.
  • Name the model that wins on MAE.
  • Name the model that minimizes illustrative cost.
  • Compare the selected model with Seasonal Naive.
  • Explain whether the improvement is operationally meaningful.
  • Weigh the selected model's interpretability against its maintenance burden.
  • State your final champion recommendation.

Your recommendation should now connect accuracy, asymmetric cost, interpretability, maintenance burden, and operational significance. That combination gives a reviewer the reasoning behind your chosen model.

Create the production plan and export the evidence

A champion needs rules for refreshes, monitoring, alerts, and failure recovery. These controls keep the forecasting process useful after the backtest ends.

  • Add a new code cell beneath the executive recommendation.
  • Define and print the production-readiness plan by running this code:
production_plan = {
    "Forecast horizon": "28 days",
    "Candidate retraining cadence": "Refresh after each 28-day actuals window",
    "Monitoring": "Track rolling 28-day MAE and signed bias",
    "Alert": "Investigate when rolling MAE exceeds backtest champion MAE by 20%",
    "Fallback": "Use seasonal naive if the selected model fails or inputs are unavailable",
}

print("\nProduction-readiness plan")
for item, recommendation in production_plan.items():
    print(f"  {item}: {recommendation}")

What does the production plan cover?

  • The 28-day horizon stays aligned with the business decision and evaluation window.
  • The refresh cadence waits for another complete 28-day actuals window.
  • Rolling MAE measures error magnitude while signed bias reveals persistent overforecasting or underforecasting.
  • The alert threshold triggers investigation when rolling MAE rises 20 percent above the backtest champion.
  • Seasonal naive provides a maintainable fallback when the selected model or its inputs fail.
  • Run the production-plan cell.
  • Review each printed recommendation against your executive recommendation.

You should see five plan entries covering the horizon, refresh cadence, monitoring, alert threshold, and fallback.

Is a production-plan line missing?

  • Confirm that all five entries sit inside the production_plan dictionary.
  • Check that the closing brace appears before the print loop.
  • Ask for help if the loop skips an entry. Help me inspect the production plan.

The final exports preserve the holdout evidence and scorecard as CSV files. These files also provide the input for the optional cost-ratio stress test.

  • Add one final code cell beneath the production plan.
  • Save both result files to the temporary runtime by running this code:
predictions.to_csv("model_predictions.csv", index=True)
comparison.to_csv("model_scorecard.csv", index=False)
print("\nSaved model_predictions.csv and model_scorecard.csv")

What gets saved?

  • The model_predictions.csv file keeps the date index plus actual and forecast values.
  • The model_scorecard.csv file keeps the MAE-sorted model scorecard without a dataframe index.
  • The files live in the temporary Colab runtime until you download them or the runtime is reset.
  • Run the export cell.
  • Select the folder icon in the Colab sidebar.
  • Confirm that the notebook has saved to Google Drive.

You should see model_predictions.csv and model_scorecard.csv in the temporary runtime files. The cell output should also name both saved files.

Cannot find the exported files?

  • Refresh the temporary runtime file list after running the export cell.
  • Confirm that the export cell completed without an exception.
  • Ask for help if the files still do not appear. Help me locate the exported CSV files.

✔️ Awesome, I've got everything!

Your notebook now contains the scorecard, both champion names, the executive recommendation, the production plan, and the two exported result files.

ⓧ I'd like to double check the full code

Compare your notebook cells with this complete case-study code. The final section begins at the predictions dataframe.

from pathlib import Path
from urllib.request import urlretrieve
from zipfile import ZipFile
import sys

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import prophet
import statsmodels
from prophet import Prophet
from statsmodels.tsa.statespace.sarimax import SARIMAX

EXPECTED_VERSIONS = {
    "python": "3.12.13",
    "numpy": "2.5.4",
    "pandas": "3.0.6",
    "matplotlib": "3.11.2",
    "statsmodels": "0.15.0",
    "prophet": "1.5.0",
}

ACTUAL_VERSIONS = {
    "python": ".".join(map(str, sys.version_info[:3])),
    "numpy": np.__version__,
    "pandas": pd.__version__,
    "matplotlib": matplotlib.__version__,
    "statsmodels": statsmodels.__version__,
    "prophet": prophet.__version__,
}

print("Runtime versions:")
for package_name, version in ACTUAL_VERSIONS.items():
    print(f"  {package_name}: {version}")
    assert version == EXPECTED_VERSIONS[package_name], (
        f"Expected {package_name} {EXPECTED_VERSIONS[package_name]}, got {version}"
    )

DATA_URL = "https://zenodo.org/records/12636070/files/m5-forecasting-accuracy.zip"
ZIP_PATH = Path("m5-forecasting-accuracy.zip")
HORIZON = 28
WEEK_LENGTH = 7
UNDERAGE_COST = 8.0
OVERAGE_COST = 2.0


def find_archive_member(members, filename):
    matches = [member for member in members if member.endswith(filename)]
    if len(matches) != 1:
        raise RuntimeError(f"Expected one {filename} member, found {matches}")
    return matches[0]


def load_m5_total():
    if not ZIP_PATH.exists():
        print(f"Downloading {DATA_URL}")
        urlretrieve(DATA_URL, str(ZIP_PATH))

    with ZipFile(ZIP_PATH) as archive:
        members = archive.namelist()
        sales_member = find_archive_member(members, "sales_train_evaluation.csv")
        calendar_member = find_archive_member(members, "calendar.csv")

        with archive.open(sales_member) as sales_file:
            header = pd.read_csv(sales_file, nrows=0)
        day_columns = [column for column in header.columns if column.startswith("d_")]
        if not day_columns:
            raise RuntimeError("No M5 day columns beginning with d_ were found")

        total_units = np.zeros(len(day_columns), dtype="float64")
        with archive.open(sales_member) as sales_file:
            chunks = pd.read_csv(
                sales_file,
                usecols=day_columns,
                dtype=np.float32,
                chunksize=4000,
            )
            for chunk in chunks:
                total_units += chunk.sum(axis=0).to_numpy(dtype="float64")

        with archive.open(calendar_member) as calendar_file:
            calendar = pd.read_csv(
                calendar_file,
                usecols=["d", "date"],
                parse_dates=["date"],
            )

    date_lookup = calendar.set_index("d").loc[day_columns, "date"]
    demand = pd.Series(
        total_units,
        index=pd.DatetimeIndex(date_lookup.to_numpy()),
        name="units",
    ).sort_index()
    demand.index.name = "date"
    demand = demand.asfreq("D")

    assert demand.index.is_unique
    assert demand.index.is_monotonic_increasing
    assert not demand.isna().any()
    assert len(demand) == len(day_columns)
    return demand


def mae(actual, forecast):
    actual_array = np.asarray(actual, dtype="float64")
    forecast_array = np.asarray(forecast, dtype="float64")
    return float(np.mean(np.abs(actual_array - forecast_array)))


def rmse(actual, forecast):
    actual_array = np.asarray(actual, dtype="float64")
    forecast_array = np.asarray(forecast, dtype="float64")
    return float(np.sqrt(np.mean((actual_array - forecast_array) ** 2)))


def asymmetric_cost(actual, forecast, underage_cost, overage_cost):
    actual_array = np.asarray(actual, dtype="float64")
    forecast_array = np.asarray(forecast, dtype="float64")
    shortage = np.maximum(actual_array - forecast_array, 0.0)
    excess = np.maximum(forecast_array - actual_array, 0.0)
    return float(np.sum(underage_cost * shortage + overage_cost * excess))


def plot_eda(demand):
    rolling_median = demand.rolling(
        window=28,
        center=True,
        min_periods=14,
    ).median()
    deviation = (demand - rolling_median).abs()
    anomaly_threshold = deviation.quantile(0.99)
    anomalies = demand[deviation >= anomaly_threshold]

    weekday_labels = [
        "Monday",
        "Tuesday",
        "Wednesday",
        "Thursday",
        "Friday",
        "Saturday",
        "Sunday",
    ]
    weekday_average = demand.groupby(demand.index.dayofweek).mean().reindex(range(7))

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

    axes[0, 0].plot(demand.index, demand.values, color="#4C78A8", alpha=0.55, label="Daily units")
    axes[0, 0].plot(rolling_median.index, rolling_median.values, color="#F58518", linewidth=2, label="28-day rolling median")
    axes[0, 0].set_title("M5 total daily unit demand")
    axes[0, 0].legend()

    axes[0, 1].bar(weekday_labels, weekday_average.values, color="#72B7B2")
    axes[0, 1].set_title("Average demand by weekday")
    axes[0, 1].tick_params(axis="x", rotation=35)

    recent = demand.iloc[-180:]
    axes[1, 0].plot(recent.index, recent.values, color="#54A24B")
    axes[1, 0].set_title("Most recent 180 days")

    axes[1, 1].plot(demand.index, demand.values, color="#B8B8B8", alpha=0.55)
    axes[1, 1].scatter(anomalies.index, anomalies.values, color="#E45756", s=24, label="Top 1% deviations")
    axes[1, 1].set_title("Large deviations from rolling median")
    axes[1, 1].legend()

    fig.autofmt_xdate()
    fig.tight_layout()
    plt.show()


demand = load_m5_total()
print(f"Loaded {len(demand):,} daily observations")
print(f"Date range: {demand.index.min().date()} to {demand.index.max().date()}")
plot_eda(demand)

train = demand.iloc[:-HORIZON]
test = demand.iloc[-HORIZON:]

seasonal_naive_values = np.tile(
    train.iloc[-WEEK_LENGTH:].to_numpy(),
    HORIZON // WEEK_LENGTH,
)
seasonal_naive = pd.Series(
    seasonal_naive_values,
    index=test.index,
    name="Seasonal Naive",
)

baseline_mae = mae(test, seasonal_naive)
baseline_rmse = rmse(test, seasonal_naive)
print(f"Seasonal naive MAE: {baseline_mae:,.2f}")
print(f"Seasonal naive RMSE: {baseline_rmse:,.2f}")

fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(test.index, test.values, color="black", linewidth=2, label="Actual")
ax.plot(seasonal_naive.index, seasonal_naive.values, color="#4C78A8", linestyle="--", label="Seasonal naive")
ax.set_title("Designed obstacle: can an advanced model beat the repeated week?")
ax.legend()
fig.tight_layout()
plt.show()

sarima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 7),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
sarima_result = sarima_model.fit()
sarima_forecast = sarima_result.get_forecast(steps=HORIZON)
sarima_frame = sarima_forecast.summary_frame(alpha=0.05)
sarima_prediction = sarima_frame["mean"].clip(lower=0.0)
sarima_prediction.index = test.index
sarima_lower = sarima_frame["mean_ci_lower"].clip(lower=0.0)
sarima_lower.index = test.index
sarima_upper = sarima_frame["mean_ci_upper"].clip(lower=0.0)
sarima_upper.index = test.index

sarima_result.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

prophet_train = train.rename_axis("ds").reset_index(name="y")
prophet_future = pd.DataFrame({"ds": test.index})
prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False,
    interval_width=0.95,
)
prophet_model.fit(prophet_train)
prophet_forecast = prophet_model.predict(prophet_future)
prophet_prediction = pd.Series(
    np.maximum(prophet_forecast["yhat"].to_numpy(), 0.0),
    index=test.index,
    name="Prophet",
)
prophet_lower = pd.Series(
    np.maximum(prophet_forecast["yhat_lower"].to_numpy(), 0.0),
    index=test.index,
)
prophet_upper = pd.Series(
    np.maximum(prophet_forecast["yhat_upper"].to_numpy(), 0.0),
    index=test.index,
)

prophet_model.plot_components(prophet_forecast, figsize=(12, 8))
plt.tight_layout()
plt.show()

predictions = pd.DataFrame(
    {
        "actual": test,
        "seasonal_naive": seasonal_naive,
        "sarima": sarima_prediction,
        "prophet": prophet_prediction,
    }
)

comparison = pd.DataFrame(
    {
        "Model": ["Seasonal Naive", "SARIMA", "Prophet"],
        "MAE": [
            mae(test, seasonal_naive),
            mae(test, sarima_prediction),
            mae(test, prophet_prediction),
        ],
        "RMSE": [
            rmse(test, seasonal_naive),
            rmse(test, sarima_prediction),
            rmse(test, prophet_prediction),
        ],
        "Illustrative Cost": [
            asymmetric_cost(test, seasonal_naive, UNDERAGE_COST, OVERAGE_COST),
            asymmetric_cost(test, sarima_prediction, UNDERAGE_COST, OVERAGE_COST),
            asymmetric_cost(test, prophet_prediction, UNDERAGE_COST, OVERAGE_COST),
        ],
    }
).sort_values("MAE").reset_index(drop=True)

fig, axes = plt.subplots(2, 1, figsize=(13, 10))
axes[0].plot(test.index, test.values, color="black", linewidth=2.5, label="Actual")
axes[0].plot(seasonal_naive.index, seasonal_naive.values, linestyle="--", label="Seasonal naive")
axes[0].plot(sarima_prediction.index, sarima_prediction.values, label="SARIMA")
axes[0].plot(prophet_prediction.index, prophet_prediction.values, label="Prophet")
axes[0].fill_between(test.index, sarima_lower.values, sarima_upper.values, alpha=0.12, label="SARIMA 95% interval")
axes[0].fill_between(test.index, prophet_lower.values, prophet_upper.values, alpha=0.12, label="Prophet 95% interval")
axes[0].set_title("28-day M5 holdout forecast comparison")
axes[0].legend(ncol=2)

axes[1].bar(comparison["Model"], comparison["MAE"], color=["#4C78A8", "#F58518", "#54A24B"])
axes[1].set_title("Model comparison by MAE")
axes[1].set_ylabel("MAE in daily units")

fig.tight_layout()
plt.show()

print("\nModel scorecard")
print(comparison.to_string(index=False, float_format=lambda value: f"{value:,.2f}"))

mae_champion = comparison.loc[comparison["MAE"].idxmin(), "Model"]
cost_champion = comparison.loc[comparison["Illustrative Cost"].idxmin(), "Model"]
print(f"\nMAE champion: {mae_champion}")
print(f"Illustrative-cost champion: {cost_champion}")
print(
    "Cost assumptions: underforecast unit = $8, overforecast unit = $2. "
    "These are editable scenario assumptions, not historical M5 costs."
)

production_plan = {
    "Forecast horizon": "28 days",
    "Candidate retraining cadence": "Refresh after each 28-day actuals window",
    "Monitoring": "Track rolling 28-day MAE and signed bias",
    "Alert": "Investigate when rolling MAE exceeds backtest champion MAE by 20%",
    "Fallback": "Use seasonal naive if the selected model fails or inputs are unavailable",
}

print("\nProduction-readiness plan")
for item, recommendation in production_plan.items():
    print(f"  {item}: {recommendation}")

predictions.to_csv("model_predictions.csv", index=True)
comparison.to_csv("model_scorecard.csv", index=False)
print("\nSaved model_predictions.csv and model_scorecard.csv")

Before your final check, think about whether the same model wins both the accuracy and illustrative-cost comparisons.

  • Rerun the new cells from top to bottom.
  • Review the printed scorecard and champion lines.
  • Review the five production-plan entries.
  • Confirm that both result files remain in the temporary runtime.

You should see an MAE-sorted scorecard, both champion names, the editable cost assumptions, the production-readiness plan, and the saved-file confirmation. Your notebook now connects forecast performance to an operational recommendation.

Secret mission

Stress-Test the Decision Under Different Cost Ratios

Your original champion depends on one inventory-cost ratio. Build a three-scenario sensitivity test to reveal whether the recommendation survives balanced, moderate-stockout, and high-stockout assumptions without refitting any model.

Clean Up Your Resources

Clean Up Your Resources

Google Colab's free-of-charge tier covers this project. No paid cloud infrastructure remains active.

Choose how much of the project you want to keep. The tabs below cover each cleanup choice.

Resources you used:

  • The Google Colab notebook stored in Google Drive.
  • The temporary M5 archive named m5-forecasting-accuracy.zip.
  • The reproducibility file named requirements.txt.
  • The cost sensitivity script named cost_stress_test.py.
  • The holdout predictions saved in model_predictions.csv.
  • The model scorecard saved in model_scorecard.csv.

Keep everything running

No action is needed for the notebook. Choose this option if you plan to revisit the forecasting analysis.

  • Keep the notebook in Google Drive.
  • Download model_predictions.csv from Colab before the runtime ends.
  • Download model_scorecard.csv from Colab before the runtime ends.
  • Download cost_stress_test.py from Colab before the runtime ends.
  • Download requirements.txt from Colab before the runtime ends.

These downloads preserve your results outside the temporary runtime.

Pause - I'll come back to this later

This option keeps the notebook in Google Drive. Temporary runtime files can be recreated by rerunning the notebook.

  • Confirm the notebook's latest changes are saved in Google Drive.
  • Download model_predictions.csv if you want to preserve the holdout predictions.
  • Download model_scorecard.csv if you want to preserve the scorecard.
  • Download cost_stress_test.py if you want to preserve the sensitivity analysis.
  • Download requirements.txt if you want a standalone record of the package pins.
  • Disconnect or close the Colab runtime.

Your notebook remains available for the next session. The temporary runtime stops using resources.

Delete - I don't want to use this again

This option removes the saved notebook plus any project files you downloaded. Closing the runtime removes its temporary files.

  • Disconnect or close the Colab runtime to remove its temporary files.

This clears m5-forecasting-accuracy.zip, requirements.txt, cost_stress_test.py, model_predictions.csv, and model_scorecard.csv from the runtime.

  • Delete the project notebook from Google Drive.
  • Delete any downloaded copy of model_predictions.csv from your computer.
  • Delete any downloaded copy of model_scorecard.csv from your computer.
  • Delete any downloaded copy of cost_stress_test.py from your computer.
  • Delete any downloaded copy of requirements.txt from your computer.
  • Search Google Drive for the project notebook.
  • Check your computer's download folder for the project files.

You should no longer see the notebook in Google Drive. You should no longer see downloaded copies of the project files.

That's the cleanup complete. No project infrastructure remains running.

Nice Work!

Nice Work!

You did it! You completed a reproducible M5 retail demand forecasting case study in Google Colab that turns public retail data into an evidence-backed model recommendation.

You've learned how to:

  • Aggregate the M5 retail hierarchy into a memory-conscious daily demand series. Your four-panel demand dashboard reveals long-term movement, weekly seasonality, recent behavior, plus unusual demand days.
  • Establish a leakage-safe seasonal-naive baseline on an unseen 28-day holdout. You challenged it with SARIMA plus Prophet using MAE plus RMSE. Residual diagnostics plus uncertainty intervals made the evidence inspectable.
  • Translate forecast errors into an editable asymmetric inventory-cost scorecard. You selected a champion from the evidence. Your production plan covers retraining, monitoring, alerting, plus seasonal-naive fallback behavior.
  • Secret Mission: Stress-test the original recommendation with a three-scenario cost-sensitivity matrix. You identified whether increasing stockout penalties changes the champion.

Ready to quiz yourself?