Demand Forecasting in Python: Baselines, Honest Backtests and the Leaks That Fake Accuracy
Most demand forecasting notebooks end with a number that is too good. The model is fine; the evaluation is not. A rolling mean that includes today, a lag_1 feature you will never have two weeks ahead, a promo flag that w
Most demand forecasting notebooks end with a number that is too good. The model is fine; the evaluation is not. A rolling mean that includes today, a lag_1 feature you will never have two weeks ahead, a promo flag that was filled in after the fact. Each one shaves a few points off the error, and the gap shows up later as empty shelves or a warehouse full of slow movers.
This is a hands-on guide to demand forecasting in Python for a supply chain setting. We will build a forecasting pipeline on a deterministic synthetic sales dataset, measure it the way a planner would, break it on purpose with three common data leaks, and then turn quantile forecasts into safety stock.
What we will build:
- A reproducible sales history: 20 SKUs Γ 5 stores, three years of daily data with weekly and yearly seasonality, promos, holidays, persistent demand shocks and intermittent demand.
- Metrics that survive zero-sales days: MAE, WAPE and bias.
- A rolling-origin backtest with a 14-day horizon.
- Baselines: naive, seasonal naive, 28-day moving average.
- A LightGBM model with lag features that respect the horizon.
- Three leaks, a βbadβ evaluation next to an honest one, and pytest checks that catch them.
- P50/P90 quantile forecasts and three ways to size safety stock.
Stack: Python 3.13, numpy 2.5.3, pandas 3.0.5, scikit-learn 1.9.1 (needed by LightGBMβs scikit-learn API), LightGBM 4.7.0, pytest 9.1.1. Every table below is real output from one run of the code, and the whole run is deterministic: rerun it and you get the same digits.
Why gradient boosting, and why baselines anyway
The M5 competition (Walmart unit sales, 42,840 hierarchical series, 28-day horizon) is the closest public benchmark to retail demand forecasting. Its organizers report that LightGBM was used by all of the top 50 competitors, and that the winning entry was an equal-weighted combination of LightGBM models. The same paper notes that simple exponential smoothing remained competitive at the product-store level. So gradient boosting is a reasonable default among demand forecasting models, but it has to earn its place against simple methods. Hyndman and Athanasopoulos put it bluntly in Forecasting: Principles and Practice: if a method does not beat the simple benchmarks, it is not worth considering.
Here is the setup we will use:
Step 1. Generate a sales history you can reason about
Real data would be better for the results, but worse for learning: you cannot know which patterns are real. A synthetic dataset lets us plant every effect on purpose and then check whether the model finds it. Three segments matter for inventory: fast movers, slow movers, and intermittent SKUs where most days have zero sales.
The generator multiplies a base level by a weekday profile (Friday/Saturday peak), a yearly sine wave, a mild trend, U.S. retail holidays (Black Friday spike, Christmas ramp, stores closed on December 25 and January 1), planned chain-wide promos, ad-hoc store markdowns, and an AR(1) shock on the log scale that makes good and bad weeks persist. Sales are drawn from a gamma-Poisson mixture, which gives overdispersed counts like real retail data.
Spoiler on dev.to: forecast/data.py (full generator)
"""Deterministic synthetic retail sales: SKUs x stores, daily, 2023-2025."""
import numpy as np
import pandas as pd
START, END = "2023-01-01", "2025-12-31"
SEGMENTS = {"fast": (20.0, 60.0), "slow": (2.0, 8.0), "intermittent": (0.15, 0.5)}
def us_holidays(years):
"""Dates with a demand spike (or a closed store) that are known a year ahead."""
days = {}
for y in years:
nov = pd.date_range(f"{y}-11-01", f"{y}-11-30")
thanksgiving = nov[nov.dayofweek == 3][3]
days[pd.Timestamp(f"{y}-07-04")] = 1.4
days[thanksgiving - pd.Timedelta(days=1)] = 1.6
days[thanksgiving + pd.Timedelta(days=1)] = 2.2 # Black Friday
for d in range(18, 25):
days[pd.Timestamp(f"{y}-12-{d}")] = 1.3 + 0.1 * (d - 18)
days[pd.Timestamp(f"{y}-12-25")] = 0.0 # stores closed
days[pd.Timestamp(f"{y}-01-01")] = 0.0
return days
def make_sales(n_skus=20, n_stores=5, seed=42):
rng = np.random.default_rng(seed)
dates = pd.date_range(START, END, freq="D")
n = len(dates)
holidays = us_holidays(sorted(set(dates.year)))
hol_mult = np.array([holidays.get(d, 1.0) for d in dates])
dow = dates.dayofweek.to_numpy()
doy = dates.dayofyear.to_numpy()
seg_names = list(SEGMENTS)
frames = []
for sku in range(n_skus):
segment = seg_names[sku % 3]
lo, hi = SEGMENTS[segment]
base = rng.uniform(lo, hi)
weekly = 1 + rng.uniform(0.1, 0.4) * np.isin(np.arange(7), [4, 5]) # Fri/Sat peak
weekly = weekly / weekly.mean()
amp, phase = rng.uniform(0.0, 0.5), rng.uniform(0, 2 * np.pi)
yearly = 1 + amp * np.sin(2 * np.pi * doy / 365.25 + phase)
trend = 1 + rng.uniform(-0.1, 0.2) * np.arange(n) / n
# Planned, chain-wide promos: a 7-day promo every 6-10 weeks, known in advance.
promo_planned = np.zeros(n, dtype=bool)
t = int(rng.integers(10, 60))
while t < n:
promo_planned[t:t + 7] = True
t += int(rng.integers(42, 71))
lift = rng.uniform(1.5, 2.0)
for store in range(n_stores):
store_mult = rng.uniform(0.6, 1.4)
# Ad-hoc store markdowns: decided on the spot, logged after the fact.
adhoc = np.zeros(n, dtype=bool)
for start in np.flatnonzero(rng.random(n) < 0.015):
adhoc[start:start + 3] = True
# Local demand shocks that persist for weeks (AR(1) on the log scale).
shock, eps = np.zeros(n), rng.normal(0, 0.08, n)
for i in range(1, n):
shock[i] = 0.95 * shock[i - 1] + eps[i]
mu = (base * store_mult * weekly[dow] * yearly * trend * hol_mult
* np.where(promo_planned, lift, 1.0) * np.where(adhoc, 1.8, 1.0) * np.exp(shock))
# Gamma-Poisson mixture = negative binomial noise (overdispersed counts).
y = rng.poisson(mu * rng.gamma(6.0, 1 / 6.0, n))
frames.append(pd.DataFrame({
"date": dates, "sku": f"SKU{sku:02d}", "store": f"S{store}",
"segment": segment, "sales": y.astype(float),
"promo_planned": promo_planned.astype(int),
"promo_actual": (promo_planned | adhoc).astype(int),
"holiday": (hol_mult != 1.0).astype(int),
}))
df = pd.concat(frames, ignore_index=True)
df["series"] = df["sku"] + "_" + df["store"]
return df.sort_values(["series", "date"], ignore_index=True)
Two details matter later. promo_planned is the promo calendar that exists weeks in advance. promo_actual also includes ad-hoc markdowns that a store manager decided on the day and that only appear in the log afterwards.
109,600 rows, 100 series, 2023-01-01 .. 2025-12-31
series mean_daily zero_days
segment
fast 35 61.26 1%
intermittent 30 0.40 70%
slow 35 6.49 5%
70% of days for intermittent SKUs have zero sales. Keep that number in mind when we pick metrics.
Step 2. Pick metrics that survive zeros
MAPE is the metric people ask for and the one you should not use here. Forecasting: Principles and Practice points out that percentage errors are infinite or undefined when an actual value is zero, and our intermittent SKUs are zero most days. We use three metrics instead, plus pinball loss for quantiles:
# ---------- metrics -------------------------------------------------------
def mae(y, f):
return float(np.mean(np.abs(f - y)))
def wape(y, f):
"""Sum of absolute errors / sum of actuals. Safe with zero-demand days."""
return float(np.sum(np.abs(f - y)) / np.sum(y))
def bias(y, f):
"""> 0 means over-forecasting (excess stock), < 0 means stockouts."""
return float(np.sum(f - y) / np.sum(y))
def pinball(y, q, alpha):
d = y - q
return float(np.mean(np.maximum(alpha * d, (alpha - 1) * d)))
- MAE is in units, easy to explain, but not comparable across SKUs of different volume.
- WAPE divides total absolute error by total demand, so it works across SKUs and is defined as long as the sum of demand is positive. It is a volume-weighted MAE, which means fast movers dominate it.
- Bias is the signed version. It tells you whether errors turn into overstock (positive) or stockouts (negative). Two models with the same WAPE can have very different bias, and inventory cares a lot about the sign.
- Pinball loss scores a quantile forecast. At P90, an under-forecast costs nine times more than an over-forecast of the same size. (FPP3 defines the quantile score with an extra factor of 2; the ranking is identical.)
Step 3. Backtest on a rolling origin
A single train/test split gives you one number from one season. A rolling-origin backtest (FPP3 calls it βevaluation on a rolling forecasting originβ) repeats the experiment: train on everything up to an origin, forecast the next 14 days, move the origin forward, repeat.
# ---------- rolling origin -------------------------------------------------
def origins(df, n_folds=6, step=28, horizon=HORIZON):
last = df["date"].max()
return [last - pd.Timedelta(days=horizon + step * k) for k in range(n_folds - 1, -1, -1)]
def split(df, origin, horizon=HORIZON):
train = df["date"] <= origin
test = (df["date"] > origin) & (df["date"] <= origin + pd.Timedelta(days=horizon))
return train.to_numpy(), test.to_numpy()
fold 0: train <= 2025-07-30 test 2025-07-31 .. 2025-08-13
fold 1: train <= 2025-08-27 test 2025-08-28 .. 2025-09-10
fold 2: train <= 2025-09-24 test 2025-09-25 .. 2025-10-08
fold 3: train <= 2025-10-22 test 2025-10-23 .. 2025-11-05
fold 4: train <= 2025-11-19 test 2025-11-20 .. 2025-12-03
fold 5: train <= 2025-12-17 test 2025-12-18 .. 2025-12-31
Rolling-origin backtest: six folds, each trains on all data up to its origin and is scored on the next 14 days
The last fold covers the Christmas ramp and the December 25 closure, which is exactly the kind of period a single random split would average away.
Step 4. Baselines first
Each baseline only sees history up to the origin:
# ---------- baselines (use only history up to the origin) ------------------
def baselines(df, origin, horizon=HORIZON):
hist = df[df["date"] <= origin]
test = df[(df["date"] > origin) & (df["date"] <= origin + pd.Timedelta(days=horizon))]
last = hist.groupby("series")["sales"].last()
ma28 = hist.groupby("series")["sales"].apply(lambda s: s.iloc[-28:].mean())
last_week = hist[hist["date"] > origin - pd.Timedelta(days=7)]
snaive = last_week.set_index(["series", last_week["date"].dt.dayofweek])["sales"]
key = list(zip(test["series"], test["date"].dt.dayofweek))
return pd.DataFrame({
"naive": test["series"].map(last).to_numpy(),
"seasonal_naive": snaive.reindex(key).to_numpy(),
"moving_avg_28": test["series"].map(ma28).to_numpy(),
}, index=test.index)
The naive forecast repeats the last observed day. Seasonal naive repeats the same weekday from the last week before the origin. The moving average uses the last 28 days. We will see their scores next to the model in the next step.
Step 5. LightGBM with lag features that respect the horizon
The single most important rule for time series forecasting in Python with tabular models: a feature is only valid if you would have it at the forecast origin. We forecast days t+1 to t+14 from origin t with one model, so the freshest sales value we can use for day t+14 is from day t. Every sales-derived feature is therefore shifted by at least 14 days:
"""Feature sets. HONEST uses only what is known at the forecast origin."""
import pandas as pd
HORIZON = 14 # we forecast days t+1 .. t+14 from origin t
def calendar(df):
out = pd.DataFrame(index=df.index)
out["dow"] = df["date"].dt.dayofweek
out["doy"] = df["date"].dt.dayofyear
out["month"] = df["date"].dt.month
xmas = pd.to_datetime(df["date"].dt.year.astype(str) + "-12-25")
out["days_to_xmas"] = (xmas - df["date"]).dt.days.clip(-10, 60)
out["holiday"] = df["holiday"]
out["promo"] = df["promo_planned"]
out["sku"] = df["sku"].astype("category")
out["store"] = df["store"].astype("category")
return out
def honest_features(df, horizon=HORIZON):
"""Every sales-derived feature is shifted by >= horizon days."""
g = df.groupby("series")["sales"]
out = calendar(df)
for lag in (horizon, horizon + 7, horizon + 14, horizon + 21, 364):
out[f"lag_{lag}"] = g.shift(lag)
past = g.shift(horizon) # the last value we can see for the furthest day
pg = past.groupby(df["series"])
for w in (7, 28):
out[f"roll_mean_{w}"] = pg.transform(lambda s: s.rolling(w).mean())
out["roll_std_28"] = pg.transform(lambda s: s.rolling(28).std())
out["zero_share_28"] = (past == 0).astype(float).groupby(df["series"]).transform(
lambda s: s.rolling(28).mean())
return out
def leaky_features(df, leaks=("short_lags", "unshifted_roll", "promo_actual")):
"""The same table with the three mistakes we see most often in forecasting code."""
out = honest_features(df)
g = df.groupby("series")["sales"]
if "short_lags" in leaks: # lag_1 is not known 14 days ahead
out["lag_1"] = g.shift(1)
out["lag_7"] = g.shift(7)
if "unshifted_roll" in leaks: # rolling window includes the target day itself
out["roll_mean_7_now"] = g.transform(lambda s: s.rolling(7).mean())
if "promo_actual" in leaks: # promo log includes ad-hoc markdowns made later
out["promo"] = df["promo_actual"]
return out
Calendar features (weekday, day of year, days to Christmas, holiday flag) and the planned promo calendar are known in advance, so they are safe without a shift. zero_share_28 gives the model a direct signal for intermittent series.
Training uses the Tweedie objective, which LightGBM documents as βTweedie regression with log-linkβ. It suits non-negative, zero-heavy targets; the M5 winner also trained its LightGBM models on a Tweedie likelihood for that reason. We pin n_jobs and deterministic=True so reruns are byte-identical:
# ---------- gradient boosting ---------------------------------------------
PARAMS = dict(n_estimators=300, learning_rate=0.05, num_leaves=31, min_child_samples=50,
subsample=0.8, subsample_freq=1, colsample_bytree=0.8,
random_state=7, n_jobs=1, deterministic=True, force_row_wise=True, verbose=-1)
def fit_predict(X, y, train, test, objective="tweedie", alpha=None):
params = dict(PARAMS, objective=objective)
if alpha is not None:
params["alpha"] = alpha
model = lgb.LGBMRegressor(**params)
model.fit(X[train], y[train])
return np.clip(model.predict(X[test]), 0, None)
def rolling_backtest(df, X, **kw):
"""Train on everything up to each origin, predict the next HORIZON days."""
y = df["sales"].to_numpy()
parts = []
for i, origin in enumerate(origins(df)):
train, test = split(df, origin)
parts.append(pd.DataFrame({"fold": i, "pred": fit_predict(X, y, train, test, **kw)},
index=df.index[test]))
return pd.concat(parts)
Results over all six folds (100 series Γ 14 days Γ 6 folds = 8,400 SKU-store-days):
MAE WAPE bias
method
naive 14.08 54.4% -18.1%
seasonal_naive 15.25 58.9% -7.4%
moving_avg_28 12.69 49.0% -5.8%
lightgbm_honest 10.90 42.1% -4.7%
LightGBM cuts WAPE from 49.0% (the best baseline) to 42.1% and has the smallest bias. Two things are worth noticing in the baselines. Seasonal naive loses to the plain moving average: one week of history is a noisy sample when daily sales are overdispersed. And the naive forecast has an 18% negative bias: all six origins fall on a Wednesday (they are 28 days apart), a below-average weekday, so βrepeat the last dayβ under-forecasts every Friday and Saturday peak.
WAPE by segment tells a less flattering story:
segment fast intermittent slow
method
lightgbm_honest 40.9 124.6 48.8
moving_avg_28 47.9 127.7 55.3
seasonal_naive 57.5 143.8 68.2
WAPE by method across six rolling-origin folds; the leaky evaluations look better than any honest model
Every method scores above 100% WAPE on intermittent SKUs. A forecast of zero for every day scores exactly 100% by definition (our test suite checks this), so on these series a point forecast is worse than predicting nothing. That is not a LightGBM bug. With 70% zero days, any forecast of a positive mean is wrong almost every day. For intermittent demand, the useful output is a quantile or a probability of selling at least one unit over the lead time, not a daily point forecast. Crostonβs method is the classic alternative, but FPP3 notes its forecasts are biased and it has no algebraic prediction intervals.
Step 6. Three leaks that make the model look better than it is
Now letβs break it. Leakage is not a beginner problem: Kapoor and Narayanan found it in 329 papers across 17 scientific fields that use machine learning. In forecasting code, three versions show up again and again:
- Short lags. lag_1 and lag_7 are great features for tomorrow. For a forecast made 14 days ahead, they reference sales that have not happened yet.
- Unshifted rolling windows. In pandas, a fixed window uses closed='right' by default, so s.rolling(7).mean() includes the current row. The target leaks into its own feature at 1/7 weight.
- Promo actuals instead of the promo plan. The promo log includes ad-hoc markdowns recorded after the fact. At the forecast origin you only have the plan.
The leaky_features function at the bottom of features.py adds all three. Then we evaluate the way many notebooks do, with a random 80/20 split on 2024β2025 rows, and compare it with the rolling-origin backtest:
# ---- 4. leakage: the 'too good to be true' evaluation ----------------------
X_leaky = leaky_features(df)
rng = np.random.default_rng(0)
recent = (df["date"] >= "2024-01-01").to_numpy()
idx = np.flatnonzero(recent)
test_rand = np.zeros(len(df), bool)
test_rand[rng.choice(idx, size=len(idx) // 5, replace=False)] = True
train_rand = recent & ~test_rand
yall = df["sales"].to_numpy()
leak_rows = [
row("random split + leaky features", yall[test_rand],
fit_predict(X_leaky, yall, train_rand, test_rand)),
row("rolling origin + leaky features", y,
rolling_backtest(df, X_leaky).loc[ev.index, "pred"].to_numpy()),
row("rolling origin + honest features", y, ev["pred"].to_numpy()),
]
report("05_leaky_vs_honest", fmt(leak_rows))
MAE WAPE bias
method
random split + leaky features 8.40 35.2% -0.8%
rolling origin + leaky features 9.21 35.6% -3.3%
rolling origin + honest features 10.90 42.1% -4.7%
The leaky setup reports 35.2% WAPE and a bias of -0.8%. The honest backtest says 42.1% and -4.7%. That is 7 points of WAPE you would promise to the business and never deliver, and a model that βlooks unbiasedβ while it systematically under-forecasts. The random split also hides the bias: with neighboring days on both sides of the split, the model is never asked to extrapolate into a December it has not seen.
Which leak does the damage? We add one at a time to the honest feature set:
MAE WAPE bias
method
honest + short_lags 10.49 40.5% -4.2%
honest + unshifted_roll 9.50 36.7% -3.4%
honest + promo_actual 10.72 41.4% -4.0%
Feature availability at the forecast origin: which features are known for days t+1 to t+14
The unshifted rolling mean is the worst offender (42.1% β 36.7%), because it contains the target itself. Short lags help less than I expected (40.5%): with a 14-day horizon, lag_1 is only βfutureβ for some of the days, and the demand shocks in this data change slowly, so lag_14 and the 7-day mean behind it probably carry most of what lag_1 would add. The promo leak is small here (41.4%) because ad-hoc markdowns cover only about 4% of days; in a business with frequent unplanned discounts it would be larger.
You should not rely on code review to catch these. The test below hides all sales after the origin and recomputes the features. Honest features for the test window must not change; leaky ones must:
def hide_future(df, origin):
"""What the warehouse actually knows at the origin: no sales after it."""
masked = df.copy()
masked.loc[masked["date"] > origin, "sales"] = np.nan
return masked
def test_data_is_deterministic(df):
again = make_sales(n_skus=6, n_stores=2)
pd.testing.assert_frame_equal(df, again)
assert not make_sales(n_skus=6, n_stores=2, seed=1)["sales"].equals(df["sales"])
@pytest.mark.parametrize("origin_idx", range(6))
def test_honest_features_do_not_change_when_future_is_hidden(df, origin_idx):
origin = origins(df)[origin_idx]
_, test = split(df, origin)
full = honest_features(df)[test]
blind = honest_features(hide_future(df, origin))[test]
pd.testing.assert_frame_equal(full, blind)
@pytest.mark.parametrize("leak", ["short_lags", "unshifted_roll"])
def test_leaky_features_do_change(df, leak):
origin = origins(df)[0]
_, test = split(df, origin)
full = leaky_features(df, leaks=(leak,))[test]
blind = leaky_features(hide_future(df, origin), leaks=(leak,))[test]
with pytest.raises(AssertionError):
pd.testing.assert_frame_equal(full, blind)
............ [100%]
12 passed in 1.61s
This βhide the future and recomputeβ test is cheap, works for any feature pipeline, and fails loudly the day someone adds rolling(7) without shift. It does not catch the promo leak, because that leak lives in the data source rather than in the code. For that you need a separate check that compares the feature with what the planning system knew at the origin.
Step 7. Quantile forecasts for safety stock
A mean forecast tells you how much to expect. Safety stock depends on how much you might need. LightGBM supports this directly with objective="quantile" and alpha, so we train P50 and P90 models on the honest features with the same rolling-origin folds:
# ---- 5. quantiles + safety stock -------------------------------------------
q = {a: rolling_backtest(df, X_honest, objective="quantile", alpha=a).loc[ev.index, "pred"]
for a in (0.5, 0.9)}
ev["p50"], ev["p90"] = q[0.5].to_numpy(), q[0.9].to_numpy()
qrows = []
for s, g in ev.groupby("segment"):
qrows.append({"segment": s,
"P90 coverage": f"{(g['sales'] <= g['p90']).mean():.1%}",
"P90 pinball": f"{pinball(g['sales'], g['p90'], 0.9):.2f}",
"P50 bias": f"{bias(g['sales'], g['p50']):+.1%}",
"mean-model bias": f"{bias(g['sales'], g['pred']):+.1%}"})
report("07_quantiles", pd.DataFrame(qrows).set_index("segment").to_string())
RESULTS["quantiles"] = qrows
P90 coverage P90 pinball P50 bias mean-model bias
segment
fast 88.2% 8.43 -15.6% -4.8%
intermittent 90.9% 0.16 -73.1% -1.9%
slow 88.3% 1.10 -15.4% -4.5%
P90 coverage lands at 88β91% on the daily level, which is close to the target. The P50 column is the trap: the median model under-forecasts by 15% for fast and slow movers and by 73% for intermittent SKUs. For right-skewed demand the median is below the mean, and for an intermittent SKU the median day is zero. If you train a βMAE-optimalβ model and feed it into replenishment as the expected demand, you get stockouts. Use the mean model (Tweedie or L2) for expected demand and explicit quantiles for buffers.
Inventory is not ordered per day, though. With a 7-day lead time, the question is: what reorder level covers 90% of 7-day demand? We compare three answers on 7-day windows, calibrating on folds 0β2 and testing on folds 3β5:
# Lead time L = 7 days, target cycle service level 90%.
# Each fold gives two 7-day windows per series (days 1-7 and 8-14 after the origin).
ev["fold"] = bt.loc[ev.index, "fold"].to_numpy()
ev["day"] = (ev["date"] - ev["fold"].map(dict(enumerate(folds)))).dt.days
ev["window"] = (ev["day"] - 1) // 7
ev["err"] = ev["sales"] - ev["pred"]
lt = ev.groupby(["series", "segment", "fold", "window"]).agg(
demand=("sales", "sum"), mean_fc=("pred", "sum"), p90_sum=("p90", "sum")).reset_index()
calib, hold = lt[lt["fold"] < 3], lt[lt["fold"] >= 3]
sigma = ev[ev["fold"] < 3].groupby("series")["err"].std()
z90 = 1.2816
resid_q = (calib["demand"] - calib["mean_fc"]).groupby(calib["segment"]).quantile(0.9)
hold = hold.assign(
normal=hold["mean_fc"] + z90 * hold["series"].map(sigma) * np.sqrt(7),
sum_p90=hold["p90_sum"],
empirical=hold["mean_fc"] + hold["segment"].map(resid_q).clip(lower=0),
)
ss_rows = []
for s, g in hold.groupby("segment"):
for m in ("normal", "sum_p90", "empirical"):
ss_rows.append({"segment": s, "method": m,
"coverage": float((g["demand"] <= g[m]).mean()),
"safety_stock": float((g[m] - g["mean_fc"]).mean())})
- normal is the textbook formula: expected lead-time demand + z Γ Ο_daily Γ βL. It assumes daily errors are independent and normal.
- sum_p90 adds up the seven daily P90 forecasts.
- empirical adds the 90th percentile of past 7-day residuals (from the calibration folds) to the 7-day mean forecast.
cycle service level (target 90%)
segment fast intermittent slow
method
empirical 84% 90% 91%
normal 74% 87% 80%
sum_p90 95% 100% 97%
mean safety stock, units per 7-day window
segment fast intermittent slow
method
empirical 189.4 3.1 25.0
normal 104.3 2.4 16.3
sum_p90 351.6 7.2 39.9
std of 7-day error / (sqrt(7) * std of daily error), median per series
segment
fast 1.46
intermittent 0.94
slow 1.25
Cycle service level versus mean safety stock for fast movers: normal formula under-covers, summed P90 overstocks
The normal formula delivers 74% instead of 90% on fast movers. The last block explains why: the standard deviation of 7-day errors is 1.46 times larger than β7 Γ daily Ο, because forecast errors on consecutive days are correlated (a bad week stays bad). The βL rule only holds when they are independent, an assumption the Wikipedia article on safety stock lists explicitly along with the warning that real demand often violates it.
Summing daily P90s fails in the other direction: 95β100% coverage with roughly twice the safety stock of the empirical method on fast movers (351.6 vs 189.4 units). Quantiles do not add up. Unless all seven days move together perfectly, the sum of seven daily P90s sits well above the P90 of the weekly total.
The empirical lead-time quantile lands closest to target: 84% / 90% / 91%. It is not perfect for fast movers, since the holdout folds include the Black Friday and Christmas peaks, which were not in the calibration folds. With 180β210 test windows per segment, these coverage numbers also carry a few points of sampling noise. In production, you would recalibrate the residual quantile every fold and per segment or SKU cluster, or train a model directly on 7-day aggregated demand.
Gotchas I hit along the way
- TimeSeriesSplit on a stacked panel. scikit-learnβs TimeSeriesSplit splits by row position and says samples must be equally spaced. Our DataFrame is sorted by series, then date, so it would train on the first SKUs and test on the last ones. Split by date, as split() does above.
- WAPE above 100% is not a bug. It means the forecast is worse than predicting zero. Report WAPE per segment, or a single WAPE will be dominated by fast movers and hide the problem.
- Bias is the metric planners feel. A -4.7% bias across all SKUs is the difference between a model that runs out on promo weeks and one that does not. Track it per fold, not just overall.
- Thread oversubscription. On a shared 2-vCPU sandbox with other jobs running, one LightGBM fit took 60 s with n_jobs=2 and 4.8 s with n_jobs=1. LightGBMβs docs recommend setting threads to the number of real CPU cores; on shared CI runners, that number is often lower than os.cpu_count() says.
- One model, fixed lags. Shifting every lag by 14 days throws away fresh information for day t+1. A direct approach (one model per horizon bucket, each with its own minimum lag) usually gains a few points on short horizons. The leakage test above works the same way for each bucket.
Production checklist for demand forecasting in supply chain
β Every feature has a documented availability time, and the pipeline is tested with future data masked.
β Promo, price and assortment features come from the plan as of the origin, not from the actuals table.
β Rolling-origin backtest with at least one fold covering a peak season.
β Naive, seasonal naive and moving average baselines in every report.
β WAPE and bias reported per segment (fast, slow, intermittent) and per fold.
β Expected demand from a mean-type objective; buffers from quantiles or empirical lead-time residuals.
β Safety stock evaluated on achieved service level, not on forecast error.
β Deterministic training (fixed seeds, pinned versions) so a changed number means a changed model.
Wrapping up
Comparing demand forecasting methods is mostly an evaluation problem. On the same data, the same LightGBM model scores 35.2% or 42.1% WAPE depending on whether the features can see the future, and the safety stock built on top of it covers 74% or 90% of lead-time demand depending on how you turn errors into buffers. The honest pipeline is about 200 lines of Python plus one pytest file, and it gives the planning team a number they can hold you to.
Which leak have you found in a forecasting pipeline that looked great in the notebook? I would like to hear about the less obvious ones, especially in retail demand forecasting with promotions.
Originally published by Dev.to WebDev. Aggregated on AIWithGhost for educational purposes β full credit and traffic to the original publisher.