Design a demand forecasting system

The problem

A grocery delivery company runs 4,000 dark stores in 40 cities. Every night it must decide how much of each product to order for each store, and every week how many pickers and riders to schedule per store and hour. Order too little and customers see "out of stock" or wait 40 minutes; order too much and fresh food is thrown away and staff stand idle.

All of those decisions rest on one system: a demand forecast per store, product and day (or hour), weeks ahead. It must be accurate where data is thin (a new store, a slow-selling product), consistent across levels (the city plan must equal the sum of its stores), aware of holidays and promotions, and honest about uncertainty, because the order quantity depends on how wrong the forecast might be, not only on its centre.

Think of a bakery owner planning tomorrow's bread. They know Saturdays are busy, the street festival doubles sales, and a rainy Tuesday is quiet. They do not bake exactly the expected number: running out loses a sale and annoys a regular, while a leftover loaf costs only flour, so they bake a bit more than the average. A demand forecasting system does this for millions of products at once, with the "bit more" computed from a prediction interval and the cost of each mistake.

Requirements

Functional

• Forecast demand for every store x product x day, 1 to 28 days ahead, refreshed daily; and orders per store x hour for staffing, 1 to 14 days ahead.
• Produce quantiles (at least P10, P50, P90), not only a point forecast.
• Keep forecasts coherent: store forecasts add up to zone, city and total forecasts used by finance and logistics.
• Take known future inputs: holidays, local events, planned promotions and prices, store openings and closures.
• Support overrides by planners, logged and measured like any other forecast.

Non-functional

• Freshness: the nightly run finishes before ordering cut-off, for every series.
• Robustness: no missing forecasts; a fallback (seasonal naive) for any series where the model fails.
• Accuracy targets per level and horizon, tracked against a naive baseline.
• Explainability: planners need to see why a forecast jumped (a holiday, a promotion), or they will override it.

Out of scope

• The ordering optimizer and staff scheduler that consume the forecasts (they take quantiles as input).
• Pricing and promotion design.
• Real-time surge decisions at the minute level, such as rider repositioning.

Estimates

Assumed numbers for a large grocery delivery business, not any company's figures:

# Back-of-envelope sizing for a demand forecasting system. Every input is an assumption; change them and re-run.
cities, zones_per_city, stores_per_zone = 40, 25, 4
products_per_store = 8_000
horizon_days = 28                 # forecasts needed for ordering and staffing plans
quantiles = 3                     # P10, P50, P90
history_days = 2 * 365
zero_share = 0.6                  # share of store-product-days with no sales (intermittent demand)

stores = cities * zones_per_city * stores_per_zone
series = stores * products_per_store
print(f"stores                      {stores:,}")
print(f"store-product series        {series / 1e6:,.0f} million")
print(f"forecast values per run     {series * horizon_days * quantiles / 1e9:,.1f} billion (daily run)")
rows = series * history_days
print(f"training rows (2 years)     {rows / 1e9:,.1f} billion, {rows * zero_share / 1e9:,.1f} billion of them zeros")
levels = {"total": 1, "city": cities, "zone": cities * zones_per_city, "store": stores, "store-product": series}
print("series per hierarchy level  " + ", ".join(f"{k} {v:,}" for k, v in levels.items()))
print(f"one model per series        {series / 1e6:,.0f} million models to fit, monitor and debug")

stores                      4,000
store-product series        32 million
forecast values per run     2.7 billion (daily run)
training rows (2 years)     23.4 billion, 14.0 billion of them zeros
series per hierarchy level  total 1, city 40, zone 1,000, store 4,000, store-product 32,000,000
one model per series        32 million models to fit, monitor and debug

• Series: 4,000 stores x 8,000 products is 32 million store-product series. A nightly run with 28 days and 3 quantiles writes 2.7 billion numbers.
• Training data: two years of daily history is 23.4 billion rows, and with an assumed 60 percent of zero-sales days, 14 billion of them are zeros: demand at the bottom of the hierarchy is intermittentDemand with many periods of zero sales and occasional small positive values, typical of slow-moving products at a single store..
• Hierarchy: 1 total, 40 cities, 1,000 zones, 4,000 stores, 32 million store-products. Plans are made at several of these levels.
• One model per series is not an option: 32 million models to fit, monitor and debug is an operations problem before it is a modelling one.

Data and labels

The label is units sold, and the first trap is that sales are not demand. When a product is out of stock, sales are zero but demand was not; training on those days teaches the model that demand drops whenever stock runs out. Mark stock-out days (from inventory logs) and either drop them or treat them as censored observations, whose true demand is at least what was sold.

