Quantile Regression for Forecasting a Range, Not a Point

A point forecast answers a question almost nobody has. Inventory needs the level demand will not exceed 95% of the time; capacity planning needs the same at a different percentile. Quantile regression produces those directly, and the loss function that makes it work is a three-line asymmetry you can verify by hand.

The pinball loss

For a target quantile τ between 0 and 1, an observation y and a prediction q:

L_τ(y, q)  =  τ · (y − q)        if  y ≥ q     (under-predicted)
              (1 − τ) · (q − y)  if  y <  q     (over-predicted)

It is absolute error with a thumb on the scale. At τ = 0.5 the two branches have equal weight and the loss is half the absolute error, whose minimiser is the median. At τ = 0.9, being one unit too low costs 0.9 and being one unit too high costs 0.1, so the estimate is pushed up until only 10% of observations remain above it. The asymmetric weight is the quantile specification; there is nothing else to it.

Roger Koenker and Gilbert Bassett established the estimator in their 1978 Econometrica paper “Regression Quantiles”, and the loss goes by several names — pinball loss, check loss, quantile loss — all of which are this expression. It is convex, so a linear quantile regression is a linear program with a global optimum, and it is differentiable everywhere except at y = q, which is why gradient-based fitting works without any special handling.

Why its minimiser is the quantile

Differentiate the total loss over a sample with respect to q. Each observation above q contributes −τ, because raising q towards it reduces that term. Each observation below contributes +(1 − τ). So for a sample of n:

d/dq  Σ L_τ  =  (1 − τ) · n_below  −  τ · n_above

set to zero:      (1 − τ) · n_below  =  τ · n_above
with n_below + n_above = n:

                  n_below / n  =  τ

The optimum sits where a fraction τ of the sample lies below it, which is the definition of the τ-quantile. No distributional assumption entered anywhere — this holds for any data, skewed, multimodal or full of zeros. That last property is why the loss is the natural choice for intermittent demand, where percentage-based losses are undefined.

The loss evaluated at five candidates

Take ten observed values and τ = 0.9. This is the whole dataset:

y = [ 2, 4, 5, 7, 8, 10, 12, 15, 20, 30 ]      τ = 0.9

Evaluate the total loss at five candidate predictions. Each figure below is 0.9 × Σ(y − q) over values above q, plus 0.1 × Σ(q − y) over values below:

q = 12    0.9 × (3+8+18)   = 26.1     0.1 × 36  = 3.6      total = 29.7
q = 15    0.9 × (5+15)     = 18.0     0.1 × 57  = 5.7      total = 23.7
q = 20    0.9 × (10)       =  9.0     0.1 × 97  = 9.7      total = 18.7
q = 25    0.9 × (5)        =  4.5     0.1 × 142 = 14.2     total = 18.7
q = 30    0.9 × (0)        =  0.0     0.1 × 187 = 18.7     total = 18.7
q = 35    0.9 × (0)        =  0.0     0.1 × 237 = 23.7     total = 23.7

The loss falls to 18.7 and then goes flat across the entire interval from 20 to 30 before rising again. That flat region is not a rounding artefact — it is what the derivation predicts. With n = 10 and τ = 0.9, the condition n_below / n = τ is satisfied by any q with exactly nine values below it, and every point in (20, 30] qualifies. The sample quantile is not unique when n·τ is a whole number, and different software resolves the tie differently.

Two things follow that matter in practice. Do not read a small difference between two candidate quantile forecasts as meaningful when the loss surface around them is flat — the optimiser is choosing arbitrarily within a region. And expect extreme quantiles to be unstable: a τ = 0.99 estimate on 200 observations is determined by the top two values, so one outlier moves it entirely. If you need a 99th percentile, either fit a distribution and read the tail off it, or get more data; the empirical loss cannot manufacture information the sample does not contain.

Quantile crossing

The usual way to produce an interval is to fit one model per quantile — separate fits for 0.1, 0.5 and 0.9. Nothing in those three optimisations knows about the others, so nothing prevents the fitted 0.1 curve from rising above the fitted 0.9 curve in some region of the input space. When it does, you have a prediction interval whose lower bound exceeds its upper bound, which is not an interval.

Crossing is not a bug in the implementation; it is what independent fitting permits. It shows up most in extrapolation, at the edges of the feature space where each fit has little data and the estimates are least stable. Three responses:

• Rearrange. Sort the predicted quantiles at each point. This is monotone rearrangement, and it is known not to make the estimates worse — a sorted set of quantiles is at least as close to the truth as the unsorted one. It is a post-hoc fix and it is usually enough.
• Fit jointly. One model with multiple quantile outputs, trained on the summed pinball losses. The Temporal Fusion Transformer does exactly this, emitting several quantiles per horizon step from one network. Shared representation reduces crossing without forbidding it.
• Constrain. Parameterise the higher quantiles as the lower one plus a non-negative increment, so crossing is structurally impossible. This is the only option that guarantees the result, and it costs you the ability to fit each quantile independently.

Quantiles do not add across horizons

This is the error that does the most damage, because the resulting number looks reasonable. A model that gives you the 95th percentile of daily demand for each of the next seven days does not let you sum those seven to get the 95th percentile of the weekly total.

The reason is that the sum of per-day 95th percentiles is the demand level that would be reached only if every one of the seven days simultaneously hit its own 95th percentile. Unless the days are perfectly correlated, the probability of that is far below 5%, so the sum is a much more extreme quantile of the weekly distribution than you asked for — you over-stock, and the size of the error grows with the number of periods you summed. The same logic applies down the product hierarchy: quantiles do not aggregate across SKUs either.

The correct approaches are to simulate whole sample paths and take the quantile of the simulated totals, or to fit the model directly at the aggregate you care about — a weekly model for a weekly decision, per choosing a forecast granularity. Both are more work than adding seven numbers, which is why adding seven numbers is so common.

Finally, validate the intervals rather than trusting them. Over a rolling-origin evaluation, count how often the actual fell below your 90th percentile forecast. If the answer is not close to 90%, the model is not calibrated, and pinball loss alone will not tell you — a forecast can have a good average loss and systematically wrong coverage. Report both the loss and the empirical coverage per quantile, and watch coverage degrade with the horizon, because it will.