Demand Forecasting With Backtesting

In the Real World · Brief · Machine learning · Advanced · 3 to 5 days

Buying is done four weeks ahead on last year plus a feeling, and stock runs out. Build a better forecast, and prove it is better. Work it yourself before reading the walkthrough.

The situation

Ashcroft Garden sells garden and outdoor goods. You have 156 weeks of weekly units for 12 SKUs, 2023-08-07 to 2026-07-27. Purchase orders go in four weeks before delivery. The buying team wants reorder points at a 90% service level.

The data

The range includes both fast moving lines and lines that sell nothing most weeks. That mix is deliberate.

What the room believes

1. A machine learning model will comfortably beat last year’s numbers.
2. MAPE is a reasonable way to compare forecast accuracy across SKUs.
3. One held-out period is enough to know how good the forecast is.
4. A good point forecast is what the buying team needs.

Definition of done

1. A verdict on each of the four beliefs, with the evidence.
2. A forecast at the horizon the business actually orders at, validated in a way that reflects how it would be used.
3. A comparison against at least one baseline that requires no model at all.
4. A defence of the accuracy metric you chose, including why you rejected the obvious one.
5. Something a buyer can put in a purchase order, and a statement of the service level it delivers.

If you want to go further

• Compute the same accuracy figure on several different test windows and look at the spread before you compare any two models.
• For each SKU, count how many test weeks your chosen metric is able to score at all.
• Work out what proportion of weeks your forecast is above actual demand, and decide whether that is the number the buyer needs.
• Look at what the stockout column does to the history the model is learning from.

In the Real World · Machine learning · Advanced · 3 to 5 days

The buying team orders on last year plus a feeling. A model beats that, by less than anyone hoped. The bigger finding is that the accuracy metric everyone reports ranks the worst forecast in the range as one of the best.

Section 1Problem Definition

No code yet. Forecasting has more ways to fool yourself than any other modelling task on this site, and three of them are decided before a model is fitted.

Business objective

Produce a 4 week ahead forecast per SKU, good enough to set reorder points at a 90% service level, and state honestly how much better it is than what the team does now.

Hypotheses

1. H1. A machine learning model will comfortably beat last year’s numbers.
2. H2. MAPE is a reasonable way to compare forecast accuracy across SKUs.
3. H3. One held-out period is enough to know how good the forecast is.
4. H4. A good point forecast is what the buying team needs.

Section 2Data Collection

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.linear_model import Ridge

SEED, HORIZON, SERVICE = 20260822, 4, 0.90

d = pd.read_csv('data/ashcroft-weekly-demand.csv').sort_values(['sku', 'week_index'])
print('rows  :', len(d))
print('skus  :', d['sku'].nunique(), ' weeks:', d['week_index'].nunique())
print('stockout weeks :', int(d['stockout'].sum()))

rows  : 1,872
skus  : 12  weeks: 156
stockout weeks : 76

Section 3Data Preprocessing

3aDuplicates and schema checks

print('duplicate rows          :', d.duplicated().sum())
print('duplicate sku and week  :', d.duplicated(subset=['sku', 'week_index']).sum())
print('missing cells           :', d.isna().sum().sum())
print('weeks per sku           :', d.groupby('sku').size().unique())

duplicate rows          : 0
duplicate sku and week  : 0
missing cells           : 0
weeks per sku           : [156]

A complete rectangular panel, which is what a lag feature needs. A SKU with missing weeks would silently shift its own history when lagged.

3bHandling categorical mess

Category and SKU are clean. Neither is one-hot encoded into the model: with 12 SKUs, dummies would let the model memorise each line’s mean and call it a forecast. The lag features carry the SKU’s own level instead, which generalises to a new line.

3cDealing with outliers

for sku, g in d.groupby('sku'):
    print(f'{sku}  {g["product_name"].iloc[0]:24s} mean {g["units_sold"].mean():7.1f}  '
          f'sd {g["units_sold"].std():6.1f}  zero weeks {(g["units_sold"] == 0).mean():6.1%}')

AG-1010  Rotary Mower 40cm        mean    48.4  sd   33.2  zero weeks   0.0%
AG-1020  Cordless Hedge Trimmer   mean    29.3  sd   17.3  zero weeks   0.0%
AG-2010  Compost 50L              mean   205.2  sd  102.9  zero weeks   0.0%
AG-2020  Tomato Feed 1L           mean   103.8  sd   77.9  zero weeks   0.0%
AG-2030  Seed Potatoes 2kg        mean    77.9  sd   75.0  zero weeks   1.3%
AG-3010  Patio Cleaner 5L         mean    54.0  sd   33.0  zero weeks   0.0%
AG-3020  Lawn Sand 20kg           mean    41.6  sd   34.5  zero weeks   3.2%
AG-4010  Teak Bench               mean     4.7  sd    5.9  zero weeks  44.2%
AG-4020  Parasol 3m               mean     7.3  sd    9.0  zero weeks  41.7%
AG-4030  Fire Pit Large           mean     3.2  sd    3.8  zero weeks  40.4%
AG-5010  Bird Feeder Deluxe       mean    24.9  sd   10.7  zero weeks   0.0%
AG-5020  Suet Balls 50pk          mean   128.3  sd   72.3  zero weeks   0.0%

