Probabilistic Forecasting & Prediction Intervals

TL;DR — A point forecast is one number ("load at 6 p.m. will be 8,400 MW"). A probabilistic forecast is a range with a probability ("90% chance load is between 7,900 and 8,900 MW"). Real decisions — how much reserve generation to hold, how many crews to staff, how much power to buy — depend on the range, not the mean. This article shows where a raw range comes from (a parametric model, a quantile-regression model trained with pinball loss, or the spread of an ensemble of models), how to score it with PICP (coverage) and MPIW (width), and why those raw ranges are often miscalibrated — which is exactly what Conformal Prediction: Coverage-Guaranteed Intervals fixes.

---

1. Simple explanation

Imagine you run the control room for a regional power grid. Tomorrow at 6 p.m. people come home, turn on ovens and air conditioners, and demand ("load") spikes. Your job: make sure supply always meets demand. If you under-supply, the lights go out. If you over-supply, you burned money and fuel for nothing.

A point forecast says: "6 p.m. load = 8,400 MW." Useful, but it hides the thing you actually need — how sure are we? If the truth could easily be 9,000 MW, you must hold extra reserve. If it's nailed down to ±50 MW, you can run lean.

A probabilistic forecast answers the real question. Instead of one number, it gives a distribution: a most-likely value and how wide the uncertainty is. From that you can read a prediction interval — "90% of the time, load lands between 7,900 and 8,900 MW."

Analogy — the weather app. A point forecast is the app saying "high of 78°F." A probabilistic forecast is "70% chance of rain, high between 74 and 82." You dress and plan differently for "70% rain" than for "5% rain," even if the headline temperature is identical. The grid operator is the same: you buy reserves against the chance of a high-load surprise, not against the single expected number. The interval is the part you actually act on.

The catch: it is easy to produce an interval, and hard to produce an honest one. A model can claim "90% interval" while the true value falls outside it 30% of the time. That gap between claimed and achieved coverage is the central problem of this whole sub-module.

---

2. Diagram

   POINT vs PROBABILISTIC FORECAST  (hourly electricity load, MW)

   load
   (MW)
   9000 |                         .--''''''--.  q95  (upper 90% band)
        |                    _.-''            ''-._
   8500 |   actual ->  o   .'      q50 (median forecast)  '.
        |            o   o'___________________________________'.___ o
   8000 |        o  '  _-''                                   ''-_ ' o
        |      o   .-''            q05  (lower 90% band)         ''-.
   7500 |   o   .-'                                                  '-.
        +----+----+----+----+----+----+----+----+----+----+----+----+---
         12   13   14   15   16   17   18   19   20   21   22   23  hour

   POINT forecast  = the q50 line only .......... "8400 MW"
   INTERVAL (90%)  = the band q05 .. q95 ......... "7900 .. 8900 MW"
   COVERAGE (PICP) = fraction of 'o' actuals that land INSIDE the band
   WIDTH   (MPIW)  = average vertical thickness of the band (q95 - q05)

   THREE PLACES A RAW BAND COMES FROM
   ----------------------------------
   (a) PARAMETRIC : assume a shape (Gaussian), band = y_hat +/- z*sigma
   (b) QUANTILE   : train models to output q05 / q50 / q95 directly
   (c) ENSEMBLE   : run N models, band from the SPREAD (std) of their forecasts
                    all three give a RAW band -> often MIScalibrated -> calibrate

---

3. How it works

3.1 Point vs probabilistic — what each one is

A point forecast is a summary of a distribution (usually its mean or median). A probabilistic forecast keeps the distribution so you can read any quantile you need.

3.2 Quantiles and prediction intervals

A quantile q_τ is the value below which a fraction τ of outcomes fall. The median is q_0.50. To build a central 90% prediction interval, take two quantiles that leave 5% in each tail:

   90% interval = [ q_0.05 , q_0.95 ]      (miss 5% low + 5% high = 10% outside)