Other data issues:

• Calendar: public holidays per country and region, school holidays, paydays, daylight saving changes.
• Events: concerts, football matches, weather alerts; they come from external feeds and planners, with a start date known in advance.
• Promotions and prices: known in advance for the forecast horizon, so they are valid features. Past promotions explain past spikes; without them the model learns noise.
• Store changes: openings, closures, remodels, a new competitor next door. A new store has no history at all.

Features

Every feature must be known at forecast time. For a forecast made today for 28 days ahead, "sales yesterday" is not available for the target day; "sales 28 days before the target day" is. Getting this wrong is the most common way to build a forecasting model that looks brilliant in backtests and fails in production.

• Calendar: day of week, day of month, week of year, days to and from each holiday.
• Known future: promotion flag and depth, price, planned events, store open hours.
• Lagged demand: sales at lags of at least the horizon (28, 35, 364 days), rolling means and their ratios, all computed per series.
• Static: store, zone and city ids, store size, product category, shelf life, price tier.
• Hierarchy aggregates: zone and city demand trends, which carry signal to sparse store-product series.

Normalize targets by a per-series scale (recent mean demand), so one model can learn from a store selling 5 units a day and one selling 500.

Model choice and training

Global models

A global modelOne forecasting model trained on many related time series together, with series identifiers and features as inputs, instead of one separate model per series. is trained on all series at once. It learns patterns that are shared ("demand rises before a holiday", "promotions add about half") from millions of series, and applies them to series with little history. Typical choices:

• Gradient boosted trees (LightGBM or similar) on the feature table above: fast, strong, the usual workhorse.
• Deep models built for many series, such as DeepAR or the Temporal Fusion Transformer ( ML and deep learning for forecasting ), which output distributions directly and handle known future inputs natively.
• Classical local models (exponential smoothing, the Theta method, seasonal naive) as baselines and fallbacks.

Uber's forecasting introduction puts it plainly: you cannot know in advance which method will perform best, so compare several on a proper backtest.

Direct vs recursive

A directA forecasting strategy that trains a model to predict h steps ahead directly from information available now, rather than feeding its own one-step forecasts back in. model predicts day $t+h$ from information at $t$, possibly one model per horizon bucket. A recursive model predicts one day ahead and feeds its prediction back in as a lag. Direct avoids compounding errors; recursive needs fewer models. Many teams average both.

Prediction intervals

Train the same model with pinball lossThe loss for a quantile q: under-predictions cost q per unit and over-predictions (1 - q) per unit, so minimizing it yields the q-th quantile. for each quantile, or use a model that outputs a distribution. The loss for quantile $q$, actual $y$ and forecast $\hat{y}$ is

Worked example for $q=0.9$: forecasting 100 when the actual is 120 costs $0.9\times 20=18$; forecasting 140 costs $0.1\times 20=2$. Over-forecasting is nine times cheaper, so the model learns to sit high: at the 90th percentile.

High-level design

Deep dives

Deep dive 1: Hierarchical forecasts

Forecasts are needed at several levels: store-product for ordering, store for staffing, city for budgets, total for finance. Forecasting each level separately gives numbers that do not add up, and every team argues about whose is right. The forecasts must be coherentForecasts across a hierarchy are coherent when every aggregate equals the sum of the forecasts beneath it.. The textbook approaches, as described in Forecasting: Principles and Practice (Hyndman and Athanasopoulos):

• Bottom-up: forecast the bottom level and sum. Captures each series' own dynamics, but bottom series are noisy.
• Top-down: forecast the total and split it by proportions. Simple and reliable at the top, but the book notes it loses information, cannot capture series-specific dynamics or events, and with historical proportions tends to be less accurate at lower levels than bottom-up.
• Middle-out: forecast a middle level (zone), sum upwards, split downwards.
• Optimal reconciliation: forecast every level, then adjust all of them together. The MinT (minimum trace) method of Wickramasuriya and colleagues chooses the adjustment that minimizes the total forecast error variance of the coherent forecasts.

The sample below builds 60 stores in 3 cities and 12 zones, two years of daily orders with weekly patterns, city-specific growth, 10 holidays a year and random store promotions. Every forecast is made 28 days ahead.

# Demand forecasting on a synthetic hierarchy: 3 cities x 4 zones x 5 stores, 2 years of daily orders.
import numpy as np
from sklearn.ensemble import HistGradientBoostingRegressor as HGB