Nothing is trimmed. What matters here is not outliers but the split between continuous and intermittent demand: three SKUs sell nothing in around 44.2% of weeks. That is not noise, it is the demand pattern, and section 7b shows what it does to the usual accuracy metric.

3dHandling missing values

None in the file. Plenty are created by the lag features, which is expected: the first 52 weeks of every SKU have no same-week-last-year value. Those rows are dropped from training rather than imputed, because imputing a lag invents history.

3eHandling skewed data

Demand is right-skewed and strongly seasonal. No log transform: the objective is units in the warehouse, and a model optimised on log units systematically under-forecasts the peaks, which are the weeks the buying team actually cares about.

3fData types and normalisation

The features, and the rule that makes them legitimate.

def build(f):
    f = f.copy()
    f['woy'] = pd.to_datetime(f['week_start']).dt.isocalendar().week.astype(int)
    for k, name in [(1, '1'), (2, '2')]:
        f['sin' + name] = np.sin(2 * np.pi * k * f['woy'] / 52.0)
        f['cos' + name] = np.cos(2 * np.pi * k * f['woy'] / 52.0)
    # every lag is at least HORIZON weeks old: a forecast made today for four
    # weeks ahead cannot use sales from next week
    for lag in (HORIZON, HORIZON + 1, HORIZON + 4, 52):
        f['lag_' + str(lag)] = f.groupby('sku')['units_sold'].shift(lag)
    f['roll8'] = (f.groupby('sku')['units_sold'].shift(HORIZON)
                   .rolling(8).mean().reset_index(level=0, drop=True))
    return f


F = build(d)
FEATS = ['sin1', 'cos1', 'sin2', 'cos2', 'on_promotion', 'week_index',
         'lag_4', 'lag_5', 'lag_8', 'lag_52', 'roll8']
print('features :', len(FEATS))
print('rows usable after lags :', F.dropna(subset=FEATS).shape[0], 'of', len(F))

features : 11
rows usable after lags : 1,248 of 1,872

Section 4Exploratory Data Analysis

4aTarget variable analysis

Weekly units, strongly seasonal, with three SKUs whose demand is intermittent rather than continuous. Those three are 25.0% of the range and they decide which accuracy metric is usable.

4bNumerical variables

Seasonality dominates everything. Growing lines peak in spring, furniture in early summer, and wildlife feed runs counter-seasonally into winter. The Fourier terms in 3f exist to give a model that shape without one dummy per week.

F['woy'] = pd.to_datetime(F['week_start']).dt.isocalendar().week.astype(int)
seasonal = (F.groupby(['category', 'woy'])['units_sold'].mean()
             .groupby('category').agg(['idxmax', 'max', 'min']))
for cat, row in seasonal.iterrows():
    print(f'{cat:10s} peaks in week {row["idxmax"][1]:2d}  '
          f'peak {row["max"]:7.1f}  trough {row["min"]:6.1f}  '
          f'ratio {row["max"] / max(row["min"], 0.1):5.1f}x')

4cCategorical variables

Promotions run in 7.4% of SKU weeks and lift demand materially, which is why on_promotion is a feature. It is also the one feature the buyer knows in advance, because the promotional calendar is set months out. That makes it usable at forecast time in a way that no other future variable is.

4dRelationships between variables

Last year’s same week is the single most useful column in the file, which is why the seasonal naive baseline in 7a is hard to beat and why any model that fails to beat it should not be deployed.

usable = F.dropna(subset=['lag_52'])
print('correlation with the same week last year :',
      round(usable['units_sold'].corr(usable['lag_52']), 4))
print('correlation with four weeks ago          :',
      round(usable['units_sold'].corr(usable['lag_4']), 4))

4eTesting our hypotheses

4fSubgroups

The furniture lines are the subgroup that breaks things. They are low volume, high value and intermittent, and they are also the lines where a stockout is most expensive. Every metric decision below is really a decision about how much those three SKUs count.

Section 5Feature Engineering

5aThe leakage trap

Covered in 3f, and worth restating because it is the one that gets published. Any lag shorter than the horizon leaks the future. So does a rolling mean computed without a shift, so does a target encoding of SKU built on the whole file, and so does scaling fitted before the split.