More generally a (1−α) interval uses [ q_(α/2) , q_(1−α/2) ]. For 90% coverage, α = 0.10, so q_0.05 and q_0.95. Quantiles are attractive because they make no shape assumption — the distribution can be skewed (load usually is; big upside spikes are more common than big downside ones) and quantiles still describe it faithfully.

3.3 Where the raw uncertainty comes from

There are several standard sources. All produce a raw band; none is automatically calibrated.

The ensemble idea is worth a sentence. Train several different models (different seeds, different features, different algorithms — this family includes deep ensembles, MC dropout, and bootstrap resampling). When they agree, you are probably in familiar territory and uncertainty is low; when they disagree, something unusual is happening and uncertainty is high. The spread of their forecasts is a cheap, general uncertainty signal. It is one standard source among several — not a special or novel construction.

3.4 Quantile regression and pinball loss

To predict a single quantile τ, you do not minimize squared error — you minimize pinball loss (a.k.a. quantile loss). It is an asymmetric absolute error: under-predicting and over-predicting cost different amounts, and the ratio of those costs is exactly what pins the prediction to quantile τ.

   L_τ(y, ŷ) = max( τ·(y − ŷ) , (τ − 1)·(y − ŷ) )

     if y > ŷ (under-predict): cost = τ·(y − ŷ)        weight τ
     if y < ŷ (over-predict) : cost = (1 − τ)·(ŷ − y)  weight (1 − τ)

For τ = 0.9, under-predicting is weighted 0.9 and over-predicting 0.1 — a 9:1 penalty ratio. Minimizing that pushes ŷ up until only 10% of actuals sit above it: that point is the 90th percentile. Fit one model at τ=0.05, one at 0.50, one at 0.95, and you have a median plus a 90% interval.

3.5 Why raw intervals are often miscalibrated

A model can claim a 90% interval yet be wrong about the width. Common causes:

• Wrong shape assumption (parametric Gaussian on skewed load).
• Optimization slack — the pinball minimizer is approximate; tree/NN models don't hit the exact quantile on unseen data.
• Distribution shift — heatwaves, holidays, and new appliances move the target away from training conditions.
• Ensemble spread is a heuristic — the std of a few models has no reason to equal the true standard deviation; it's typically too narrow (overconfident).

The symptom is measurable: your claimed-90% band covers, say, only 78% of test actuals. That is a calibration failure, and it is exactly what Conformal Prediction: Coverage-Guaranteed Intervals repairs with a distribution-free coverage guarantee. First, though, you have to measure coverage — the next tools do that.

3.6 Scoring a probabilistic forecast

PICP and MPIW trade off: you can always hit 100% coverage by making the band infinitely wide, and you can always make MPIW tiny by making the band useless. The goal is the narrowest band that still hits the target coverage. Report them together — never one alone.

---

4. The math

4.1 Pinball loss, worked

Target quantile τ = 0.9, actual y = 100.

   Case A — prediction ŷ = 90 (UNDER-predict, ŷ < y):
       y − ŷ = 100 − 90 = 10
       τ·(y − ŷ)     = 0.9 · 10  =  9.0
       (τ−1)·(y − ŷ) = -0.1 · 10 = -1.0
       L = max(9.0, -1.0) = 9.0        <- under-shooting the 90th pct hurts a LOT

   Case B — prediction ŷ = 110 (OVER-predict, ŷ > y):
       y − ŷ = 100 − 110 = -10
       τ·(y − ŷ)     = 0.9 · (-10)  = -9.0
       (τ−1)·(y − ŷ) = -0.1 · (-10) =  1.0
       L = max(-9.0, 1.0) = 1.0        <- over-shooting hurts 9x LESS

Same size error (10 MW off) costs 9.0 when low but 1.0 when high. The optimizer therefore keeps raising ŷ until being-too-low becomes rare — until only 10% of actuals exceed it. That equilibrium point is the 90th percentile. Symmetric MSE could never do this; the asymmetry is the whole trick.

4.2 PICP and MPIW, worked

