25 August 2026

Forecasting the products that barely sell

A catalogue of two thousand SKUs usually has about fifty that sell every day and a very long tail that sells nothing for three weeks and then four at once. Forecasting the fifty is easy and nobody needs an app for it. The tail is where the money is stuck, and it is where ordinary forecasting quietly falls apart.

The failure is not that the numbers come out wrong. It is that they come out plausible. Run a moving average over a SKU that sells three units every eleven days and you get 0.27 units a day — a number describing no week that ever happened, and a reorder point built on it is wrong in both directions at once.

This is what levelin runs instead, and the two places the library did something we did not expect.

First, decide what kind of thing you are looking at

Before choosing a method you have to know whether demand is intermittent, and by how much. The standard answer is Syntetos–Boylan–Croston: two numbers per SKU.

ADI, average demand interval — how many periods pass between sales. CV², the squared coefficient of variation of the order sizes.

import numpy as np

def sbc(series: np.ndarray) -> tuple[float, float]:
    """Average demand interval and squared CV of non-zero demand sizes."""
    nonzero = series[series > 0]
    if nonzero.size < 2:
        return float("nan"), float("nan")

    adi = series.size / nonzero.size
    cv2 = float((nonzero.std(ddof=0) / nonzero.mean()) ** 2)
    return adi, cv2

# 52 weeks, sells in bursts of 2-5 roughly every third week
weekly = np.array([0, 0, 3, 0, 0, 2, 0, 0, 0, 4, 0, 0, 2, 0, 0, 5, 0, 0, 3, 0,
                   0, 0, 2, 0, 0, 4, 0, 0, 3, 0, 0, 2, 0, 0, 0, 5, 0, 0, 3, 0,
                   0, 2, 0, 0, 4, 0, 0, 3, 0, 0, 2, 0])

adi, cv2 = sbc(weekly)
print(f"ADI {adi:.2f}  CV2 {cv2:.3f}")

The CV² is over the non-zero sizes only. This is the detail most implementations get wrong, and it is worth being pedantic about: including the zeros measures how often the SKU sells, which is what ADI already tells you. Doing it twice makes every intermittent series look erratic, and erratic SKUs get a much larger safety buffer. You would be paying for the same fact twice, in inventory.

The cutoffs are ADI 1.32 and CV² 0.49, which gives four quadrants:

	CV² ≤ 0.49	CV² > 0.49
ADI < 1.32	SMOOTH	ERRATIC
ADI ≥ 1.32	INTERMITTENT	LUMPY

Those constants look arbitrary because they nearly are. They come from Syntetos, Boylan and Croston (2005), and Kostenko and Hyndman later showed the exact boundary is closer to 4/3 and depends on the smoothing parameter. We kept 1.32 because it is what the literature and every other implementation uses, and a SKU sitting exactly on that line is one whose classification does not matter much either way.

What the classification is for is choosing a method — and, for one quadrant, choosing not to.

Croston, SBA, TSB

Croston’s method is the standard answer for intermittent demand. Rather than smoothing the series, it smooths two series separately: the sizes of the non-zero demands, and the intervals between them. The forecast is one divided by the other.

That decomposition is the whole idea, and it is why a moving average cannot compete: the average is dragged toward zero by the gaps, while Croston keeps the gaps as information about timing rather than about quantity.

Croston’s estimator is biased upward, though. SBA (Syntetos–Boylan Approximation) corrects it by multiplying the result by 1 − α/2.

Here is the first thing that surprised us.

import inspect
from statsforecast import models

print(inspect.getsource(models._croston_sba))

def _croston_sba(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
) -> Dict[str, np.ndarray]:
    out = _croston_classic(y=y, h=h, fitted=fitted)
    out["mean"] *= 0.95
    if fitted:
        out["fitted"] *= 0.95
    return out

The correction is a hardcoded 0.95. Not 1 - alpha/2 — the literal number.

That is not a bug, and it took reading _croston_classic to see why:

import inspect
from statsforecast import models

source = inspect.getsource(models._croston_classic)
print("\n".join(line for line in source.splitlines() if "_ses_forecast" in line))

    ydp, ydf = _ses_forecast(yd, 0.1)
    yip, yif = _ses_forecast(yi, 0.1)

α is hardcoded to 0.1 as well, and 1 − 0.1/2 is exactly 0.95. The two constants agree because neither can move.