rng = np.random.default_rng(7)
S, T, H = 60, 2 * 365 + 56, 56                          # stores, days, test days; every forecast is made 28 days ahead
city, zone = np.arange(S) // 20, np.arange(S) // 5
level = rng.lognormal(4, 0.5, S)                         # store size: median about 55 orders a day before growth
dow = np.arange(T) % 7
week = 1 + 0.25 * np.isin(dow, [4, 5]) - 0.2 * (dow == 6)
hol = np.isin(np.arange(T) % 365, [0, 45, 120, 180, 250, 300, 330, 358, 359, 364])
promo = rng.random((S, T)) < 0.05
trend = 1 + 0.0004 * np.arange(T) * (1 + city[:, None])  # cities grow at different rates
mu = level[:, None] * week * trend * (1 + 0.8 * hol) * (1 + 0.5 * promo)
y = rng.poisson(mu).astype(float)

def feats(s, t, use_hol=True):                           # only information known 28+ days ahead
    lag = y[s, t - 28]; base = y[s, t - 56:t - 28].mean()
    f = [dow[t], promo[s, t], lag / base, city[s]] + ([hol[t]] if use_hol else [])
    return f, base
def table(days, use_hol=True):
    X, B, Y = [], [], []
    for t in days:
        for s in range(S):
            f, b = feats(s, t, use_hol); X.append(f); B.append(b); Y.append(y[s, t])
    return np.array(X), np.array(B), np.array(Y)
train, test = range(56, T - H), range(T - H, T)
Xtr, Btr, Ytr = table(train); Xte, Bte, Yte = table(test)
glob = HGB(max_iter=300, random_state=0).fit(Xtr, Ytr / Btr)          # one global model, scale-free target
pred = (glob.predict(Xte) * Bte).reshape(H, S)
naive = np.array([[y[s, t - 28] for s in range(S)] for t in test])  # seasonal naive: same weekday 4 weeks ago
act = Yte.reshape(H, S)
wape = lambda f, a: np.abs(f - a).sum() / a.sum()
agg = {"store": np.eye(S), "zone": np.eye(12)[zone].T, "city": np.eye(3)[city].T, "total": np.ones((1, S))}
print("WAPE by level (bottom-up: store forecasts summed, so every level adds up)")
for k, A in agg.items():
    print(f"  {k:6s} seasonal naive {wape(naive @ A.T, act @ A.T):.3f}   global model {wape(pred @ A.T, act @ A.T):.3f}")
share = y[:, T - H - 365:T - H].sum(1) / np.bincount(city, y[:, T - H - 365:T - H].sum(1))[city]
topdown = (pred @ agg["city"].T)[:, city] * share       # city total split by last year's store shares
print(f"top-down store WAPE {wape(topdown, act):.3f} vs bottom-up {wape(pred, act):.3f} (city totals identical)")
m = np.array([hol[t] for t in test]).repeat(S)
nohol = HGB(max_iter=300, random_state=0).fit(table(train, False)[0], Ytr / Btr)
p2 = nohol.predict(table(test, False)[0]) * Bte
print(f"holiday days WAPE: with holiday feature {wape(pred.ravel()[m], Yte[m]):.3f}, without {wape(p2[m], Yte[m]):.3f}")
lo, hi = (HGB(loss="quantile", quantile=q, max_iter=300, random_state=0).fit(Xtr, Ytr / Btr).predict(Xte) * Bte for q in (0.1, 0.9))
print(f"P10-P90 interval covers {np.mean((Yte >= lo) & (Yte <= hi)):.1%} of store-days (target 80%), "
      f"mean width {np.mean(hi - lo):.1f} orders vs mean demand {Yte.mean():.1f}")

WAPE by level (bottom-up: store forecasts summed, so every level adds up)
  store  seasonal naive 0.211   global model 0.098
  zone   seasonal naive 0.148   global model 0.066
  city   seasonal naive 0.121   global model 0.058
  total  seasonal naive 0.108   global model 0.057
top-down store WAPE 0.123 vs bottom-up 0.098 (city totals identical)
holiday days WAPE: with holiday feature 0.109, without 0.399
P10-P90 interval covers 80.1% of store-days (target 80%), mean width 32.5 orders vs mean demand 97.8

• Error shrinks as you aggregate. The global model's WAPEWeighted absolute percentage error: the sum of absolute errors divided by the sum of actuals, so large series count more and zero-demand days do not divide by zero. is 0.098 per store, 0.066 per zone, 0.058 per city and 0.057 in total: errors partly cancel when added up.
• The global model beats seasonal naive at every level, 0.098 against 0.211 per store, because it learns holidays, promotions and growth from all stores at once.
• Top-down hurts the bottom. Splitting the same city totals by last year's store shares gives a store WAPE of 0.123 instead of 0.098, with identical city totals: the shares cannot know which store has a promotion next month.