5bNew features

Four Fourier terms for annual seasonality, four lags, and an eight week rolling mean shifted by the horizon. Nothing per-SKU, so a new product can be forecast the day it has enough history rather than after a model rebuild.

5cEncoding

None needed. Promotion is already binary and week of year enters through the Fourier terms rather than as 52 dummies, which would cost more parameters than there are years of data to fit them.

5dFeature selection

Eleven features, fixed before the backtest. Selecting features by backtest performance would overfit the backtest, which is the one thing standing between this project and a number nobody can trust.

Section 6Model Selection

Section 7Model Training

7aBaselines

Three, and the middle one is the one to beat.

7bComparing candidates

def wape(a, p):
    return float(np.sum(np.abs(a - p)) / np.sum(np.abs(a)))


def mape(a, p):
    m = a != 0
    return float(np.mean(np.abs((a[m] - p[m]) / a[m]))) if m.sum() else float('nan')


ORIGINS = [104, 112, 120, 128, 136, 144]


def run_origin(origin):
    """Train on everything before the origin, forecast the next four weeks."""
    tr = F[F['week_index'] < origin].dropna(subset=FEATS)
    te = F[(F['week_index'] >= origin)
           & (F['week_index'] < origin + HORIZON)].dropna(subset=FEATS)
    out = {'actual': te['units_sold'].values, 'sku': te['sku'].values}

    last = F[F['week_index'] == origin - 1].set_index('sku')['units_sold']
    out['naive'] = np.array([last.get(s, 0) for s in te['sku']])

    ly = F[(F['week_index'] >= origin - 52) & (F['week_index'] < origin - 52 + HORIZON)]
    ly_map = {(r.sku, r.week_index + 52): r.units_sold for r in ly.itertuples()}
    out['seasonal_naive'] = np.array([ly_map.get((s, w), 0)
                                      for s, w in zip(te['sku'], te['week_index'])])

    ma = (F[(F['week_index'] >= origin - 8) & (F['week_index'] < origin)]
          .groupby('sku')['units_sold'].mean())
    out['moving_average'] = np.array([ma.get(s, 0) for s in te['sku']])

    out['ridge'] = np.clip(Ridge(alpha=1.0).fit(tr[FEATS], tr['units_sold'])
                           .predict(te[FEATS]), 0, None)
    gbm = HistGradientBoostingRegressor(max_depth=4, learning_rate=0.06, max_iter=350,
                                        random_state=SEED).fit(tr[FEATS], tr['units_sold'])
    out['boosting'] = np.clip(gbm.predict(te[FEATS]), 0, None)

    qm = HistGradientBoostingRegressor(loss='quantile', quantile=SERVICE, max_depth=4,
                                       learning_rate=0.06, max_iter=350,
                                       random_state=SEED).fit(tr[FEATS], tr['units_sold'])
    out['p90'] = np.clip(qm.predict(te[FEATS]), 0, None)
    return out


METHODS = ['naive', 'seasonal_naive', 'moving_average', 'ridge', 'boosting']
folds = [run_origin(o) for o in ORIGINS]

for m in METHODS:
    a = np.concatenate([f['actual'] for f in folds])
    p = np.concatenate([f[m] for f in folds])
    print(f'{m:16s} WAPE {wape(a, p):.4f}  MAPE {mape(a, p):.4f}  '
          f'MAE {np.mean(np.abs(a - p)):6.2f}  bias {np.mean(p - a):+6.2f}')

naive            WAPE 0.3682  MAPE 0.5779  MAE  21.51  bias  -0.93
seasonal_naive   WAPE 0.3456  MAPE 0.3945  MAE  20.19  bias  +1.33
moving_average   WAPE 0.4418  MAPE 1.1966  MAE  25.82  bias  +2.19
ridge            WAPE 0.3374  MAPE 0.8805  MAE  19.72  bias  +0.61
boosting         WAPE 0.3073  MAPE 0.5678  MAE  17.96  bias  -0.10

Boosting is the most accurate at 0.3073 WAPE. Seasonal naive, which is what the team already does, reaches 0.3456.

The model wins by 0.0383 WAPE, an improvement of 11.1% over the baseline. That is real and it is smaller than most forecasting projects promise. H1 is half right.

What MAPE does to the same results

for sku in sorted(d['sku'].unique()):
    a = np.concatenate([f['actual'][f['sku'] == sku] for f in folds])
    p = np.concatenate([f['boosting'][f['sku'] == sku] for f in folds])
    used = int((a != 0).sum())
    print(f'{sku}  WAPE {wape(a, p):.3f}   MAPE {mape(a, p):.3f}   '
          f'weeks MAPE can use {used:2d} of {len(a)}')

