Predict-then-Optimize for Inventory — a worked example on real retail data

Machine Learning in Operations Research · demand distribution → inventory decision

Most forecasting write-ups stop at a prediction. But a forecast is not the goal — a decision is: how many units of each product should we stock? This notebook walks the full predict-then-optimize loop on a real, named-product retail dataset, and measures the one thing that matters — the realized cost of the decisions, in money.

The thesis. When demand is uncertain and intermittent (lots of zero-sales days), a single point forecast is the wrong object. The inventory decision needs a distribution. We will:

1. Explore the data and show demand is overwhelmingly intermittent (Sections 1–3);
2. Engineer leak-free features;
3. Forecast the distribution with quantile gradient boosting;
4. Decide order quantities with the newsvendor rule + a capacity linear program, using per-product economics derived from real prices;
5. Backtest over rolling origins and put a number on the value of modeling uncertainty.

This is the self-contained, research-style companion to a fully engineered version of the same pipeline (a tested, deployable Python package). Same idea, told as one narrative here.

Dataset — UCI Online Retail II. ~1.07M transactions from a UK-based online gift retailer, Dec 2009 – Dec 2011, with real product descriptions and unit prices (openly licensed, CC BY 4.0).

import re
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

pd.set_option("display.width", 120)
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#a0aec0","#6b46c1"

# resolve the cached data whether run from the folder or the repo root
here = Path.cwd()
DATA = next(p for p in [here/"data", here.parent/"data",
                        here/"Z-ML Operations Research-Predict-then-Optimize"/"data"] if (p/"online_retail_II.parquet").exists())
raw = pd.read_parquet(DATA/"online_retail_II.parquet")
print(f"{len(raw):,} transactions | {raw.InvoiceDate.min().date()} -> {raw.InvoiceDate.max().date()}")
raw.head()

1,067,371 transactions | 2009-12-01 -> 2011-12-09

	Invoice	StockCode	Description	Quantity	InvoiceDate	Price	Customer_ID	Country
0	489434	85048	15CM CHRISTMAS GLASS BALL 20 LIGHTS	12	2009-12-01 07:45:00	6.95	13085.0	United Kingdom
1	489434	79323P	PINK CHERRY LIGHTS	12	2009-12-01 07:45:00	6.75	13085.0	United Kingdom
2	489434	79323W	WHITE CHERRY LIGHTS	12	2009-12-01 07:45:00	6.75	13085.0	United Kingdom
3	489434	22041	RECORD FRAME 7" SINGLE SIZE	48	2009-12-01 07:45:00	2.10	13085.0	United Kingdom
4	489434	21232	STRAWBERRY CERAMIC TRINKET BOX	24	2009-12-01 07:45:00	1.25	13085.0	United Kingdom

1 · The raw data, warts and all

One row per line item on an invoice. Before modelling we have to reckon with what real transactional data contains — none of this is error, it is the texture of a live retail system:

column	meaning
Invoice	invoice number; a leading C marks a cancellation/return
StockCode	product code — mostly 5-digit (e.g. 85123A), but some are service codes (POST, DOT, M, BANK CHARGES, AMAZONFEE, ADJUST)
Description	free-text product name
Quantity	units on the line; negative on returns
InvoiceDate	timestamp
Price	unit price (GBP); occasionally ≤ 0 for adjustments
Customer_ID	customer (missing on ~some rows)
Country	ship-to country; ~92% United Kingdom

Known quirks we will handle explicitly: returns (negative quantity), non-product service codes, non-positive prices, extreme bulk/adjustment quantities, and the fact that the firm does not process orders on Saturdays (so daily demand has a structural weekly zero).

# quantify the quirks
print("cancellations (Invoice starts 'C'):", f"{raw.Invoice.str.startswith('C').mean()*100:.1f}%")
print("negative quantity rows            :", f"{(raw.Quantity < 0).mean()*100:.1f}%")
print("non-positive price rows           :", f"{(raw.Price <= 0).mean()*100:.1f}%")
print("United Kingdom share              :", f"{(raw.Country=='United Kingdom').mean()*100:.1f}%")
svc = raw.loc[~raw.StockCode.str.match(r'^\d{5}[A-Za-z]?$'), 'StockCode'].value_counts().head(8)
print("\nexamples of non-product StockCodes:\n", svc.to_string())
print("\nweekday of invoices (0=Mon..6=Sun):",
      raw.InvoiceDate.dt.dayofweek.value_counts().sort_index().to_dict())

cancellations (Invoice starts 'C'): 1.8%
negative quantity rows            : 2.2%
non-positive price rows           : 0.6%
United Kingdom share              : 91.9%

examples of non-product StockCodes:
 StockCode
POST       2122
DOT        1446
M          1421
15056BL     923
C2          282
79323LP     232
D           177
79323GR     123

weekday of invoices (0=Mon..6=Sun): {0: 189084, 1: 196626, 2: 185051, 3: 203149, 4: 153803, 5: 402, 6: 139256}

Note the weekday counts: Saturday (5) is essentially absent — the Saturday closure is real. Now the cleaning.

2 · Cleaning → a demand panel

We reduce the transaction log to a daily demand panel: for each product and calendar day, how many units were demanded. Steps, each reported so nothing is silently dropped:

1. keep United Kingdom (one clean market);
2. drop cancellations and non-positive quantity (we model demand, i.e. gross sales);
3. keep only real product codes (^\d{5}[A-Z]?$), removing postage/fees/adjustments;
4. drop non-positive price;
5. aggregate to (product × day) demand, take the top ~200 products by volume (enough history and signal), and reindex to every calendar day so days with no sales become explicit zeros.

def clean_transactions(df):
    steps = [("raw", len(df))]
    df = df[df.Country == "United Kingdom"].copy();           steps.append(("UK only", len(df)))
    df = df[~df.Invoice.str.startswith("C")];                 steps.append(("drop cancellations", len(df)))
    df = df[df.Quantity > 0];                                 steps.append(("qty > 0", len(df)))
    df = df[df.StockCode.str.match(r"^\d{5}[A-Za-z]?$")];     steps.append(("real product codes", len(df)))
    df = df[df.Price > 0];                                    steps.append(("price > 0", len(df)))
    return df, pd.DataFrame(steps, columns=["step", "rows"])

clean, log = clean_transactions(raw)
log["kept_%"] = (log["rows"] / log["rows"].iloc[0] * 100).round(1)
print(log.to_string(index=False))
print(f"\nremaining: {clean.StockCode.nunique():,} products over "
      f"{clean.InvoiceDate.dt.normalize().nunique()} active days")

              step    rows  kept_%
               raw 1067371   100.0
           UK only  981330    91.9
drop cancellations  964680    90.4
           qty > 0  961223    90.1
real product codes  957393    89.7
         price > 0  954729    89.4

remaining: 4,862 products over 604 active days

N_PRODUCTS = 200
clean["date"] = clean["InvoiceDate"].dt.normalize()
daily = clean.groupby(["StockCode", "date"], as_index=False)["Quantity"].sum()

top = daily.groupby("StockCode")["Quantity"].sum().sort_values(ascending=False).head(N_PRODUCTS).index
daily = daily[daily.StockCode.isin(top)]

all_days = pd.date_range(clean.date.min(), clean.date.max(), freq="D")
panel = (daily.set_index(["StockCode", "date"])["Quantity"]
         .reindex(pd.MultiIndex.from_product([top, all_days], names=["StockCode", "date"]), fill_value=0)
         .rename("demand").reset_index())

desc = clean.groupby("StockCode")["Description"].agg(lambda s: s.mode().iloc[0])
price = clean.groupby("StockCode")["Price"].median()
panel["description"] = panel.StockCode.map(desc)
panel["price"] = panel.StockCode.map(price)

print(f"panel: {len(panel):,} rows | {panel.StockCode.nunique()} products x {panel.date.nunique()} days")
print(f"zero-demand share: {(panel.demand==0).mean():.0%} | median unit price GBP {panel.price.median():.2f}")
panel.head()

panel: 147,800 rows | 200 products x 739 days
zero-demand share: 48% | median unit price GBP 1.25

	StockCode	date	demand	description	price
0	84077	2009-12-01	48	WORLD WAR 2 GLIDERS ASSTD DESIGNS	0.29
1	84077	2009-12-02	435	WORLD WAR 2 GLIDERS ASSTD DESIGNS	0.29
2	84077	2009-12-03	56	WORLD WAR 2 GLIDERS ASSTD DESIGNS	0.29
3	84077	2009-12-04	48	WORLD WAR 2 GLIDERS ASSTD DESIGNS	0.29
4	84077	2009-12-05	48	WORLD WAR 2 GLIDERS ASSTD DESIGNS	0.29

3 · Demand character — the Syntetos–Boylan lens

Before modelling we classify each product's demand pattern, because the pattern dictates what kind of forecast is even sensible. The framework comes from the intermittent-demand forecasting literature.

The model behind it — Croston's method

Classical forecasts (simple exponential smoothing, ARIMA) assume demand occurs every period, so they smear a slow mover's demand into a meaningless fractional average. Croston (1972) fixed this by splitting demand into two separate streams and smoothing each on its own:

• the size $z_t$ of a sale when one occurs, and
• the interval $p_t$ between successive sales.

Each is updated by exponential smoothing (only on periods where demand occurs), and the per-period demand-rate forecast is their ratio $\hat z_t / \hat p_t$. Syntetos & Boylan (2005) later showed Croston's estimator is biased and proposed the bias-corrected SBA (Syntetos–Boylan Approximation — multiply by $1-\alpha/2$, where $\alpha$ is the interval smoothing constant). Croston and SBA are the workhorses for intermittent demand, and the basis for classifying it.

Two orthogonal statistics

Each product is summarized by two numbers — one for sparsity in time, one for variability in quantity:

• ADI — average demand interval $=\dfrac{\#\,\text{periods}}{\#\,\text{nonzero-demand periods}}$. How often it sells (ADI = 1 → sells every day; ADI = 4 → on average one sale every 4 days).
• CV² — squared coefficient of variation of the nonzero sale sizes, $(\sigma/\mu)^2$. How variable the size is when a sale happens (0 → every sale is identical).

The 2×2 map — and where 1.32 and 0.49 come from

	CV² < 0.49 (steady size)	CV² ≥ 0.49 (variable size)
ADI < 1.32 (frequent)	smooth — plain SES is fine	erratic — sells often, wild sizes
ADI ≥ 1.32 (sporadic)	intermittent — Croston / SBA	lumpy — sparse and wild (hardest)

The two cut-offs are not arbitrary. Syntetos, Boylan & Croston (2005) wrote down the theoretical mean-squared error of the competing estimators (SES vs Croston vs SBA) as functions of ADI and CV², and solved for the points where the best estimator switches:

• ADI = 1.32 is where the SBA estimator's accuracy overtakes the alternatives — i.e. where demand becomes intermittent enough that the bias-corrected intermittent method pays off. Below it, sales are frequent and ordinary smoothing is competitive.
• CV² = 0.49 is where the size-variance term begins to dominate the error, separating "predictable size" from "erratic size."

Both numbers fall straight out of equating those MSE expressions — they are approximations (later work, e.g. Kostenko & Hyndman 2006, proposed refined boundaries), so treat them as well-grounded rules of thumb, not physical constants. We now compute ADI and CV² for every product, classify, and plot the map.

def sb_stats(df):
    def one(g):
        y = g["demand"].to_numpy(); nz = y[y > 0]
        adi = len(y) / max(len(nz), 1)
        cv2 = (nz.std() / nz.mean())**2 if len(nz) > 1 else 0.0
        return pd.Series({"mean": y.mean(), "zero_share": (y == 0).mean(), "adi": adi, "cv2": cv2})
    return df.groupby("StockCode").apply(one, include_groups=False)

def sb_class(adi, cv2):
    if adi < 1.32 and cv2 < 0.49: return "smooth"
    if adi >= 1.32 and cv2 < 0.49: return "intermittent"
    if adi < 1.32 and cv2 >= 0.49: return "erratic"
    return "lumpy"

ps = sb_stats(panel)
ps["class"] = [sb_class(a, c) for a, c in zip(ps.adi, ps.cv2)]
print("top-200 panel, share of products by class (%):")
print((ps["class"].value_counts(normalize=True)*100).round(0).to_string())

# The same lens over the WHOLE cleaned catalogue, not just the top 200. ADI and CV-squared have
# closed forms on a zero-filled panel -- every series is exactly len(all_days) long -- so the
# counts need no 4,862 x 739 reindex.
_all = (clean.groupby(["StockCode", "date"], as_index=False)["Quantity"].sum()
        .groupby("StockCode")["Quantity"])
_pos = _all.size()
_adi = len(all_days) / _pos.clip(lower=1)
_cv2 = ((_all.std(ddof=0) / _all.mean())**2).fillna(0.0).where(_pos > 1, 0.0)
_cls = pd.Series([sb_class(a, c) for a, c in zip(_adi, _cv2)], index=_adi.index)
_vc = _cls.value_counts().reindex(["smooth", "intermittent", "erratic", "lumpy"]).fillna(0).astype(int)
print(f"\nfull UK catalogue, {len(_cls):,} products, count by class:")
print(_vc.to_string())

cls_col = {"smooth": GREEN, "intermittent": BLUE, "erratic": ORANGE, "lumpy": RED}
fig, ax = plt.subplots(2, 2, figsize=(13, 9))
ax[0,0].hist(panel.demand.clip(upper=panel.demand.quantile(.995)), bins=60, color=BLUE)
ax[0,0].set_yscale("log"); ax[0,0].set_title("Daily demand is zero-inflated & right-skewed")
ax[0,0].set_xlabel("units/day"); ax[0,0].set_ylabel("count (log)")
ax[0,1].hist(ps.zero_share, bins=30, color=PURP)
ax[0,1].set_title("Share of zero-demand days per product"); ax[0,1].set_xlabel("zero-share")
for cl, c in cls_col.items():
    m = ps["class"] == cl
    ax[1,0].scatter(ps.adi[m], ps.cv2[m], s=14, color=c, alpha=.6, label=cl)
ax[1,0].axvline(1.32, color=GREY, ls="--"); ax[1,0].axhline(0.49, color=GREY, ls="--")
ax[1,0].set_xlim(1, min(ps.adi.quantile(.98), 15)); ax[1,0].set_ylim(0, min(ps.cv2.quantile(.98), 4))
ax[1,0].set_xlabel("ADI"); ax[1,0].set_ylabel("CV² of nonzero demand")
ax[1,0].set_title("Syntetos–Boylan classification"); ax[1,0].legend(fontsize=8)
vc = ps["class"].value_counts()
ax[1,1].bar(vc.index, vc.values, color=[cls_col[k] for k in vc.index])
ax[1,1].set_title("Products by demand class"); ax[1,1].set_ylabel("# products")
fig.suptitle("Demand character of the top-200 UK products", fontsize=13)
fig.tight_layout(); plt.show()

top-200 panel, share of products by class (%):
class
lumpy           94.0
erratic          6.0
intermittent     0.0

full UK catalogue, 4,862 products, count by class:
smooth             0
intermittent     835
erratic           12
lumpy           4015

Even among the best-selling 200 products, almost all land in the lumpy quadrant — Syntetos–Boylan's hardest class: they sell on many days but in wildly variable quantities (a few units one day, a bulk order of hundreds the next). That combination of frequent-but-erratic demand is precisely where a point forecast fails and a full distribution is required — the empirical motivation for everything that follows.

A gallery of demand patterns — named products

The classes are abstract until you see them. Across the full UK catalogue (not just the top-200) the counts are themselves telling — 0 smooth, 835 intermittent, 12 erratic, 4,015 lumpy across all 4,862 products: a gift retailer has essentially no product that sells steadily every day, and four in five are the hardest class of all. Here is one real product from each of the three classes that actually occur:

• Intermittent — WRAP RED APPLES (gift wrap): bought in occasional but fairly consistent batches.
• Erratic — WHITE HANGING HEART T-LIGHT HOLDER (85123A, the shop's signature item): sells almost every open day, but the quantity swings wildly.
• Lumpy — MEDIUM CERAMIC TOP STORAGE JAR: long dead spells punctuated by occasional huge wholesale orders — the hardest to forecast.

gallery = [("22704", "intermittent"), ("85123A", "erratic"), ("23166", "lumpy")]
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
for (sc, label), a in zip(gallery, ax):
    s = clean[clean.StockCode == sc].groupby("date")["Quantity"].sum().reindex(all_days, fill_value=0)
    nz = s[s > 0]
    adi = len(s) / max(len(nz), 1); cv2 = (nz.std() / nz.mean())**2
    a.plot(s.index, s.values, color=GREY, lw=.5)
    a.plot(s.index, s.rolling(14).mean(), color=cls_col[label], lw=1.9)
    a.set_title(f"{label.upper()} · {sc}\n{desc[sc][:30].strip()}\nADI {adi:.2f}, CV² {cv2:.2f}", fontsize=9)
    a.set_ylabel("units/day")
fig.suptitle("Three products, three demand patterns (raw + 14-day mean)", fontsize=13)
fig.tight_layout(); plt.show()

Same catalogue, three very different forecasting problems. The erratic signature item and the lumpy long-tail — which together are the bulk of the business — are exactly where a point forecast collapses and the full demand distribution earns its keep.

Temporal structure and a named example

Retail demand has strong, learnable structure: a Christmas surge, a weekly rhythm with the Saturday closure, and product-level idiosyncrasy. Here is the aggregate, and one recognizable product up close.

fig, ax = plt.subplots(2, 2, figsize=(13, 8))
tot = panel.groupby("date")["demand"].sum()
ax[0,0].plot(tot.index, tot.values, color=GREY, lw=.5)
ax[0,0].plot(tot.index, tot.rolling(28).mean(), color=RED, lw=2)
ax[0,0].set_title("Total daily demand (28-day mean): Christmas surge"); ax[0,0].set_ylabel("units/day")
dow = panel.assign(dow=panel.date.dt.dayofweek).groupby("dow")["demand"].mean()
names = ["Mon","Tue","Wed","Thu","Fri","Sat","Sun"]
ax[0,1].bar(range(7), dow.values, color=BLUE); ax[0,1].set_xticks(range(7)); ax[0,1].set_xticklabels(names)
ax[0,1].set_title("Day-of-week (note the Saturday closure)"); ax[0,1].set_ylabel("mean units")
mo = panel.assign(m=panel.date.dt.month).groupby("m")["demand"].mean()
ax[1,0].plot(mo.index, mo.values, "o-", color=GREEN); ax[1,0].set_xticks(range(1,13))
ax[1,0].set_title("Month-of-year seasonality"); ax[1,0].set_xlabel("month"); ax[1,0].set_ylabel("mean units")
star = ps["mean"].idxmax()
one = panel[panel.StockCode == star].set_index("date")["demand"]
ax[1,1].plot(one.index, one.rolling(7).mean(), color=PURP, lw=1.3)
ax[1,1].set_title(f"{star} — {desc[star][:34]}"); ax[1,1].set_ylabel("units/day (7d mean)")
fig.suptitle("Temporal structure and demand drivers", fontsize=13)
fig.tight_layout(); plt.show()
print("busiest product:", star, "|", desc[star])

busiest product: 84077 | WORLD WAR 2 GLIDERS ASSTD DESIGNS

Section 3 takeaways

After honest cleaning, the top-200 UK products form a daily demand panel that is overwhelmingly lumpy (~48% zero-days and highly variable sale sizes), with a strong Christmas peak and a hard Saturday zero. The distribution — not the mean — is the object to forecast. Cached artifacts for the next sections: the panel (product × day demand + price) and the per-product Syntetos–Boylan table ps.

Next — Section 4: features. We build a leak-free model matrix (calendar, lags, rolling windows, price, product identity), taking care that no feature can peek past the forecast origin.

panel.to_parquet(DATA/"panel.parquet", index=False)
ps.to_parquet(DATA/"sb_stats.parquet")
print("saved panel.parquet and sb_stats.parquet ->", DATA)

saved panel.parquet and sb_stats.parquet -> C:\Users\user\dev\personal-website\Z-ML Operations Research-Predict-then-Optimize\data

4 · Feature engineering — a leak-free model matrix

We will forecast with a gradient-boosted model, which needs a tabular feature matrix: one row per product-day, with columns describing what we know at decision time. The single most important discipline is no leakage — a feature for a target day may use only information available when the order is placed.

The forecast set-up: direct H-step

We make a direct H-step-ahead forecast. To decide how much of a product to stock for day $d$, we stand at the origin $o = d - H$ (here $H = 7$ days, a weekly review) and may use only data up to $o$. Every feature derived from past demand is therefore built from one leak-safe primitive, computed within each product:

$$\texttt{avail} = \texttt{demand.shift}(H).$$

avail on day $d$ equals the demand on day $d-H$ — the most recent value known at the origin. Lags and rolling windows are measured backward from avail, so no feature can peek past the origin. Calendar and price columns describe the (known-in-advance) target day and need no shift. This shift(H) trick is the whole game: get it right and the backtest in Act 5 is trustworthy; get it wrong and every result is optimistic fiction.

The features we build

group	features	why it is known at decision time
calendar	dow, is_saturday (the closure), month, day, weekofyear, dayofyear, sin/cos_year, days_to_xmas	the target date is known in advance
price	price, log_price	a stable product attribute in this data
history (from avail)	lag_7/14/28, roll_mean/std/zero_{7,28}	shifted by $H$ → cannot see past the origin
intermittency	days_since_last_sale	Croston-style recency, from avail
identity / regime	StockCode, sb_class (native categoricals)	product attributes

days_to_xmas (circular distance to 25 Dec) hands the model the seasonal surge directly; sb_class tells it each product's demand regime. StockCode as a native LightGBM category lets a single global model learn per-product levels — no one-hot, no target encoding, no leakage.

HORIZON = 7                 # forecast/lead horizon (days) -> weekly review
LAGS = [7, 14, 28]
WINDOWS = [7, 28]
SERIES = "StockCode"

panel = pd.read_parquet(DATA/"panel.parquet"); panel["date"] = pd.to_datetime(panel["date"])
ps = pd.read_parquet(DATA/"sb_stats.parquet")

def _days_since_last(a):
    out = np.empty(len(a)); since = np.nan
    for i, v in enumerate(a):
        if np.isfinite(v) and v > 0: since = 0.0
        elif np.isfinite(since):     since += 1.0
        out[i] = since
    return out

def build_features(panel, ps, H, lags, windows):
    df = panel.sort_values([SERIES, "date"]).reset_index(drop=True).copy()
    dt = df["date"].dt
    # --- calendar (target day, known in advance) ---
    df["dow"] = dt.dayofweek
    df["is_saturday"] = (df["dow"] == 5).astype("int8")     # structural closure
    df["month"] = dt.month; df["day"] = dt.day
    df["weekofyear"] = dt.isocalendar().week.astype("int32")
    df["dayofyear"] = dt.dayofyear
    ang = 2*np.pi*df["dayofyear"]/365.25
    df["sin_year"] = np.sin(ang); df["cos_year"] = np.cos(ang)
    dxm = (df["dayofyear"] - 359).abs()                      # Dec 25 ~ doy 359
    df["days_to_xmas"] = np.minimum(dxm, 365 - dxm)
    feats = ["dow","is_saturday","month","day","weekofyear","dayofyear","sin_year","cos_year","days_to_xmas"]
    # --- price (stable product attribute) ---
    df["log_price"] = np.log(df["price"].clip(lower=1e-3))
    feats += ["price", "log_price"]
    # --- leak-safe availability: demand as of the origin ---
    avail = df.groupby(SERIES, observed=True)["demand"].shift(H)
    gav = avail.groupby(df[SERIES], observed=True)
    df["lag_0"] = avail; feats.append("lag_0")
    for k in lags:
        df[f"lag_{k}"] = gav.shift(k); feats.append(f"lag_{k}")
    for w in windows:
        df[f"roll_mean_{w}"] = gav.transform(lambda s: s.rolling(w, min_periods=1).mean())
        df[f"roll_std_{w}"]  = gav.transform(lambda s: s.rolling(w, min_periods=2).std())
        df[f"roll_zero_{w}"] = gav.transform(lambda s: s.eq(0).rolling(w, min_periods=1).mean())
        feats += [f"roll_mean_{w}", f"roll_std_{w}", f"roll_zero_{w}"]
    df["days_since_last_sale"] = gav.transform(lambda s: pd.Series(_days_since_last(s.to_numpy()), index=s.index))
    feats.append("days_since_last_sale")
    # --- identity + demand regime (native categoricals) ---
    df["sb_class"] = df[SERIES].map(ps["class"]).astype("category")
    df[SERIES] = df[SERIES].astype("category")
    cat_cols = [SERIES, "sb_class"]
    feats += cat_cols
    return df, feats, cat_cols

feat_df, FEATURES, CAT_COLS = build_features(panel, ps, HORIZON, LAGS, WINDOWS)
print(f"feature matrix: {len(feat_df):,} rows x {len(FEATURES)} features "
      f"({len(CAT_COLS)} categorical), H={HORIZON}")
feat_df[["date", SERIES, "description", "demand", "lag_0", "roll_mean_28",
         "days_since_last_sale", "sb_class"]].head()

feature matrix: 147,800 rows x 24 features (2 categorical), H=7

	date	StockCode	description	demand	lag_0	roll_mean_28	days_since_last_sale	sb_class
0	2009-12-01	15034	PAPER POCKET TRAVELING FAN	3	NaN	NaN	NaN	lumpy
1	2009-12-02	15034	PAPER POCKET TRAVELING FAN	0	NaN	NaN	NaN	lumpy
2	2009-12-03	15034	PAPER POCKET TRAVELING FAN	0	NaN	NaN	NaN	lumpy
3	2009-12-04	15034	PAPER POCKET TRAVELING FAN	0	NaN	NaN	NaN	lumpy
4	2009-12-05	15034	PAPER POCKET TRAVELING FAN	0	NaN	NaN	NaN	lumpy

The no-leakage check

lag_0 on day $d$ is defined as demand on day $d-H$. So plotting a product's actual demand and its lag_0 should show the lag_0 curve as the demand curve shifted right by exactly $H$ days — the features always trail reality, never anticipate it. If they lined up (or lag_0 led), we would have a leak. Here is the erratic signature item.

sc = "85123A"
one = feat_df[feat_df[SERIES] == sc].set_index("date").loc["2011-04-01":"2011-07-01"]
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(one.index, one["demand"], color=GREY, lw=.9, label="demand (actual, target day)")
ax.plot(one.index, one["lag_0"], color=RED, lw=1.4, label=f"lag_0 = demand at origin (d−{HORIZON})")
ax.plot(one.index, one["roll_mean_28"], color=BLUE, lw=2, label="roll_mean_28 (trailing, from origin)")
i = 30
ax.annotate("", xy=(one.index[i], one["demand"].iloc[i]),
            xytext=(one.index[i-HORIZON], one["demand"].iloc[i]),
            arrowprops=dict(arrowstyle="<->", color="black"))
ax.text(one.index[i-HORIZON//2], one["demand"].iloc[i]*1.05, f"H={HORIZON}d", ha="center", fontsize=9)
ax.set_title(f"No leakage: features trail actual demand by H days  ({sc})")
ax.set_ylabel("units/day"); ax.legend(fontsize=8); fig.tight_layout(); plt.show()

The red lag_0 is the grey demand moved right by $H$ days — visual proof that on any target day the model sees only origin-and-earlier demand.

Intermittency signals — handing the model the pattern

For lumpy products the gaps are informative. days_since_last_sale (recency) and roll_zero_28 (share of recent zero-days) encode them directly, both measured from the origin.

_lumpy = ps[ps["class"] == "lumpy"].sort_values("mean")
sc2 = _lumpy.index[len(_lumpy)//2]          # a mid-volume lumpy product
seg = feat_df[feat_df[SERIES] == sc2].set_index("date").loc["2011-01-01":"2011-12-01"]
fig, ax = plt.subplots(2, 1, figsize=(12, 6), sharex=True)
ax[0].bar(seg.index, seg["demand"], color=BLUE, width=1.0)
ax[0].set_ylabel("demand"); ax[0].set_title(f"Lumpy product {sc2} — {desc.get(sc2, '')[:34]}")
ax[1].plot(seg.index, seg["days_since_last_sale"], color=ORANGE, lw=1.2, label="days_since_last_sale")
axb = ax[1].twinx()
axb.plot(seg.index, seg["roll_zero_28"], color=PURP, lw=1.4, label="roll_zero_28")
ax[1].set_ylabel("days since sale", color=ORANGE); axb.set_ylabel("zero-share (28d)", color=PURP)
ax[1].legend(loc="upper left", fontsize=8); axb.legend(loc="upper right", fontsize=8)
fig.tight_layout(); plt.show()

When demand goes quiet, days_since_last_sale ramps up and roll_zero_28 climbs; a sale resets both. These are exactly the quantities Croston's method tracks — now ordinary columns a tree model can split on.

Missingness is structural — in two different ways

The NaNs above are expected rather than dirty, and they arise for two distinct reasons that are worth keeping apart.

No history yet. Every feature is built from demand available $H$ days before the target, so the opening rows of each product have nothing to look back at — $H$ rows for the shift itself, and more for a feature that then needs a further lag or a wider window. Each of these is a clean leading block: the NaN count is exactly the number of products times the length of that block.

No sale yet. days_since_last_sale stays undefined until a product has sold at least once in the available window, which for a slow mover can be hundreds of rows in. This NaN carries information — the product has never sold — which is not the same as a value being absent.

LightGBM handles NaN natively, learning at each split which way to send missing values, so both kinds are usable signal rather than rows to drop. The backtest trains only where history exists.

na = feat_df[FEATURES].isna().mean()
na = na[na > 0].sort_values()

# Two kinds of missing, and they are not the same thing. Warm-up NaNs form a leading block in every
# product's series -- the horizon shift, plus whatever extra a longer lag or wider window needs -- so
# the count is exactly n_products x block length. Anything failing that test is scattered through the
# series instead: a product that has not sold yet, which is informative rather than absent.
_n_prod  = feat_df[SERIES].nunique()
_row     = feat_df.sort_values([SERIES, "date"]).groupby(SERIES, observed=True).cumcount()
_deepest = {c: int(_row[feat_df[c].isna()].max()) for c in na.index}
_warmup  = {c: feat_df[c].isna().sum() == _n_prod * (_deepest[c] + 1) for c in na.index}

fig, ax = plt.subplots(figsize=(9, max(3, .35*len(na))))
ax.barh(na.index, na.values*100,
        color=[PURP if _warmup[c] else "#c2703d" for c in na.index])
ax.set_xlabel("% of rows missing")
ax.set_title("Structural missingness: warm-up (purple) vs. no sale yet (orange)")
fig.tight_layout(); plt.show()

for c in na.index[::-1]:
    kind = (f"warm-up: {_n_prod} x {_deepest[c]+1} leading rows" if _warmup[c]
            else "no sale yet: scattered through the series")
    print(f"  {c:22} {feat_df[c].isna().sum():>6,} NaN ({na[c]*100:4.1f}%)  "
          f"deepest at row {_deepest[c]:>3}  [{kind}]")
feat_df.to_parquet(DATA/"features.parquet", index=False)
print(f"saved features.parquet | {len(FEATURES)} features, {len(CAT_COLS)} categorical -> {DATA}")

  days_since_last_sale   15,892 NaN (10.8%)  deepest at row 738  [no sale yet: scattered through the series]
  lag_28                  7,000 NaN ( 4.7%)  deepest at row  34  [warm-up: 200 x 35 leading rows]
  lag_14                  4,200 NaN ( 2.8%)  deepest at row  20  [warm-up: 200 x 21 leading rows]
  lag_7                   2,800 NaN ( 1.9%)  deepest at row  13  [warm-up: 200 x 14 leading rows]
  roll_std_28             1,600 NaN ( 1.1%)  deepest at row   7  [warm-up: 200 x 8 leading rows]
  roll_std_7              1,600 NaN ( 1.1%)  deepest at row   7  [warm-up: 200 x 8 leading rows]
  roll_mean_28            1,400 NaN ( 0.9%)  deepest at row   6  [warm-up: 200 x 7 leading rows]
  roll_mean_7             1,400 NaN ( 0.9%)  deepest at row   6  [warm-up: 200 x 7 leading rows]
  lag_0                   1,400 NaN ( 0.9%)  deepest at row   6  [warm-up: 200 x 7 leading rows]
saved features.parquet | 24 features, 2 categorical -> C:\Users\user\dev\personal-website\Z-ML Operations Research-Predict-then-Optimize\data

Section 4 takeaways

One primitive — avail = demand.shift(H) — makes every history feature leak-safe, proven visually. The matrix mixes known-in-advance calendar/price signal, leak-safe history, explicit intermittency features, and native categoricals (product + demand regime).

Next — Section 5: the quantile forecast. We train LightGBM to predict the distribution of demand (a set of quantiles) by minimizing the pinball loss, and check its calibration.

5 · Forecasting the demand distribution

We do not predict a single number. For each product-day we predict a set of quantiles of demand — the conditional distribution — because the inventory decision in Section 6 needs a specific tail quantile, not the mean.

Quantile regression via the pinball loss

Ordinary regression minimizes squared error and recovers the conditional mean. To target the $\tau$-quantile instead, we minimize the pinball (quantile) loss

$$L_\tau(y,\hat y)=\begin{cases}\tau\,(y-\hat y) & y\ge \hat y\\[2pt](1-\tau)(\hat y-y) & y<\hat y\end{cases}$$

which penalizes under- and over-prediction asymmetrically. For $\tau=0.9$, an under-forecast costs 9× an over-forecast, so the minimizer sits where only 10% of demand exceeds it — exactly the 0.9-quantile. $\tau=0.5$ is symmetric → the median. No distribution is assumed; the data set the shape.

Gradient boosting (LightGBM) as the learner

The model is an additive ensemble of shallow trees, each fit to the gradient of the loss, so the ensemble descends the loss surface. Because the objective is pluggable, objective="quantile", alpha=τ makes it minimize the pinball loss directly. LightGBM also (a) splits our category-dtype StockCode/sb_class natively — no one-hot, no target-encoding leakage; (b) routes the structural warm-up NaNs from Section 4. We fit one model per quantile — a single global model across all products (borrowing strength across the catalogue), not one model per product.

A valid distribution out — non-crossing

Independently-fit quantiles can occasionally cross; we sort each row's predictions ascending and clip at 0, so every row is a non-crossing, non-negative demand distribution ready for the optimizer.

Honest evaluation — pinball and calibration

We use a temporal split (never shuffle time) with early stopping on validation pinball, and judge the model by average pinball loss and calibration: for a well-calibrated $\tau$-quantile, empirically $P(\text{demand}\le \hat q_\tau)\approx\tau$. That coverage table is what makes the downstream newsvendor decision trustworthy. (The full rolling-origin backtest is Section 7.)

Why not a point or parametric model? A point forecast rounds lumpy demand toward 0 → chronic stockouts; a parametric model (Poisson/NegBin) forces a shape the erratic data may not follow. Quantile boosting is nonparametric and targets exactly the quantile the decision needs.

import os
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")   # avoid duplicate-OpenMP kernel death (Windows)
import lightgbm as lgb

QUANTILES = [0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975, 0.99]
LGB_PARAMS = dict(objective="quantile", n_estimators=400, learning_rate=0.05, num_leaves=63,
                  min_child_samples=50, subsample=0.8, subsample_freq=1, colsample_bytree=0.8,
                  verbosity=-1, n_jobs=-1)

def pinball(y, q, tau):
    d = np.asarray(y, float) - np.asarray(q, float)
    return float(np.mean(np.maximum(tau*d, (tau-1)*d)))

class QuantileGBM:
    """One LightGBM per quantile; predictions are made non-crossing and non-negative."""
    def __init__(self, quantiles, params):
        self.quantiles = sorted(quantiles); self.params = params; self.models = {}
    def fit(self, X, y, Xv, yv):
        for tau in self.quantiles:
            m = lgb.LGBMRegressor(alpha=tau, random_state=7, **self.params)
            m.fit(X, y, eval_set=[(Xv, yv)], eval_metric="quantile",
                  callbacks=[lgb.early_stopping(50, verbose=False), lgb.log_evaluation(0)])
            self.models[tau] = m
        return self
    def predict(self, X):
        P = np.column_stack([self.models[t].predict(X) for t in self.quantiles])
        P = np.clip(P, 0, None); P.sort(axis=1)
        return pd.DataFrame(P, columns=self.quantiles, index=X.index)

# temporal split: train on the past, validate on the last 28 days
cutoff = feat_df["date"].max() - pd.Timedelta(days=28)
tr, va = feat_df[feat_df["date"] <= cutoff], feat_df[feat_df["date"] > cutoff]
model = QuantileGBM(QUANTILES, LGB_PARAMS).fit(tr[FEATURES], tr["demand"], va[FEATURES], va["demand"])

def evaluate(model, frame):
    y = frame["demand"].to_numpy(); P = model.predict(frame[FEATURES])
    return pd.DataFrame([dict(quantile=t, pinball=pinball(y, P[t].to_numpy(), t),
                              coverage=float(np.mean(y <= P[t].to_numpy()))) for t in model.quantiles])

val = evaluate(model, va)
print(f"trained {len(QUANTILES)} quantile models | validation after {cutoff.date()} "
      f"| mean pinball {val['pinball'].mean():.3f}")
val.round(3)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\lightgbm\sklearn.py:1106: LGBMDeprecationWarning: The argument 'eval_set' is deprecated, use 'eval_X' and 'eval_y' instead.
  eval_set = _validate_eval_set_Xy(eval_set=eval_set, eval_X=eval_X, eval_y=eval_y)

trained 9 quantile models | validation after 2011-11-11 | mean pinball 17.787

	quantile	pinball	coverage
0	0.050	2.703	0.296
1	0.100	5.333	0.306
2	0.250	12.392	0.462
3	0.500	22.049	0.676
4	0.750	27.700	0.837
5	0.900	26.112	0.932
6	0.950	23.609	0.963
7	0.975	21.367	0.977
8	0.990	18.814	0.990

Reading the validation table

Three columns, and the middle one is the most important diagnostic in the whole forecast:

• quantile ($\tau$) — the level each of the 9 models targets (0.5 = median, 0.99 = 99th percentile).
• pinball — the mean pinball loss on the held-out days, in demand units (lower = better); the quantity each model minimizes.
• coverage — the empirical hit rate $P(\text{demand}\le\hat q_\tau)$. The calibration test: for a well-calibrated $\tau$-quantile, coverage should $\approx\tau$.

Reading coverage against quantile reveals three regimes:

$\tau$	coverage	what it means
0.05, 0.10	≈ 0.30	pinned at the zero-share — the atom at zero
0.25 – 0.75	above $\tau$	mildly conservative (predicts a touch high)
0.90 – 0.99	≈ $\tau$	calibrated, essentially exact at 0.975–0.99

1. Low quantiles pinned at the zero-share (~0.30). About 30% of validation days have zero demand, so the true 5th/10th percentiles genuinely are 0. The model predicts $\hat q=0$ and its coverage is $P(\text{demand}\le 0)=$ the zero-share. Any quantile below the zero-share collapses onto 0 — correct behaviour, not a defect (a probability atom at zero).
2. The middle runs slightly high. The median model covers ~0.68 of days rather than 0.50: mild over-forecasting. For a newsvendor that just means marginally more safety stock than nominal — a safe direction to err.
3. The upper tail is calibrated — and that is what matters. The inventory decision reads a high quantile, and there coverage tracks $\tau$ almost exactly (0.975 → ~0.98, 0.99 → ~0.99). The model is trustworthy precisely where we will spend it.

Why pinball is largest in the middle (~28 at $\tau=0.75$ vs ~3 at $\tau=0.05$): the low quantiles predict ≈0 and accrue only small errors across the many zero days, while the mid-quantiles sit where demand is most variable and hardest to pin down. The absolute scale is large only because demand here is in the hundreds of units — pinball is measured in demand units. The reliability diagram below is the same story, drawn.

Calibration — the headline

For a well-calibrated $\tau$-quantile the empirical coverage $P(\text{demand}\le\hat q_\tau)$ should sit on the 45° line. Watch what the ~48% zero-share does to the low quantiles.

zero_share = float((va["demand"] == 0).mean())
fig, ax = plt.subplots(1, 2, figsize=(13, 5))
ax[0].plot([0, 1], [0, 1], "--", color=GREY, label="perfect calibration")
ax[0].plot(val["quantile"], val["coverage"], "o-", color=BLUE, label="empirical coverage")
ax[0].axhline(zero_share, color=RED, ls=":", lw=1.5, label=f"zero-share = {zero_share:.2f}")
ax[0].set_xlabel("target quantile τ"); ax[0].set_ylabel("coverage  P(demand ≤ q̂)")
ax[0].set_title("Reliability diagram"); ax[0].legend(fontsize=8)
ax[1].bar(val["quantile"].astype(str), val["pinball"], color=PURP, width=.7)
ax[1].set_xlabel("quantile τ"); ax[1].set_ylabel("mean pinball loss")
ax[1].set_title("Pinball loss by quantile"); ax[1].tick_params(axis="x", rotation=45)
fig.tight_layout(); plt.show()
print("upper-tail coverage:", val[val["quantile"] >= 0.75].set_index("quantile")["coverage"].round(3).to_dict())

upper-tail coverage: {0.75: 0.837, 0.9: 0.932, 0.95: 0.963, 0.975: 0.977, 0.99: 0.99}

Reading it. Every quantile below the zero-share collapses onto $\hat q=0$, so its coverage is $P(\text{demand}\le 0)=$ the zero-share — the flat left arm (≈0.30 here, lower than the ~48% overall because this validation window is peak season). That is not miscalibration: when a product is zero on a large share of days, its true low percentiles genuinely are 0 (a probability atom at zero). Above that, coverage climbs and at the extreme upper tail (0.975–0.99), where the inventory decision lives, it is essentially exact. In the mid-range the model runs slightly conservative (coverage a touch above target) — for a newsvendor that just means marginally higher service than nominal, a safe direction to err. The model is trustworthy exactly where we will spend it.

A fan chart — the distribution, per day

For the signature product over the last months: actual demand with the predicted quantile bands. The bands widen and narrow with the model's day-to-day uncertainty.

sc = "85123A"
P_all = model.predict(feat_df[FEATURES])
J = feat_df[["date", SERIES, "description", "demand"]].join(P_all)
one = J[J[SERIES] == sc].set_index("date").sort_index().iloc[-120:]
fig, ax = plt.subplots(figsize=(12, 4.5))
ax.fill_between(one.index, one[0.05], one[0.95], color=BLUE, alpha=.15, label="q0.05–q0.95")
ax.fill_between(one.index, one[0.25], one[0.75], color=BLUE, alpha=.30, label="q0.25–q0.75")
ax.plot(one.index, one[0.5], color=BLUE, lw=1.5, label="median (q0.5)")
ax.plot(one.index, one[0.9], color=GREEN, lw=1.4, ls="--", label="q0.9 (an upper quantile)")
ax.plot(one.index, one["demand"], color=RED, lw=.9, marker=".", ms=3, label="actual")
ax.axvline(cutoff, color="#718096", ls=":", lw=1); ax.text(cutoff, ax.get_ylim()[1]*.9, " validation →", fontsize=8)
ax.set_title(f"Predicted demand distribution vs actual — {sc} {one['description'].iloc[0][:30]}")
ax.set_ylabel("units/day"); ax.legend(fontsize=8, ncol=2); fig.tight_layout(); plt.show()

A note on that spike. The one actual bar that pierces the top band (~1,300 units) is not an error — it is a genuine bulk wholesale order. This item's all-time peak, 3,114 units on 2011-01-11, came from just two invoices (1,930 and 1,010 units) placed by two different customers at normal prices, against a routine level of ~80/day. Because this is a wholesaler, a single retailer stocking up can order a thousand units at once. (The one true data anomaly — a 4,000-unit line at £0.00 with no customer — was already removed by the Price > 0 filter in Section 2.) Such one-off bulk orders are the essence of lumpy/erratic demand: no daily model can anticipate a single customer's decision, so they remain tail risk even at q0.95. This is exactly why we forecast a distribution and carry safety stock for the routine variation, while accepting that extreme orders sit beyond what any SKU-day forecast can catch.

What drives the forecast

Gain-based importance for the median model. Recent-history and the intermittency/seasonal features should dominate.

imp = pd.Series(model.models[0.5].booster_.feature_importance(importance_type="gain"),
                index=FEATURES).sort_values().tail(15)
fig, ax = plt.subplots(figsize=(9, 6))
ax.barh(imp.index, imp.values, color=GREEN)
ax.set_xlabel("gain"); ax.set_title("Top features (median model, gain importance)")
fig.tight_layout(); plt.show()

Interpretability — SHAP values (exact TreeSHAP)

Gain importance ranks features globally, but it is biased toward high-cardinality splits (like the 200-level StockCode) and says nothing about direction. SHAP (SHapley Additive exPlanations) fixes both: it distributes each prediction among its features so the contributions sum to the prediction, giving a signed, per-row attribution grounded in cooperative game theory. LightGBM computes exact TreeSHAP itself via predict(pred_contrib=True), so we get it with no extra dependency. Below is the median model on a sample of rows.

The beeswarm plots one dot per row for each feature: horizontal position = its SHAP value (impact on the predicted median demand, in units), colour = the feature's value (low → high).

booster = model.models[0.5].booster_
Xs = tr[FEATURES].sample(min(4000, len(tr)), random_state=7)
contrib = booster.predict(Xs, pred_contrib=True)     # exact TreeSHAP: (n, n_features + 1)
sv = contrib[:, :-1]                                  # per-feature SHAP; last col = base value
order = np.argsort(np.abs(sv).mean(0))[::-1][:12]     # top-12 by mean |SHAP|
rng = np.random.RandomState(7)

fig, ax = plt.subplots(figsize=(10, 7))
for rank, fi in enumerate(order):
    fv = Xs.iloc[:, fi]
    if str(fv.dtype) == "category":
        col = fv.cat.codes.to_numpy().astype(float)
    else:
        col = pd.to_numeric(fv, errors="coerce").to_numpy()
        col = np.where(np.isnan(col), np.nanmedian(col), col)
    cnorm = np.argsort(np.argsort(col)) / max(len(col)-1, 1)     # rank-normalized colour
    y = (len(order)-1-rank) + (rng.rand(len(fv)) - 0.5) * 0.6
    scat = ax.scatter(sv[:, fi], y, c=cnorm, cmap="coolwarm", s=6, alpha=.5, linewidths=0)
ax.set_yticks(range(len(order))); ax.set_yticklabels([FEATURES[i] for i in order][::-1])
ax.axvline(0, color=GREY, lw=.8)
ax.set_xlabel("SHAP value  (impact on predicted median demand, units)")
ax.set_title("SHAP summary (beeswarm) — median model")
cb = plt.colorbar(scat, ax=ax); cb.set_label("feature value (low → high)")
fig.tight_layout(); plt.show()

Read a row like this: dots to the right push the forecast up, left push it down; red means a high feature value, blue low. The recent-history features (lag_0, roll_mean_28/7, lag_14/7) show red on the right — high recent demand raises the forecast, as it should. dow captures the weekly rhythm (high day-of-week values — the weekend and the Saturday closure — push the forecast down); roll_zero_28 pushes down when recent zero-days pile up; dayofyear/sin_year carry the seasonal ramp; and StockCode encodes the per-product level. Nothing here is spurious — which is the audit SHAP buys us. The dependence plots isolate the shape of two features.

fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
for a, feat in zip(ax, ["days_to_xmas", "roll_mean_28"]):
    fi = FEATURES.index(feat)
    a.scatter(pd.to_numeric(Xs[feat], errors="coerce"), sv[:, fi], s=8, alpha=.35, color=PURP)
    a.axhline(0, color=GREY, lw=.8)
    a.set_xlabel(feat); a.set_ylabel(f"SHAP value ({feat}, units)")
    a.set_title(f"Dependence: {feat}")
fig.tight_layout(); plt.show()

roll_mean_28 is essentially monotone — a higher recent average lifts the forecast (with diminishing returns at the very top). days_to_xmas is subtler, and more honest than a naive "closer = higher": its contribution turns positive in the weeks before Christmas (the wholesale surge, roughly 2–4 weeks out) and negative right at and after 25 Dec, when this B2B gift wholesaler's orders collapse. The model recovered the true pre-Christmas timing of the peak rather than a spike on the day itself — exactly the structure SHAP lets us verify before a forecast drives a decision.

Section 5 takeaways

Quantile gradient boosting gives a full demand distribution per product-day with no distributional assumption — the right object for lumpy, zero-inflated demand. Calibration is exact in the upper tail (the service-level region the decision uses); the flat low-quantile arm is the correct signature of the zero-atom, not a defect.

Next — Section 6: the decision. We turn each day's distribution into an order with the newsvendor rule, using per-product economics derived from the real GBP prices — a richer setup than a single fixed service level.

6 · From distribution to decision — the newsvendor

The forecaster gives a demand distribution per product-day. The inventory decision reads a single quantile off it. Which quantile, and why, is the classic newsvendor model.

The rule

Order $Q$ before demand $D\sim F$ is realized. Two asymmetric costs:

• overstock $C_o$ per unit ordered but unsold (holding);
• understock $C_u$ per unit of demand unmet (lost margin).

Expected cost $C(Q)=C_o\,\mathbb{E}[(Q-D)^+]+C_u\,\mathbb{E}[(D-Q)^+]$. One more unit helps iff demand exceeds $Q$ (probability $1-F(Q)$, saving $C_u$) and is wasted otherwise (probability $F(Q)$, costing $C_o$). Order more while benefit ≥ cost:

$$C_u\big(1-F(Q)\big)\ge C_o\,F(Q)\;\Longrightarrow\; \boxed{\,Q^*=F^{-1}\!\left(\frac{C_u}{C_u+C_o}\right).}$$

The optimal order is the demand quantile at the critical fractile $C_u/(C_u+C_o)$ — which is the target service level. (Fractile is a synonym for quantile/percentile: the $\tau$-fractile is the demand level that actual demand falls at or below with probability $\tau$ — so a critical fractile of 0.86 means an 86% service level, i.e. a 14% chance of a stockout.) We read it straight off the predicted quantiles (interpolating the grid).

Predict and optimize are the same objective. With $\tau=C_u/(C_u+C_o)$, the pinball loss the model minimized equals the newsvendor cost up to the constant $(C_u+C_o)$: $(C_u+C_o)\,L_\tau(D,Q)=C_u(D-Q)^+ + C_o(Q-D)^+$. Training at the critical fractile is minimizing expected inventory cost — smart predict-then-optimize.

Per-product economics from the real prices

The dataset has selling price but not cost/margin/holding, so we infer $C_u,C_o$ with two stated assumptions (any real deployment supplies these from accounting):

• $C_u = \text{margin\_rate}\times\text{price}$ — a stockout forgoes the gross margin, not the whole price. We take margin_rate = 0.5 (gift items carry 40–60% margins; 50% is a defensible midpoint). We omit goodwill/lost-future-sales, so this is conservative.
• $C_o = £0.10$ per unit (flat) — a leftover non-perishable just carries to next period, so the period cost is holding. Physical storage scales with an item's size, not price, and for these cheap items (median £1.25) the cost-of-capital piece is negligible — so a flat per-unit charge is the honest model.

The consequence is the whole point: since $C_u$ scales with price but $C_o$ is flat,

$$\text{fractile}=\frac{0.5\,p}{0.5\,p+0.10}$$

rises with price — expensive items are stocked to a higher service level (a stockout wastes more margin, while holding costs the same), cheap items are allowed to stock out more. A single fixed fractile (as in a textbook example) cannot express that; real prices can.

import warnings, pulp
warnings.filterwarnings("ignore")

MARGIN_RATE, HOLDING = 0.5, 0.10
price_by_prod = panel.groupby("StockCode")["price"].first()
cu_prod = MARGIN_RATE * price_by_prod          # understock cost per unit (lost margin)
co_prod = HOLDING                              # overstock cost per unit (flat holding)
fractile = cu_prod / (cu_prod + co_prod)       # per-product critical fractile

summary = pd.DataFrame({"price": price_by_prod.round(2), "Cu": cu_prod.round(2),
                        "fractile": fractile.round(3)})
print(f"critical fractile spans {fractile.min():.2f} – {fractile.max():.2f} "
      f"(median {fractile.median():.2f}) across the 200 products")
summary.sort_values("price").iloc[[0, 100, -1]]

critical fractile spans 0.38 – 0.98 (median 0.86) across the 200 products

	price	Cu	fractile
StockCode
20668	0.12	0.06	0.375
84947	1.25	0.62	0.862
22423	12.75	6.38	0.985

pp = np.linspace(0.1, price_by_prod.quantile(.99), 300)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
for m, c in zip([0.4, 0.5, 0.6], [GREY, BLUE, GREEN]):
    ax[0].plot(pp, m*pp/(m*pp + HOLDING), color=c, label=f"margin {m:.0%}")
ax[0].scatter(price_by_prod.clip(upper=pp.max()), fractile, s=10, color=RED, alpha=.4, label="products (margin 50%)")
ax[0].set_xlabel("unit price (£)"); ax[0].set_ylabel("critical fractile = service level")
ax[0].set_title("Service level rises with price (+ margin sensitivity)"); ax[0].legend(fontsize=8)
ax[1].hist(fractile, bins=30, color=PURP)
ax[1].set_xlabel("critical fractile"); ax[1].set_ylabel("# products")
ax[1].set_title("Distribution of per-product service levels")
fig.tight_layout(); plt.show()

The order and its safety stock

For a decision snapshot — every product on the last available date — we read each product's order at its own fractile, and compare to the median. The gap is safety stock; note it grows with price.

def scenario_weights(quantiles):
    q = np.array(sorted(quantiles)); edges = np.concatenate([[0], (q[:-1]+q[1:])/2, [1]]); return np.diff(edges)

def order_at_fractile(P, fractiles, quantiles):
    q = np.array(sorted(quantiles)); V = P[sorted(quantiles)].to_numpy()
    return np.array([np.interp(fractiles[i], q, V[i]) for i in range(len(fractiles))])

snap = feat_df[feat_df["date"] == feat_df["date"].max()].copy()
Psnap = model.predict(snap[FEATURES]); Psnap.index = snap["StockCode"].to_numpy()
frac_snap = fractile.reindex(Psnap.index).to_numpy()
order = pd.Series(order_at_fractile(Psnap, frac_snap, QUANTILES), index=Psnap.index)
median = Psnap[0.5]
price_snap = price_by_prod.reindex(Psnap.index)

fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
sc = ax[0].scatter(median, order, c=np.log10(price_snap), cmap="viridis", s=16, alpha=.7)
lim = max(order.max(), median.max())*1.05
ax[0].plot([0, lim], [0, lim], "--", color=GREY, label="order = median")
ax[0].set_xlabel("median demand q0.5"); ax[0].set_ylabel("newsvendor order")
ax[0].set_title("Order sits above median — more so for pricier items"); ax[0].legend(fontsize=8)
plt.colorbar(sc, ax=ax[0], label="log10 price")
ax[1].scatter(price_snap, (order-median).clip(lower=0), s=14, color=GREEN, alpha=.6)
ax[1].set_xlabel("unit price (£)"); ax[1].set_ylabel("safety stock (order − median)")
ax[1].set_title("Safety stock grows with price"); fig.tight_layout(); plt.show()
print(f"snapshot {snap['date'].iloc[0].date()} | total order {order.sum():.0f} vs total median {median.sum():.0f} "
      f"(+{(order.sum()/median.sum()-1):.0%} safety stock)")

snapshot 2011-12-09 | total order 17148 vs total median 5318 (+222% safety stock)

A worked example — the buyer's order sheet

How a business actually uses this. Every review cycle (here, weekly) a purchasing planner runs the pipeline and receives an order sheet: for each product, the number of units to have in stock for the coming period. That number is the newsvendor order-up-to level $Q^*$ — the demand quantile at the product's service level. In a live system the purchase order placed is order-up-to − on-hand stock (you top up to the target); here we show the target level itself. When the warehouse or the week's budget can't fit every product's ideal order, the planner applies the capacity LP (next subsection) to allocate the limited space. The loop then repeats next cycle with fresh data — this is exactly how a replenishment system runs in practice.

Below is the sheet for the snapshot day, for a selection of recognizable products (cheap → pricey). In the chart the grey bar is the plausible demand range (q0.05–q0.95), the blue band the interquartile range, the blue dot the median (what a naive point forecast would order), and the ★ is the newsvendor order. The gap between dot and star is the safety stock the economics justify.

cands = ["85123A","84077","85099B","84879","47566","22086","23084","20725","22910","21212","23203","22386"]
picks = [sc for sc in cands if sc in Psnap.index]
rows = []
for sc in picks:
    r = Psnap.loc[sc]; f = float(fractile[sc])
    od = float(np.interp(f, sorted(QUANTILES), r[sorted(QUANTILES)].to_numpy()))
    rows.append(dict(code=sc, product=str(desc[sc]).strip()[:30], price=float(price_by_prod[sc]),
                     service=f, q05=float(r[0.05]), q25=float(r[0.25]), median=float(r[0.5]),
                     q75=float(r[0.75]), q95=float(r[0.95]), order=od))
T = pd.DataFrame(rows).sort_values("price").reset_index(drop=True)
T["safety"] = T["order"] - T["median"]
print(f"recommended inventory for {snap['date'].iloc[0].date()} (order for demand {HORIZON} days out):\n")
print(T[["code","product","price","service","median","order","safety"]].to_string(index=False,
      formatters={"price":"£{:.2f}".format, "service":"{:.0%}".format, "median":"{:.0f}".format,
                  "order":"{:.0f}".format, "safety":"{:.0f}".format}))

recommended inventory for 2011-12-09 (order for demand 7 days out):

  code                        product price service median order safety
 84077 WORLD WAR 2 GLIDERS ASSTD DESI £0.29     59%     96   117     21
 21212 PACK OF 72 RETROSPOT CAKE CASE £0.55     73%     77   124     48
 20725        LUNCH BAG RED RETROSPOT £1.65     89%     53   164    111
 84879  ASSORTED COLOUR BIRD ORNAMENT £1.69     89%    127   342    215
 22386        JUMBO BAG PINK POLKADOT £1.95     91%     36   169    133
85099B        JUMBO BAG RED RETROSPOT £1.95     91%     99   363    263
 23084             RABBIT NIGHT LIGHT £2.08     91%     49   318    269
 23203        JUMBO BAG VINTAGE DOILY £2.08     91%     74   250    176
85123A WHITE HANGING HEART T-LIGHT HO £2.95     94%    129   372    243
 22086 PAPER CHAIN KIT 50'S CHRISTMAS £2.95     94%    154   350    196
 22910 PAPER CHAIN KIT VINTAGE CHRIST £2.95     94%    109   245    136
 47566                  PARTY BUNTING £4.95     96%     22   158    136

Tc = T.head(8); y = np.arange(len(Tc))
lo = np.maximum(Tc.q05.values, 0.5); q25 = np.maximum(Tc.q25.values, 0.5)
med = np.maximum(Tc["median"].values, 0.5); od = Tc["order"].values
fig, ax = plt.subplots(figsize=(11, 5.5))
ax.hlines(y, lo, Tc.q95.values, color=GREY, lw=2.5, alpha=.6, label="q0.05–q0.95")
ax.hlines(y, q25, Tc.q75.values, color=BLUE, lw=8, alpha=.4, label="interquartile")
ax.scatter(med, y, color=BLUE, s=45, zorder=3, label="median (point forecast)")
ax.scatter(od, y, color=RED, marker="*", s=210, zorder=4, label="newsvendor order")
for i in range(len(Tc)):
    ax.annotate(f"{od[i]:.0f} ({Tc.service.iloc[i]:.0%})", (od[i], i), xytext=(8, 7),
                textcoords="offset points", fontsize=8, color=RED)
ax.set_yticks(y); ax.set_yticklabels([f"{p}  £{pr:.2f}" for p, pr in zip(Tc["product"], Tc.price)], fontsize=8)
ax.set_xscale("log"); ax.set_xlabel("units of demand / order  (log scale)")
ax.set_title(f"Inventory decision per product — {snap['date'].iloc[0].date()}")
ax.legend(fontsize=8, loc="lower right"); fig.tight_layout(); plt.show()

A planner reads this sheet directly: the order column is the target stock level (net off on-hand to get the purchase quantity), service is the intended in-stock probability, and safety is the buffer above the expected (median) demand. Two things stand out. First, the order (★) always sits at or above the median, and for the lumpier, pricier items it lands far out in the right tail — their high service level and heavy-tailed demand both demand a large buffer. Second, products with similar median forecasts can warrant very different orders once the distribution and per-product economics are taken into account — exactly the information a single point forecast discards. This is the deliverable: not a forecast, but a decision a person can execute.

The capacity linear program

The per-product $Q^*$ above is unconstrained. A shared warehouse capacity couples products — you cannot grant each its ideal order. We represent each product's predicted quantiles as a scenario distribution and solve, with PuLP,

$$\min_{Q,u,o}\ \sum_p\sum_i w_i\big(C_{u,p}\,u_{p,i}+C_o\,o_{p,i}\big)\ \text{ s.t. }\ u_{p,i}\ge d_{p,i}-Q_p,\ o_{p,i}\ge Q_p-d_{p,i},\ \textstyle\sum_p Q_p\le \text{capacity},\ Q,u,o\ge 0.$$

Now the per-product $C_{u,p}$ makes the LP protect high-margin (pricey) items when space is scarce. Slack capacity → the closed form; binding → the constraint's shadow price values a unit of space.

What is a linear program?

A linear program (LP) optimizes a linear objective subject to linear constraints. It has three ingredients:

• decision variables — the quantities you choose (here each product's order $Q_p$, plus auxiliary $u_{p,i}, o_{p,i}$);
• a linear objective to minimize or maximize (here total expected cost);
• linear constraints defining what is feasible (the shortfall/overage definitions, $\sum_p Q_p\le$ capacity, and non-negativity).

Geometrically the constraints carve out a convex feasible region (a polytope), and a linear objective always attains its optimum at a vertex (corner) of it. That structure is why LPs solve exactly and globally — no local minima: the simplex method walks the vertices, interior-point methods cut through the middle. We use the open-source CBC solver via PuLP.

Why the newsvendor cost fits an LP. The per-product expected cost $\sum_i w_i\big[C_u(d_i-Q)^+ + C_o(Q-d_i)^+\big]$ is piecewise-linear and convex in $Q$, but the $(\cdot)^+$ (“max with 0”) is not itself linear. The standard trick linearizes it: introduce a shortfall $u_i\ge d_i-Q$ and an overage $o_i\ge Q-d_i$ (both $\ge 0$) and let the objective minimize them — at the optimum $u_i=(d_i-Q)^+$ and $o_i=(Q-d_i)^+$ automatically, because nothing pushes them above their lower bounds. The convex cost becomes a pure LP, coupled across products only by the single shared capacity row.

What is a shadow price?

Every constraint in an LP has a shadow price (its dual value): the rate at which the optimal objective would improve if you relaxed that constraint by one unit. For the capacity row, the shadow price is exactly the marginal value of one more unit of warehouse space — the £ of expected cost you would save by storing one extra unit, evaluated at the current optimum.

Three properties, all visible in the capacity sweep above:

• Complementary slackness. If the constraint is slack (you have more space than you'd use anyway), its shadow price is 0 — extra space is worthless. That is the flat right-hand end of the curve, past ~17,000 units.
• Binding → positive, and decreasing as capacity loosens. The shadow price is a step function (each drop is where the next product stops being rationed): £0.30 per unit when space is tight (6,000 units) down to £0 when it is ample. Scarcer space → each unit is worth more.
• Bounded by $C_u$. One extra unit can save at most a full understock cost, never more.

Economically the shadow price is a decision number: it is the most you should pay to rent or build one more unit of capacity. If a unit of shelf space costs less than the shadow price, expanding pays; if more, it does not. That "price of a constraint" falls straight out of LP duality — a manager gets it for free alongside the orders.

def capacity_lp(P, cu_vec, co, capacity, quantiles):
    skus = list(P.index); n, K = len(skus), len(quantiles)
    w = scenario_weights(quantiles); D = P[sorted(quantiles)].to_numpy()
    prob = pulp.LpProblem("cap", pulp.LpMinimize)
    Q = pulp.LpVariable.dicts("Q", range(n), lowBound=0)
    u = pulp.LpVariable.dicts("u", (range(n), range(K)), lowBound=0)
    o = pulp.LpVariable.dicts("o", (range(n), range(K)), lowBound=0)
    prob += pulp.lpSum(w[i]*(cu_vec[s]*u[s][i] + co*o[s][i]) for s in range(n) for i in range(K))
    for s in range(n):
        for i in range(K):
            prob += u[s][i] >= float(D[s, i]) - Q[s]
            prob += o[s][i] >= Q[s] - float(D[s, i])
    cap = pulp.lpSum(Q[s] for s in range(n)) <= capacity; prob += cap
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    orders = pd.Series([Q[s].value() or 0.0 for s in range(n)], index=skus)
    return orders, (abs(cap.pi) if cap.pi is not None else 0.0)

def exp_cost(orders, P, cu_vec, co, quantiles):
    w = scenario_weights(quantiles); D = P[sorted(quantiles)].to_numpy()
    diff = D - orders.reindex(P.index).to_numpy().reshape(-1, 1)
    return float(((cu_vec.reshape(-1, 1)*np.maximum(diff, 0) + co*np.maximum(-diff, 0)) * w).sum())

cu_vec = cu_prod.reindex(Psnap.index).to_numpy()
U = float(order.sum())
caps = np.linspace(0.35, 1.1, 10) * U
rows = [(c,) + (lambda oc, sh: (exp_cost(oc, Psnap, cu_vec, HOLDING, QUANTILES), oc.sum(), sh))(
            *capacity_lp(Psnap, cu_vec, HOLDING, c, QUANTILES)) for c in caps]
res = pd.DataFrame(rows, columns=["capacity", "exp_cost", "ordered", "shadow"])

fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
ax[0].plot(res.capacity, res.exp_cost, "o-", color=RED); ax[0].axvline(U, ls=":", color=GREEN, label=f"unconstrained {U:.0f}")
ax[0].set_xlabel("capacity (units)"); ax[0].set_ylabel("expected cost (£)"); ax[0].set_title("Expected cost vs capacity"); ax[0].legend(fontsize=8)
ax[1].plot(res.capacity, res.shadow, "s-", color=PURP); ax[1].axhline(0, color=GREY, lw=.8)
ax[1].set_xlabel("capacity (units)"); ax[1].set_ylabel("shadow price (£/unit)"); ax[1].set_title("Marginal value of capacity")
fig.tight_layout(); plt.show()

Rationing under scarcity — and sensitivity

At a binding capacity the LP chooses which products to cut. With per-product economics it protects the high-margin items: points below the 45° line are rationed, coloured by price.

cap = 0.6 * U
alloc, shadow = capacity_lp(Psnap, cu_vec, HOLDING, cap, QUANTILES)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
sc = ax[0].scatter(order, alloc.reindex(order.index), c=np.log10(price_snap), cmap="viridis", s=16, alpha=.7)
lim = order.max()*1.05; ax[0].plot([0, lim], [0, lim], "--", color=GREY, label="fully served")
ax[0].set_xlabel("unconstrained order"); ax[0].set_ylabel(f"allocated (capacity {cap:.0f})")
ax[0].set_title(f"Rationing at 60% capacity (shadow £{shadow:.2f})"); ax[0].legend(fontsize=8)
plt.colorbar(sc, ax=ax[0], label="log10 price")
# sensitivity: total unconstrained order under different assumptions
grid = []
for m in [0.4, 0.5, 0.6]:
    for h in [0.05, 0.10, 0.20]:
        fr = (m*price_snap)/(m*price_snap + h)
        tot = order_at_fractile(Psnap, fr.reindex(Psnap.index).to_numpy(), QUANTILES).sum()
        grid.append((m, h, tot))
G = pd.DataFrame(grid, columns=["margin", "holding", "order"]).pivot(index="margin", columns="holding", values="order")
im = ax[1].imshow(G.values, cmap="Blues", aspect="auto")
ax[1].set_xticks(range(3)); ax[1].set_xticklabels([f"£{h}" for h in [0.05,0.10,0.20]])
ax[1].set_yticks(range(3)); ax[1].set_yticklabels([f"{m:.0%}" for m in [0.4,0.5,0.6]])
ax[1].set_xlabel("holding cost"); ax[1].set_ylabel("margin rate"); ax[1].set_title("Total order vs assumptions (sensitivity)")
for i in range(3):
    for j in range(3): ax[1].text(j, i, f"{G.values[i,j]:.0f}", ha="center", va="center", fontsize=9)
fig.tight_layout(); plt.show()

The sensitivity grid shows the total order moving sensibly with the assumptions — higher margin or lower holding → higher service → more stock — with no cliff or pathology. Nothing hinges on a magic number; the shape (service rising with price, cost convex in capacity) is what carries the analysis.

Section 6 takeaways

• The newsvendor rule converts each predicted distribution into an order at that product's critical fractile, carrying more safety stock for pricier items — the payoff of real per-product economics.
• The capacity LP generalizes it to a shared constraint, pricing space via the shadow price and rationing to protect high-margin products.
• Predict-then-optimize, end to end: distribution → decision.

Next — Section 7: the backtest. We replay the whole loop over rolling origins and measure the realized cost of this policy against a point-forecast policy — the dollar value of modeling uncertainty.

7 · Backtest — the value of modeling uncertainty

Everything so far is a claim. The backtest tests it on held-out data by replaying the whole forecast→decision loop over rolling origins and scoring the realized cost in £.

Design. We walk the forecast origin forward over the last N_FOLDS windows of FOLD_H days. For each fold we train only on data before the window, forecast the window, turn the forecast into orders under three policies, and charge the actual newsvendor cost against realized demand — using the per-product economics from Section 6 ($C_{u,p}=0.5\times\text{price}$, $C_o=£0.10$):

• point (mean) — a good L2 point forecaster; order the expected demand;
• median (q0.5) — order the predicted median;
• newsvendor — order each product's quantile at its critical fractile.

Same features, same model family — only the decision differs. The £ gap between point and newsvendor is the value of modeling uncertainty. (Only the median and per-product fractile drive the orders, so we fit a compact quantile grid per fold.)

BT_QUANTILES = [0.25, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99]
N_FOLDS, FOLD_H = 3, 28

def fold_windows(max_date, n, h):
    out = [(max_date - pd.Timedelta(days=k*h) - pd.Timedelta(days=h-1),
            max_date - pd.Timedelta(days=k*h)) for k in range(n)]
    return list(reversed(out))

def run_backtest():
    parts = []
    for fold, (ts, te) in enumerate(fold_windows(feat_df["date"].max(), N_FOLDS, FOLD_H)):
        train = feat_df[feat_df["date"] < ts]
        test = feat_df[(feat_df["date"] >= ts) & (feat_df["date"] <= te)].copy()
        cut = train["date"].max() - pd.Timedelta(days=FOLD_H)
        tr2, va2 = train[train["date"] <= cut], train[train["date"] > cut]
        qm = QuantileGBM(BT_QUANTILES, LGB_PARAMS).fit(tr2[FEATURES], tr2["demand"], va2[FEATURES], va2["demand"])
        mp = dict(LGB_PARAMS); mp["objective"] = "regression"
        mm = lgb.LGBMRegressor(random_state=7, **mp)
        mm.fit(tr2[FEATURES], tr2["demand"], eval_set=[(va2[FEATURES], va2["demand"])], eval_metric="l2",
               callbacks=[lgb.early_stopping(50, verbose=False), lgb.log_evaluation(0)])
        P = qm.predict(test[FEATURES])
        frac = fractile.reindex(test["StockCode"]).to_numpy()
        test["Q_newsvendor"] = order_at_fractile(P, frac, BT_QUANTILES)
        test["Q_median"] = P[0.5].to_numpy()
        test["Q_point"] = np.clip(mm.predict(test[FEATURES]), 0, None)
        test["cu"] = cu_prod.reindex(test["StockCode"]).to_numpy(); test["fold"] = fold
        parts.append(test[["fold","date","StockCode","demand","cu","Q_point","Q_median","Q_newsvendor"]])
    return pd.concat(parts, ignore_index=True)

bt = run_backtest()
D, CU = bt["demand"].to_numpy(float), bt["cu"].to_numpy(float)
POL = {"point (mean)":"Q_point", "median (q0.5)":"Q_median", "newsvendor":"Q_newsvendor"}
for name, qc in POL.items():
    Q = bt[qc].to_numpy(float)
    bt[qc.replace("Q_","cost_")] = CU*np.maximum(D-Q,0) + HOLDING*np.maximum(Q-D,0)
print(f"backtest: {N_FOLDS} folds x {FOLD_H}d | {len(bt):,} product-days scored")

backtest: 3 folds x 28d | 16,800 product-days scored

The headline — cumulative realized cost

Every product-day in the test windows adds its realized £ cost. The policy whose line stays lowest wins.

COST = {n: qc.replace("Q_","cost_") for n, qc in POL.items()}
COL = {"point (mean)":ORANGE, "median (q0.5)":PURP, "newsvendor":GREEN}
daily = bt.groupby("date")[list(COST.values())].sum()
cum = daily.cumsum()
fig, ax = plt.subplots(figsize=(12, 5))
for n in POL: ax.plot(cum.index, cum[COST[n]], lw=2, color=COL[n], label=n)
ax.set_ylabel("cumulative realized cost (£)"); ax.set_title("Realized cost over the backtest windows")
ax.legend(fontsize=9); fig.tight_layout(); plt.show()

tot = {n: bt[COST[n]].sum() for n in POL}
voi = (tot["point (mean)"] - tot["newsvendor"]) / tot["point (mean)"]
print(f"total realized cost:  point £{tot['point (mean)']:,.0f} | median £{tot['median (q0.5)']:,.0f} | "
      f"newsvendor £{tot['newsvendor']:,.0f}")
print(f"value of modeling uncertainty: newsvendor cuts realized cost {voi:+.1%} vs the point-forecast policy")

total realized cost:  point £327,126 | median £359,802 | newsvendor £273,949
value of modeling uncertainty: newsvendor cuts realized cost +16.3% vs the point-forecast policy

Why newsvendor wins — the service/holding trade-off

The newsvendor policy deliberately holds more stock to buy a much higher fill rate; because a stockout costs far more than a unit of holding (per-product $C_u\gg C_o$), that trade is a net win. The median policy shows the opposite failure — the median of lumpy demand under-stocks.

def fill(Q): return np.minimum(Q, D).sum() / D.sum()
order = ["point (mean)", "median (q0.5)", "newsvendor"]; cols = [COL[n] for n in order]
fig, ax = plt.subplots(1, 3, figsize=(14, 4.2))
ax[0].bar(order, [tot[n] for n in order], color=cols); ax[0].set_title("Total realized cost (£)"); ax[0].tick_params(axis="x", rotation=12)
for i, n in enumerate(order): ax[0].text(i, tot[n], f"£{tot[n]:,.0f}", ha="center", va="bottom", fontsize=8)
ax[1].bar(order, [fill(bt[POL[n]].to_numpy(float)) for n in order], color=cols)
ax[1].set_title("Fill rate (demand met)"); ax[1].set_ylim(0, 1); ax[1].tick_params(axis="x", rotation=12)
x = np.arange(len(order)); w = 0.38
held = [np.maximum(bt[POL[n]].to_numpy(float)-D, 0).sum() for n in order]
short = [np.maximum(D-bt[POL[n]].to_numpy(float), 0).sum() for n in order]
ax[2].bar(x-w/2, held, w, label="held (overage)", color=GREY); ax[2].bar(x+w/2, short, w, label="short (underage)", color=RED)
ax[2].set_xticks(x); ax[2].set_xticklabels(order, rotation=12); ax[2].set_title("Units held vs short"); ax[2].legend(fontsize=8)
fig.tight_layout(); plt.show()

Robustness across folds

The advantage should hold in every rolling window, not just on average.

pf = bt.groupby("fold")[[COST[n] for n in order]].sum()
x = np.arange(len(pf)); w = 0.26
fig, ax = plt.subplots(figsize=(9, 4.3))
for i, n in enumerate(order): ax.bar(x + (i-1)*w, pf[COST[n]], w, label=n, color=COL[n])
ax.set_xticks(x); ax.set_xticklabels([f"fold {i}" for i in pf.index]); ax.set_ylabel("realized cost (£)")
ax.set_title("Realized cost by fold (rolling origin)"); ax.legend(fontsize=8); fig.tight_layout(); plt.show()
wins = int((pf[COST["newsvendor"]] < pf[[COST["point (mean)"], COST["median (q0.5)"]]].min(axis=1)).sum())
print(f"newsvendor is the cheapest policy in {wins}/{len(pf)} folds")

newsvendor is the cheapest policy in 3/3 folds

Section 7 takeaways

• On held-out data, the distribution-driven newsvendor policy cuts realized inventory cost versus an equally-good point forecaster — the dollar value of modeling uncertainty, measured in £.
• The win is robust across folds and mechanistically clear (higher fill rate, far fewer costly stockouts), not a fluke of one window.
• Predict-then-optimize, validated end to end: distribution → decision → measured £ saved.

Next — Section 8: wrap-up, tying the thread from lumpy demand to a calibrated distribution to a per-product decision to money saved.

8 · Wrap-up — from lumpy demand to money saved

One thread runs through the notebook:

1. The demand is lumpy (§3): ~94% of the top products sell often but in wildly variable sizes, with ~48% zero-days. A point forecast is the wrong object.
2. So we forecast the distribution (§5): quantile gradient boosting, calibrated exactly in the upper tail where the decision lives (the flat low end is the honest zero-atom, not a defect).
3. The distribution drives a decision (§6): the newsvendor order at each product's price-driven critical fractile — an order sheet a buyer can execute, holding more safety stock for pricier items.
4. And it saves money (§7): on held-out data that policy cut realized cost 16.3% versus an equally-good point forecaster, in every fold.

Same features, same model family — the gain came from using the distribution and the economics, not from a cleverer forecast. That is predict-then-optimize: the pinball loss the model minimized is the newsvendor cost at the critical fractile, so predicting and optimizing are a single objective.

What carried the analysis: the leak-safe shift(H) discipline (§4) that makes the backtest trustworthy; per-product economics from real prices (§6), richer than a fixed service level; and honest calibration plus a SHAP audit (§5) — earning trust before the forecast drives a decision.

Alternative approaches — and when they would win

Quantile GBM is the pragmatic choice here (it won the real M5 competition), but not the only one — and several sit squarely in this portfolio's other arcs:

• Bayesian hierarchical hurdle / zero-inflated NegBin — the statistician's answer, and arguably a more principled match to intermittent demand: a hurdle model factorizes demand into P(any sale) × size | sale — exactly Croston's decomposition — while product-level random effects pool information the way the global GBM does. It returns a coherent posterior predictive distribution (full uncertainty, not independently-fit quantiles) and handles counts/zeros natively. It would likely win when data is scarce (few products, short history, cold-start items), where priors and pooling matter most; it loses on scale (MCMC over 200×739 vs seconds for GBM) and on the heavy bulk-order tail, which no NegBin shape captures. (Builds directly on the ZIP/ZINB and hierarchical-Poisson/NegBin work in the Bayesian arc.)
• Croston / SBA — the classical intermittent-demand baseline (the method behind §3's classification). Cheap and robust, but it yields a demand rate, not a distribution — you must bolt on a distributional assumption to reach a newsvendor quantile. Best as a sanity benchmark.
• DeepAR-style probabilistic RNN — the modern deep-learning competitor (a global recurrent net emitting a per-step distribution). It can shine on long, many-series panels, but on short daily retail series GBMs usually win (as M5 showed) — consistent with the portfolio's finding that deep learning beats classical on structured, not tabular/short, data.
• Gaussian-process / BNP density regression — beautifully calibrated uncertainty, but the $O(n^3)$ cost is prohibitive at 147k rows and the smooth-latent assumption fights lumpy spikes.

Verdict: for this problem — many series, ample history, a newsvendor that needs only a marginal tail quantile — quantile GBM is hard to beat, because it targets that exact quantile cheaply. The Bayesian hurdle model is the one worth racing it against, and the natural hybrid is to let the Bayesian model handle sparse/new products (pooling + priors) and GBM handle the data-rich tail.

The bigger picture

This notebook is the research telling of a predict-then-optimize pipeline that also exists as a tested, deployable Python package — the same idea, engineered. Together they make the point a portfolio should: the method is sound and it ships.