The practical consequence is the part worth knowing. CrostonClassic and CrostonSBA are both fixed at α = 0.1. CrostonOptimized fits α per series — and there is no SBA variant of it. So you get a tuned smoothing parameter with no bias correction, or the bias correction with a fixed one. You cannot have both, and no signature in the API hints at the trade.

We route to CrostonSBA by default and let backtesting overrule it per SKU, which sidesteps the argument: if the tuned α is genuinely better for a series, CrostonOptimized wins on that series’ own history and gets used.

TSB (Teunter–Syntetos–Babai) changes the question. Instead of smoothing the interval between demands, it smooths the probability of demand, and it updates that probability every period, including the empty ones.

That difference matters for exactly one thing, and it is the thing that made us keep it.

import numpy as np
import pandas as pd
from statsforecast import StatsForecast
from statsforecast.models import CrostonSBA, TSB

# Sold steadily, then stopped a year ago. A discontinued product.
history = [4, 0, 3, 5, 0, 4, 6, 0, 3, 4, 5, 0, 4, 3, 0, 5, 4, 0, 3, 4]
dead = [0] * 32

df = pd.DataFrame({
    "unique_id": "GFT-BOX-LRG",
    "ds": pd.date_range("2025-01-06", periods=len(history + dead), freq="W-MON"),
    "y": history + dead,
})

sf = StatsForecast(models=[CrostonSBA(), TSB(alpha_d=0.2, alpha_p=0.2)], freq="W-MON")
print(sf.forecast(df=df, h=1).round(4).to_string(index=False))

  unique_id         ds  CrostonSBA    TSB
GFT-BOX-LRG 2026-01-05      2.8669 0.0023

SBA still forecasts 2.87 units a week for a product that has not sold in over seven months — about 1,200 times what TSB says. It cannot do otherwise: it only ever looks at the non-zero demands and the intervals between them, and from that viewpoint nothing has changed — the last sale was still four units, the typical gap was still short. The dead period is invisible to it because there is nothing in it to smooth.

TSB sees thirty-two consecutive periods of no demand and decays the probability toward zero. It is looking at the same data and reaching the opposite conclusion, because it is asking a different question.

We do not use TSB as the default, because that same responsiveness makes it jumpy on live seasonal products. We use the ratio of the two as an obsolescence signal: when TSB has decayed far below SBA, the series is telling you the product is dying, whatever its recent sales suggest. That number feeds the dead-stock report rather than the reorder point.

Safety stock, and the term everyone forgets

Given a forecast, the reorder point is the demand you expect over the lead time plus a buffer for the times it runs high:

SS  = z * sqrt( L * sigma_d^2  +  d_mean^2 * sigma_L^2 )
ROP = d_mean * L  +  SS

Two things about this in practice.

L is the lead time plus the review period. If your supplier takes 14 days and you place orders fortnightly, stock has to cover 28 days, not 14 — you are not going to reorder tomorrow just because you dipped below the line today. Leaving the review period out is the most common way a reorder point comes out too low, and it fails in the least visible way: everything looks fine until the one time demand runs high near the end of a cycle.

σ_L is not optional. A supplier who is usually 14 days and occasionally 28 needs a much larger buffer than one who is always 14, and if you set σ_L to zero you under-buffer precisely the unreliable suppliers that need it most.

No merchant knows their supplier’s lead-time standard deviation. They do know “usually two weeks, worst I have seen is four”, so that is what we ask for, and treat the worst case as roughly a 95th percentile:

import math

def safety_stock(mean_daily, sigma_daily, lead_days, review_days, sigma_lead, z=1.65):
    horizon = lead_days + review_days          # NOT just the lead time
    variance = horizon * sigma_daily ** 2 + (mean_daily ** 2) * sigma_lead ** 2
    return z * math.sqrt(variance)

def sigma_lead_from(typical_days, worst_case_days):
    """Treat the worst case as ~2 standard deviations out."""
    return max(0.0, (worst_case_days - typical_days) / 2)

reliable = safety_stock(2.0, 1.4, 14, 14, sigma_lead_from(14, 14))
flaky    = safety_stock(2.0, 1.4, 14, 14, sigma_lead_from(14, 28))

print(f"reliable supplier: {reliable:5.1f} units")
print(f"flaky supplier:    {flaky:5.1f} units")

Same demand, same lead time, same service level. The unreliable supplier costs roughly twice the buffer, and that is the honest answer rather than a modelling artefact.

The estimate is deliberately crude. It is a stated assumption a merchant can argue with, which beats a precise-looking number derived from a figure they guessed.

Do not score this with MAPE

MAPE divides by the actual value. Intermittent demand is mostly zeros. That is the whole objection, and it is fatal rather than inconvenient:

import numpy as np

actual = np.array([0, 0, 3, 0, 0, 2, 0, 4, 0, 0])
predicted = np.array([0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8])

with np.errstate(divide="ignore", invalid="ignore"):
    mape = np.abs((actual - predicted) / actual)

print("per-period MAPE:", mape)
print("mean MAPE:", np.nanmean(mape[np.isfinite(mape)]))

Seven of ten periods are undefined. Dropping them, as most implementations silently do, scores the model only on the periods where something sold — exactly the periods an intermittent forecaster finds easiest, and it throws away every period where it predicted demand that never came.

We use MASE, which scales the error by what a naive forecast would have achieved on the same series, and wMAPE (nd in utilsforecast), which divides total absolute error by total demand rather than dividing period by period.

import numpy as np
import pandas as pd
from utilsforecast.losses import mase, nd

df = pd.DataFrame({
    "unique_id": ["SKU-1"] * 10,
    "ds": pd.date_range("2026-01-05", periods=10, freq="W-MON"),
    "y": [0, 0, 3, 0, 0, 2, 0, 4, 0, 0],
    "model": [0.8] * 10,
})
train = pd.DataFrame({
    "unique_id": ["SKU-1"] * 10,
    "ds": pd.date_range("2025-10-27", periods=10, freq="W-MON"),
    "y": [0, 2, 0, 0, 3, 0, 0, 2, 0, 4],
})

print(mase(df, models=["model"], seasonality=1, train_df=train).to_string(index=False))
print(nd(df, models=["model"]).to_string(index=False))

MASE below 1 means you beat the naive baseline. Above 1 means you did not, and for genuinely lumpy SKUs you often will not — which is information a merchant should be given rather than shielded from.

What the numbers actually came out as, and what we got wrong

We ran the whole routing table over the UCI Online Retail II dataset — a real transactional catalogue, not a simulation — and published the result: 55.7% beat the naive baseline, median MASE 0.878.

Then we went back and interrogated our own number, and it does not mean what we said it meant. Three things are wrong with it.

The denominator is not the catalogue. 55.7% is the share of the 3,314 SKUs we could score. The catalogue is 4,984. A third of it never gets scored at all — too few sales to backtest — and falls through to a default with nothing measured behind it. As a share of the catalogue, the figure is 37.0%.

Dead stock was carrying it. 396 of those scored SKUs had zero actual demand across the entire evaluation window. They score a median MASE of 0.028 and 99% of them “beat naive” — they are not being forecast, they are being scored against an empty test set. Strip them out and the number is 49.8%, which is below the bar we set ourselves.

The champion graded its own homework. The winning method per SKU was chosen by lowest error across the same three windows it was then scored on. Use the routing default with no per-SKU selection and the same run reads 43.8%.

Here is the part that survives. Select the champion on earlier windows and score it on one the model has never seen, restricted to the SKUs a merchant actually orders for: median MASE 0.877, against a four-week moving average’s 0.938. Per SKU, that wins 42%, ties 18% and loses 40%.

So we do not beat a moving average on most of your products. What we do is lose less badly when we lose: about 3–5% lower average error, and materially fewer catastrophic misses. That is a smaller claim than the one we started with, and it is the one we can defend.

One more, because it is the least flattering and the most useful. Split the catalogue by what demand is doing rather than by its SBC class, and the engine is near-perfect on stock that is dying (median MASE 0.27) and close to useless on stock that is taking off (1.46, with only 35% beating naive). A stockout on a ramping product is the expensive failure, and it is precisely the one this family of methods cannot see: Croston anchors its level on a mostly-zero history, and every model in the table emits a flat path. No model we tested, in the routing table or out of it, fixes that. Detecting it and saying so is the only honest response, and that is what we are building next.

The rule we are most confident about

For LUMPY SKUs — rare and wildly variable — levelin produces no suggested quantity at all.

The case that settled it: a SKU averaging about three units a week whose single largest order had been forty-five. The arithmetic wants a thirty-unit safety stock. The shelf already held ten weeks of cover. Every step of that calculation is correct and the answer is a purchase order for cash you cannot get back.

So those SKUs go to a separate list, with their history, for a human to decide. A merchant knows about the wholesale enquiry that produced the forty-five. The statistics only know it happened once.

Being honest about what a model cannot do turns out to be a feature you can ship.

levelin is demand forecasting for Shopify stores with long tails. The forecasting engine is Python and statsforecast; the app is Remix and Postgres.