Deep dive 2: Global models, holidays and events

A holiday happens once a year per series. A local model sees two examples; a global model sees two per store, thousands in total, and learns that holidays raise demand by a similar factor everywhere. In the sample, removing the holiday feature raises WAPE on holiday days from 0.109 to 0.399: the model without it forecasts an ordinary day and misses the spike entirely.

Events need the same treatment: a feature that marks "match at the stadium within 2 km" lets the model learn the effect from every past match at every nearby store. Uber's forecasting introduction makes the same point about its own data: daily trips show trend and weekly seasonality, and December has more peak dates because of the many holidays in that month; it lists exogenous variables such as weather and concerts among the factors in choosing a method.

For a brand-new store, a global model with static features (size, city, zone demand) gives a sensible forecast from day one; a local model gives nothing.

Deep dive 3: Prediction intervals drive the plan

The ordering system does not want the median. For a perishable product where a lost sale costs 4 times as much as a wasted unit, the cost-minimizing order is the quantile

the classic newsvendorThe single-period inventory problem: order the quantity at which the probability of selling out equals the cost of a lost sale divided by the sum of that cost and the cost of a leftover unit. result. Staffing follows the same logic with the cost of an idle hour against the cost of slow deliveries. So the forecasting system's job is to deliver calibrated quantiles; each consumer picks its own.

Calibration matters as much as the centre. In the sample, the P10 to P90 interval from quantile models covers 80.1 percent of store-days, against a target of 80, with an average width of 32.5 orders for an average demand of 97.8. An interval that is too narrow looks confident and causes stock-outs; one that is too wide wastes stock. Uber's forecasting post shows two forecasts with the same point values and very different intervals, and notes the wider one calls for much higher capacity reserves.

Try it

1. At P50 with a shortage cost of 4, the plan runs short on 50 of 84 days, 967 units short against 834 left over, for a total cost of 4,701. Move to P80, the quantile the formula recommends: 19 days short, 249 units short, 2,023 left over and the lowest cost, 3,019. Go further to P99: no shortages, but 5,137 units left over and a cost of 5,137.
2. At P80, set holidays in the model to no: total cost rises to 3,648, and 213 of the 447 units short fall on the 4 holidays. Four days out of 84 cause about half the shortage.
3. With holidays back on, set the interval width to too narrow: at P80 the plan is now short on 30 days instead of 19 and the cost rises to 3,415. The plan is only as good as the interval's calibration.

Going deeper

Evaluation

Backtest the whole pipeline with rolling origins, as in Backtesting and forecast evaluation : pretend it is a date in the past, build features with only what was known then, forecast, and compare. Uber's post describes sliding and expanding windows, and suggests starting with an expanding window and switching to a sliding one once there is enough data. Then:

• Point accuracy with WAPE or MASE per level and horizon, always against a seasonal naive baseline (Uber's post recommends comparing against the naive forecast).
• Quantile accuracy with pinball loss and coverage of each interval, per level.
• Bias: the average signed error. A small persistent bias at the store-product level becomes a large one when summed.
• Decision metrics: stock-outs, waste and staff utilization in a simulation that replays the ordering logic on the backtest. A forecast that wins on WAPE but loses on waste is not better.
• Planner overrides: measure them like a model. If overrides make forecasts worse on average, show planners the evidence.

Monitoring

Track accuracy and coverage per level, horizon and category weekly; alert on bias that persists for several days; check that inputs arrived (a missing promotion feed silently removes all promotion uplifts); and watch the share of series falling back to the naive forecast.

Bottlenecks and failure modes

Trade-offs we made

Extensions

• Hourly forecasts for staffing, reconciled with the daily forecast (temporal hierarchies).
• Causal promotion effects: estimate uplift with experiments so the model does not confuse promotions with seasonality.
• Foundation models for time series as zero-shot baselines for new series ( ML and deep learning for forecasting ).
• Spatio-temporal demand for ride-hailing or delivery zones, where neighbouring zones share demand.

In the wild

Where you will see this: retailers ordering stock per store, grocery and food delivery planning pickers and riders, ride-hailing marketplaces positioning drivers, energy companies forecasting load per region, and cloud providers planning capacity all run a global forecasting model with quantiles and a reconciliation step.

Check yourself

Recap