Forecasting many series at once

In practice you rarely forecast one thing. You forecast every product in every store, including the ones launched last week with almost no history. This section covers the choice that follows: a separate model per series, which is simple but learns nothing across them, or one global model that shares what it learns between them. It also covers hierarchies that have to add up, and demand that is zero on most days.

One model per series, or one for all#

With ten thousand products, fitting a separate classical model to each is simple but learns nothing across products. A single global model, such as gradient boosting with product features and lags, learns shared patterns and helps products with short histories. Global models have done very well in forecasting competitions for this reason.

                      Local: one model per series    Global: one model, all series
short-history series  fits badly or not at all       borrows shape from similar ones
per-series quirks     captured exactly               captured only if a feature says so
training time         10,000 small fits              one fit
deployment            10,000 artefacts to version    one artefact
a new product         no model until it has history  predicted from day one
diagnosis             read one model                 read feature importances, once
scales to             hundreds of series             millions

The global model needs the series identity supplied as features rather than as an identifier, because an id is meaningless to a tree. Give it what actually distinguishes the series: category, price band, store size, how long it has existed, its average volume, its typical weekday shape. Then a product launched last week is predicted from products that resemble it, which is the whole point and is exactly what a local model cannot do.

Hierarchies must add up#

If store forecasts do not sum to the region forecast, planners get two different answers. Reconciliation methods adjust forecasts at all levels so they add up, often improving accuracy at the same time. Libraries in the statsforecast family implement them.

The two simple approaches are worth knowing by hand, because they encode different beliefs about where the signal is. Bottom-up forecasts each leaf and adds, which is right when the leaves are stable and noisy in independent ways. Top-down forecasts the total, which is smoother and easier to predict, then splits it by historical shares. Neither is generally better, and the honest answer is to check which wins at every level of your own hierarchy.

import numpy as np

def bottom_up(leaf_forecasts):
    return float(np.sum(leaf_forecasts)), np.asarray(leaf_forecasts, dtype=float)

def top_down(total_forecast, historical_shares):
    shares = np.asarray(historical_shares, dtype=float)
    shares = shares / shares.sum()
    return float(total_forecast), total_forecast * shares

leaves = np.array([100.0, 60.0, 40.0])     # per-store forecasts
total, per_store = bottom_up(leaves)
print(total, per_store)                    # 200.0 [100.  60.  40.]

total, per_store = top_down(210.0, leaves) # the total model says 210
print(total, per_store.round(1))           # 210.0 [105.  63.  42.]

Mostly zeros#

Many products sell nothing on most days. Standard models forecast a small fraction every day, which is never what happens. Methods for intermittent demand, such as Croston's, forecast the size and timing of sales separately. Evaluate them on the decision they support, such as stock-outs, not only on error.

Croston's method is two exponential smoothers: one over the sizes of the non-zero sales, one over the gaps between them. The forecast rate is the smoothed size divided by the smoothed gap. What it buys you is not a better average — the average is easy — but the two numbers a stock decision actually needs: how much when it sells, and how long you wait.

import numpy as np

def croston(series, alpha=0.1):
    """Returns the demand rate per period, from smoothed size and smoothed gap."""
    size = interval = None
    gap = 0
    for x in series:
        gap += 1
        if x > 0:
            size = x if size is None else alpha * x + (1 - alpha) * size
            interval = gap if interval is None else alpha * gap + (1 - alpha) * interval
            gap = 0
    return size / interval, size, interval

d = np.array([0, 0, 0, 5, 0, 0, 0, 0, 7, 0, 0, 6, 0, 0, 0, 0, 0, 4.0])
rate, size, interval = croston(d)
print(round(float(d.mean()), 3))     # 1.222  a flat model forecasts this every day
print(round(rate, 3))                # 1.229  the same rate, arrived at differently
print(round(size, 2), round(interval, 2))   # 5.15 4.19

Which is also why MASE and RMSE are poor judges here. A model that forecasts 1.22 every day has a low absolute error and orders the wrong amount at the wrong time. Score these on the decision: how often the item was out of stock, and how much capital sat in a warehouse. If those two numbers are not available, forecasting intermittent demand accurately is not something anyone can confirm you have done.