Ten test hours. Model outputs a 90% band [q05, q95] per hour; we check each actual.

   hour  q05   q95   actual y   inside?   width (q95-q05)
    1    7900  8900   8400       yes         1000
    2    8100  9100   8600       yes         1000
    3    7600  8600   8550       yes         1000
    4    8200  9300   9500       NO          1100   <- actual above q95
    5    7400  8400   7800       yes         1000
    6    7700  8700   8300       yes         1000
    7    8000  9000   8200       yes         1000
    8    7500  8500   8100       yes         1000
    9    7800  8800   8700       yes          1000
   10    8300  9400   8900       yes         1100

   PICP = (# inside) / N = 9 / 10 = 0.90     -> matches target 1-α = 0.90  (well calibrated)
   MPIW = mean width = (8*1000 + 2*1100)/10 = 10200/10 = 1020 MW

Here PICP = 0.90 exactly on target — good. Suppose a different, overconfident model had only 7 of 10 inside: PICP = 0.70, far below the promised 0.90. Its bands are too narrow — a calibration failure you would only catch by measuring. To compare two models fairly, look at both: at equal PICP, the one with smaller MPIW wins.

4.3 The parametric band (for contrast)

If you assume Gaussian residuals with standard deviation σ, the (1−α) band is ŷ ± z_(1−α/2)·σ. For 90%, z_0.95 ≈ 1.645. With ŷ = 8400 and σ = 300:

   band = 8400 ± 1.645·300 = 8400 ± 493.5 = [7906.5 , 8893.5]  MW

Clean and cheap — but only correct if load residuals really are Gaussian. When they are skewed (fat upper tail on hot evenings), this band is systematically wrong, which is why quantile models and calibration exist.

---

5. Real code

Quantile regression with gradient-boosted trees (LightGBM if available, else scikit-learn's GradientBoostingRegressor, which supports quantile loss natively), plus a small ensemble-spread band, all scored with PICP / MPIW / pinball.

"""Probabilistic hourly load forecasting: quantile regression + ensemble spread,
scored with PICP / MPIW / pinball loss.

Two raw interval sources are built and compared:
  (a) quantile regression at tau = 0.05 / 0.50 / 0.95  -> a direct 90% band
  (b) an ensemble of models -> a Gaussian-shaped band from the forecast SPREAD
Both are RAW (uncalibrated). Calibration is handled in the conformal article.
"""
import numpy as np

rng = np.random.default_rng(0)

# ---- 1. A small, synthetic-but-realistic hourly load problem -----------------
# Features: hour-of-day, temperature, is_weekend. Target: load in MW, skewed up.
N = 4000
hour = rng.integers(0, 24, N)
temp = rng.normal(20, 8, N)                         # deg C
weekend = rng.integers(0, 2, N)
# daily double-peak shape + temperature (AC) effect + weekend dip
base = (7000
        + 900 * np.sin((hour - 8) / 24 * 2 * np.pi)   # morning/evening shape
        + 25 * (temp - 20) ** 2 / 3                    # hot -> AC load (skewed up)
        - 400 * weekend)
noise = rng.normal(0, 200, N) + rng.exponential(120, N)  # asymmetric (fat upper tail)
y = base + noise
X = np.column_stack([hour, temp, weekend])

# chronological-style split (train / test); no shuffling of the tail
ntr = 3000
Xtr, ytr, Xte, yte = X[:ntr], y[:ntr], X[ntr:], y[ntr:]

# ---- 2. Quantile regression at tau = 0.05, 0.50, 0.95 ------------------------
try:
    from lightgbm import LGBMRegressor
    def make_q(tau):
        return LGBMRegressor(objective="quantile", alpha=tau,
                             n_estimators=300, learning_rate=0.05,
                             num_leaves=31, verbosity=-1)
except Exception:                                   # fallback: sklearn
    from sklearn.ensemble import GradientBoostingRegressor
    def make_q(tau):
        return GradientBoostingRegressor(loss="quantile", alpha=tau,
                                         n_estimators=300, learning_rate=0.05,
                                         max_depth=3)

q_lo = make_q(0.05).fit(Xtr, ytr)
q_md = make_q(0.50).fit(Xtr, ytr)
q_hi = make_q(0.95).fit(Xtr, ytr)

lo = q_lo.predict(Xte)
md = q_md.predict(Xte)
hi = q_hi.predict(Xte)
lo, hi = np.minimum(lo, hi), np.maximum(lo, hi)     # guard against quantile CROSSING

# ---- 3. Ensemble-spread band (one standard source among several) -------------
# Train several point models on bootstrap resamples; band = mean +/- 1.645*std.
from sklearn.ensemble import ExtraTreesRegressor
members = []
for s in range(8):
    idx = rng.integers(0, ntr, ntr)                 # bootstrap sample
    m = ExtraTreesRegressor(n_estimators=100, max_depth=8, random_state=s)
    m.fit(Xtr[idx], ytr[idx])
    members.append(m)
preds = np.stack([m.predict(Xte) for m in members])  # shape (8, n_test)
ens_mean = preds.mean(0)
ens_std = preds.std(0)                               # DISAGREEMENT = uncertainty
z = 1.645                                            # 90% under a Gaussian assumption
ens_lo, ens_hi = ens_mean - z * ens_std, ens_mean + z * ens_std

# ---- 4. Scoring: PICP, MPIW, pinball -----------------------------------------
def picp(y, lo, hi):
    """Coverage: fraction of actuals inside [lo, hi]. Target ~ 1 - alpha."""
    return np.mean((y >= lo) & (y <= hi))

def mpiw(lo, hi):
    """Mean prediction-interval width. Smaller is better AT the target coverage."""
    return np.mean(hi - lo)

def pinball(y, yhat, tau):
    """Average pinball (quantile) loss for a single quantile tau."""
    d = y - yhat
    return np.mean(np.maximum(tau * d, (tau - 1) * d))

print("Quantile regression  PICP=%.3f  MPIW=%.0f MW" % (picp(yte, lo, hi), mpiw(lo, hi)))
print("Ensemble spread      PICP=%.3f  MPIW=%.0f MW" % (picp(yte, ens_lo, ens_hi),
                                                        mpiw(ens_lo, ens_hi)))
print("pinball(q05)=%.1f  pinball(q50)=%.1f  pinball(q95)=%.1f"
      % (pinball(yte, lo, 0.05), pinball(yte, md, 0.50), pinball(yte, hi, 0.95)))

# Typical output: BOTH raw bands MISS the 0.90 target. The quantile band lands a bit
# under (~0.87). The bootstrap-ensemble band undercovers BADLY (often ~0.2), because
# its spread captures only model DISAGREEMENT (epistemic), not the irreducible noise
# (aleatoric) in load -- a vivid reminder that raw spread is a heuristic, not a sigma.
# Closing both gaps is exactly what conformal calibration does -> conformal-prediction.md

The load-bearing lines: np.minimum/np.maximum fixes quantile crossing (a real failure where q05 > q95 on some rows); ens_std is the disagreement signal; and the printed PICP will usually land below 0.90 — the raw miscalibration this article warns about.

---

6. Real-world example

Scenario: day-ahead reserve procurement for a regional operator.

You forecast tomorrow's 6 p.m. load. Two teams hand you outputs:

• Team Point: "6 p.m. load = 8,400 MW." You have no idea how much reserve to buy. If you guess and load hits 9,000, you shed load (a blackout event) — catastrophic and headline-making.
• Team Probabilistic: "median 8,400 MW; 90% interval [7,900, 8,900]." Now you can procure reserves to cover the upper end you're willing to insure against.

Put numbers on it. Reserves cost about $40 / MW / hour to hold. The cost of unserved energy (a shortfall) is roughly $9,000 / MWh (a standard order-of-magnitude "value of lost load").

   Strategy A — cover to the point forecast (8400 MW):
       If true load = 8900 (inside the 90% band, entirely plausible):
       shortfall = 500 MW for 1 hour  ->  0.5 GWh? no: 500 MW*1h = 500 MWh
       expected shortfall cost ~ 500 MWh * $9000 = $4,500,000 exposure event

   Strategy B — cover to q95 = 8900 MW:
       extra reserve held = 8900 - 8400 = 500 MW
       reserve cost = 500 MW * $40 = $20,000 for the hour
       shortfall risk in the covered range ~ 0 (you insured the upper tail)

For $20,000 of reserve you removed a $4.5M tail exposure. That trade is only visible because you had an interval — the point forecast literally cannot express "cover to the 95th percentile." And it only works if the interval is honest: if your "90%" band actually covers 78% of evenings, you are under-reserved and the tail bites anyway. That is why the next step after producing a band is measuring PICP and, if it misses, calibrating it.

---

7. Interview questions companies actually ask

Q1 [easy] (Amazon energy/retail forecasting, utilities) "Point vs probabilistic forecast —
   what's the difference and why prefer probabilistic?"
  A A point forecast is one number (the mean/median); a probabilistic forecast is a full
    distribution / quantiles / an interval. Decisions under uncertainty — reserves, staffing,
    inventory, bidding — need the RANGE and the tail, not just the center. The interval is the
    part you actually act on; the point forecast hides your risk.

Q2 [easy] (Google, Meta forecasting) "What is a 90% prediction interval, in quantiles?"
  A [q0.05, q0.95]: the band from the 5th to the 95th percentile, leaving 5% of outcomes in
    each tail so 10% total fall outside. Generally a (1-alpha) interval is
    [q_(alpha/2), q_(1-alpha/2)].

Q3 [medium] (grid operators, forecasting teams) "Why train with pinball loss instead of MSE
   for a quantile?"
  A MSE is symmetric and its minimizer is the MEAN. Pinball loss is asymmetric: for quantile
    tau it weights under-prediction by tau and over-prediction by (1-tau). Minimizing it pushes
    the prediction to the point where exactly tau of the mass lies below — i.e. the tau-th
    quantile. The asymmetry is what targets a specific quantile.

Q4 [medium] (Amazon, utilities) "Your quantile model outputs q05 > q95 for some rows. What
   happened and how do you fix it?"
  A Quantile CROSSING: the three quantiles are fit as independent models, so nothing enforces
    monotonicity. Cheap fix is to sort/clip per row (min for lower, max for upper). Better fixes:
    joint/monotone quantile models, or isotonic post-processing. Always mention it — it signals
    hands-on experience.

Q5 [medium] (Meta, Google) "What do PICP and MPIW measure, and why report both?"
  A PICP = coverage = fraction of actuals inside the interval; it should ~ 1-alpha. MPIW = mean
    interval width; smaller is better. They trade off: infinite width gives 100% coverage but
    zero information. The goal is the NARROWEST band that still HITS the target coverage, so you
    must report them together.

Q6 [medium] (forecasting platforms) "Where does the raw uncertainty even come from?"
  A Three standard sources: (a) PARAMETRIC — assume a shape (Gaussian) and use yhat +/- z*sigma;
    (b) QUANTILE models trained with pinball loss that output q05/q50/q95 directly; (c) ENSEMBLE
    / model-disagreement — run several models and use the SPREAD (std) of their forecasts. All
    three are RAW and usually need calibration.

Q7 [hard] (Amazon, Google, utilities) "You built a '90%' interval but it covers only 78% of
   test points. What's wrong and what do you do?"
  A It's MIScalibrated / overconfident — too narrow. Causes: wrong shape assumption, optimizer
    slack, distribution shift, or (for ensembles) the spread being a heuristic not a true sigma.
    Fix: recalibrate. The distribution-free, guaranteed way is conformal prediction — hold out a
    calibration set, measure the residual quantile, and widen the band to hit 90%.

Q8 [hard] (grid ops, energy trading) "Load has a fat UPPER tail on hot evenings. Why does that
   break a Gaussian parametric band, and what beats it?"
  A A Gaussian is symmetric, so yhat +/- z*sigma under-covers the heavy upper tail and over-covers
    the light lower tail — wrong on BOTH sides. Quantile regression makes no shape assumption and
    can put q95 far above q50 while keeping q05 close, capturing the skew. Conformal methods (esp.
    CQR) then guarantee the coverage.

Q9 [hard] (Meta, Amazon) "How does an ENSEMBLE give uncertainty, and why can't you trust its
   spread directly?"
  A Train several models (seeds, features, algorithms; deep ensembles / MC dropout / bootstrap).
    Where they AGREE, uncertainty is low; where they DISAGREE, it's high — so the std of their
    forecasts is a cheap uncertainty signal. But that std has no reason to equal the true sigma;
    it's typically too NARROW (overconfident). Treat it as a raw SCORE and calibrate it (e.g.
    conformally) before promising coverage.

Q10 [medium] (utilities, staffing/ops) "Give a decision that needs the interval, not the mean."
  A Reserve procurement: you buy reserves to cover the UPPER quantile of load you're willing to
    insure against (e.g. q95), not the mean. The mean can't express 'cover the 95th percentile.'
    Same logic for on-call staffing, safety stock, and energy-purchase hedging.

---

8. When to use / tradeoffs

   USE probabilistic forecasting when:
     ✓ the decision is asymmetric — a shortfall costs far more (or less) than a surplus
     ✓ you must size reserves / staff / inventory / bids against a TAIL, not a mean
     ✓ downstream consumers can act on a range (risk limits, cover-to-quantile rules)

   PICK the source by the situation:
     • parametric (yhat +/- z*sigma) — fast, fine when residuals really are near-Gaussian
     • quantile regression           — when the distribution is SKEWED / heavy-tailed (load is)
     • ensemble spread               — when you already run several models and want a cheap signal

   HONEST LIMITS:
     ✗ RAW bands are usually miscalibrated — measure PICP before you trust any '90%' claim
     ✗ quantile CROSSING can produce q05 > q95 — sort/clip or use monotone models
     ✗ a parametric Gaussian band is wrong for skewed load — under/over-covers the tails
     ✗ ensemble spread is a HEURISTIC, not a true sigma — typically overconfident
     ✗ wider isn't 'safer' for free — MPIW is a cost; a useless-wide band helps no decision
     ✗ i.i.d. assumptions are shaky in time series (autocorrelation, drift) — see conformal caveat

   RULE OF THUMB: produce a band, MEASURE PICP/MPIW, and if PICP misses the target, CALIBRATE
   it (conformal) rather than hand-tuning the width.

---

9. Summary + related articles

• A point forecast is one number; a probabilistic forecast is a distribution / quantiles / an interval — and the interval is what decisions act on.
• A 90% interval is [q_0.05, q_0.95]; generally [q_(α/2), q_(1−α/2)] for (1−α) coverage.
• Quantile regression trains directly for a quantile by minimizing pinball loss, whose asymmetry (weight τ vs 1−τ) pins the prediction to that quantile.
• Raw uncertainty has several standard sources: parametric shape, quantile models, and ensemble / model-disagreement spread — the last is one option among several, useful and cheap but heuristic.
• Score bands with PICP (coverage ≈ 1−α) and MPIW (width, smaller-is-better) — always together — plus pinball loss for quantile accuracy.
• Raw bands are commonly miscalibrated (often overconfident); that gap is what calibration fixes.

Related: Conformal Prediction: Coverage-Guaranteed Intervals · Forecast Evaluation & Backtesting · Time-Series Forecasting & Uncertainty — Interview Questions

Resources

• Koenker & Bassett — Quantile Regression (foundational) — https://www.jstor.org/stable/1913643
• Gneiting & Raftery — Strictly Proper Scoring Rules (CRPS, pinball) — https://doi.org/10.1198/016214506000001437
• LightGBM quantile objective — https://lightgbm.readthedocs.io/
• scikit-learn quantile GradientBoosting — https://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_quantile.html
• Lakshminarayanan et al. — Deep Ensembles for predictive uncertainty — https://arxiv.org/abs/1612.01474
• GEFCom2014 probabilistic load forecasting competition — https://doi.org/10.1016/j.ijforecast.2016.02.001