AG-1010   WAPE 0.239   MAPE 0.251   weeks MAPE can use 24 of 24
AG-1020   WAPE 0.312   MAPE 0.285   weeks MAPE can use 24 of 24
AG-2010   WAPE 0.293   MAPE 0.477   weeks MAPE can use 24 of 24
AG-2020   WAPE 0.377   MAPE 0.493   weeks MAPE can use 24 of 24
AG-2030   WAPE 0.307   MAPE 0.517   weeks MAPE can use 23 of 24
AG-3010   WAPE 0.196   MAPE 0.291   weeks MAPE can use 24 of 24
AG-3020   WAPE 0.410   MAPE 1.970   weeks MAPE can use 23 of 24
AG-4010   WAPE 1.069   MAPE 0.972   weeks MAPE can use 11 of 24
AG-4020   WAPE 0.811   MAPE 0.661   weeks MAPE can use 15 of 24
AG-4030   WAPE 1.129   MAPE 0.519   weeks MAPE can use 11 of 24
AG-5010   WAPE 0.330   MAPE 0.380   weeks MAPE can use 24 of 24
AG-5020   WAPE 0.218   MAPE 0.281   weeks MAPE can use 24 of 24

The two metrics rank the range differently. Fire Pit Large has a WAPE of 1.1294, among the worst in the range, and a MAPE of 0.5192, among the best.

MASE puts every SKU on the same scale by dividing by what a naive forecast would have managed. 5 of 12 lines beat the benchmark, and the ones that do not are the intermittent lines.

7cHyperparameter tuning

Barely any, and for a reason the backtest makes visible.

fold_wape = [wape(f['actual'], f['boosting']) for f in folds]
for origin, w in zip(ORIGINS, fold_wape):
    print(f'origin week {origin}  WAPE {w:.4f}')
print('spread across folds :', round(max(fold_wape) - min(fold_wape), 4))
print('last origin alone   :', round(fold_wape[-1], 4))
print('average of all six  :', round(wape(np.concatenate([f['actual'] for f in folds]),
                                          np.concatenate([f['boosting'] for f in folds])), 4))

origin week 104  WAPE 0.3645
origin week 112  WAPE 0.2495
origin week 120  WAPE 0.3183
origin week 128  WAPE 0.3111
origin week 136  WAPE 0.3560
origin week 144  WAPE 0.2542
spread across folds : 0.115
last origin alone   : 0.2542
average of all six  : 0.3073

The same model scores between 0.2495 and 0.3645 depending on which four weeks you test it on.

7dFinal evaluation

The buyer does not order the forecast. They order the forecast plus enough cover to hit a service level, and that is a different quantity.

actual = np.concatenate([f['actual'] for f in folds])
point  = np.concatenate([f['boosting'] for f in folds])
p90    = np.concatenate([f['p90'] for f in folds])

print('point forecast covers demand in', round((point >= actual).mean(), 4), 'of weeks')
print('p90 forecast covers demand in  ', round((p90 >= actual).mean(), 4), 'of weeks')
print('extra units carried            ', round((p90 - point).mean(), 2))

point forecast covers demand in 0.5625 of weeks
p90 forecast covers demand in   0.8264 of weeks
extra units carried             18.89

Section 8Documentation and Handoff

What to do, and who owns it

What not to do

• Do not report MAPE. It ignores every week with zero demand, which on the intermittent lines is most of them.
• Do not tune on one holdout window. The folds here differ by more than the model beats the baseline by.
• Do not order to the point forecast. It covers demand in 56.2% of weeks.
• Do not use a lag shorter than the horizon. It produces an excellent backtest and a forecast that cannot be produced on the day it is needed.
• Do not claim the 90% service level is met. The quantile forecast reaches 82.6% on this data, and saying so is the difference between a model people trust and one they stop believing in March.

Reproducibility

What to take from this

• Beat the seasonal naive baseline or do not deploy. Last year’s same week reaches 0.3456 here, and plenty of models never get past it.
• MAPE is not a general purpose metric. It cannot divide by zero, so it quietly drops the weeks it cannot score, and those are the weeks you got wrong.
• One backtest window is one sample. Report the spread across origins, not the best number you found.
• Every lag must be at least as old as the horizon. This is the forecasting version of leaking the target and it is far easier to do by accident.
• Buyers order quantiles. A point forecast is right half the time by construction, which is not a service level.
• Sales are not demand when the shelf was empty. Until that is fixed, the model learns to reproduce your stockouts.

The buying team wanted a forecast that beats last year. They got one, by 11.1%, and two things worth more than that: a metric that does not lie about the intermittent lines, and a reorder point that reflects how uncertain the forecast actually is.