Continuous Probabilistic Forecasting

A continuous probabilistic forecast predicts a continuous distribution for a future outcome rather than a value or a few cut points. Instead of “next week’s demand is 500,” the model reports \(Y_{t+1} \sim \mathcal{N}(500, 30^2)\), from which any quantile, interval or event probability can be computed.

The ladder is short. A point forecast says 500 units. A quantile forecast says the 10th is 450, the median 500, the 90th 560. A continuous probabilistic forecast names the distribution, and the quantiles follow from it.

Three worked cases show the range. Tomorrow’s temperature as \(N(25, 2^2)\) gives roughly 68% between 23 and 27 °C and 95% between 21 and 29. Weekly demand as LogNormal(6.2, 0.3) on the log scale gives a median near 493 units, a 90th percentile of about 724, and \(P(\text{demand} > 700) \approx 0.12\). Weekly returns as Student-\(t\) with 5 degrees of freedom let a risk desk compute value-at-risk and expected shortfall. That third case is where the whole approach earns its keep, and where it is most easily got wrong.

The family is a decision, and the tail is where it shows

Naming a distribution means naming a family. Fix the mean and the standard deviation and the family is still free, so it is fair to ask how much of the answer the remaining freedom controls.

import numpy as np
from scipy.stats import norm, t, laplace

m, sd = 0.01, 0.02
"""One week of returns, three families, all with the same mean and the same
   standard deviation. The families differ only in shape."""

nu = 5
fams = (("Normal", norm(loc=m, scale=sd)),
        ("Student-t, 5 df", t(nu, loc=m, scale=sd / np.sqrt(nu / (nu - 2)))),
        ("Laplace", laplace(loc=m, scale=sd / np.sqrt(2))))

N = 4_000_000
u = (np.arange(N) + 0.5) / N
print(f"  mean {m}, standard deviation {sd} in every row:")
print(f"{'family':>18}{'sd (check)':>12}{'5th pct':>11}{'1st pct':>10}"
      f"{'99% VaR':>10}{'P(< -5%)':>11}{'ES at 1%':>11}")
for name, D in fams:
    q = D.ppf(u)
    print(f"{name:>18}{q.std():>12.4f}{D.ppf(0.05):>11.4f}{D.ppf(0.01):>10.4f}"
          f"{-D.ppf(0.01):>10.4f}{D.cdf(-0.05):>11.6f}{q[q < D.ppf(0.01)].mean():>11.4f}")
print("  the 5th percentile is within half a point across all three. The 1st")
print("  percentile already separates them, and the probability of a 5% loss")
print(f"  differs by a factor of {laplace(loc=m, scale=sd / np.sqrt(2)).cdf(-0.05) / norm(loc=m, scale=sd).cdf(-0.05):.0f}"
      " between the Normal and the Laplace")

print(f"\n  the same comparison as a central interval, where they agree:")
print(f"{'family':>18}{'80% interval':>26}{'98% interval':>28}")
for name, D in fams:
    lo80, hi80 = D.ppf(0.10), D.ppf(0.90)
    lo98, hi98 = D.ppf(0.01), D.ppf(0.99)
    print(f"{name:>18}{'[' + format(lo80, '.4f') + ', ' + format(hi80, '.4f') + ']':>26}"
          f"{'[' + format(lo98, '.4f') + ', ' + format(hi98, '.4f') + ']':>28}")
print("  central intervals differ less than tail probabilities here. The choice of")
print("  family can matter centrally too; central coverage alone does not verify")
print("  the tail probabilities needed for the decision")

rng = np.random.default_rng(0)
truth = t(nu, loc=m, scale=sd / np.sqrt(nu / (nu - 2)))
sample = truth.ppf(rng.random(2000))
fit_n = norm(loc=sample.mean(), scale=sample.std(ddof=1))
print(f"\n  and 2,000 observations do not settle it. Fitting a Normal to data")
print(f"  that came from the t:")
print(f"    sample mean {sample.mean():>9.6f}   sample sd {sample.std(ddof=1):>9.6f}")
print(f"    fitted Normal says P(< -5%) = {fit_n.cdf(-0.05):.6f}")
print(f"    the truth is                 {truth.cdf(-0.05):.6f}")
print(f"    observed in the sample:      {np.mean(sample < -0.05):.6f}"
      f"   ({int(np.sum(sample < -0.05))} of {len(sample):,} days)")
print("  the events that decide the answer are the ones there are fewest of.")
print("  A distributional forecast is only as good as the family it assumed,")
print("  and the data argues about that family least where it matters most")
#   mean 0.01, standard deviation 0.02 in every row:
#             family  sd (check)    5th pct   1st pct   99% VaR   P(< -5%)   ES at 1%
#             Normal      0.0200    -0.0229   -0.0365    0.0365   0.001350    -0.0433
#    Student-t, 5 df      0.0200    -0.0212   -0.0421    0.0421   0.005862    -0.0590
#            Laplace      0.0200    -0.0226   -0.0453    0.0453   0.007185    -0.0595
#   the 5th percentile is within half a point across all three. The 1st
#   percentile already separates them, and the probability of a 5% loss
#   differs by a factor of 5 between the Normal and the Laplace
#
#   the same comparison as a central interval, where they agree:
#             family              80% interval                98% interval
#             Normal         [-0.0156, 0.0356]           [-0.0365, 0.0565]
#    Student-t, 5 df         [-0.0129, 0.0329]           [-0.0421, 0.0621]
#            Laplace         [-0.0128, 0.0328]           [-0.0453, 0.0653]
#   central intervals differ less than tail probabilities here. The choice of
#   family can matter centrally too; central coverage alone does not verify
#   the tail probabilities needed for the decision
#
#   and 2,000 observations do not settle it. Fitting a Normal to data
#   that came from the t:
#     sample mean  0.009722   sample sd  0.020329
#     fitted Normal says P(< -5%) = 0.001653
#     the truth is                 0.005862
#     observed in the sample:      0.007000   (14 of 2,000 days)
#   the events that decide the answer are the ones there are fewest of.
#   A distributional forecast is only as good as the family it assumed,
#   and the data argues about that family least where it matters most

What that means for the methods

The usual toolkit divides by how the family gets chosen. Parametric models — ARIMA with Gaussian errors, GARCH for volatility — assume one outright, so their tails are an assumption rather than a finding. Bayesian models produce posteriors naturally and inherit the same question one level up, in the likelihood. Simulation and ensemble methods generate sample paths and read the distribution off them, which trades the assumption for a sample size: the empirical distribution can only describe a tail as far as it has draws. Among neural models, DeepAR is an autoregressive network with likelihood outputs, temporal fusion transformers can produce predictive distributions, and mixture density networks output the parameters of a mixture — all of which pick a family too, and mixtures mainly widen the choice rather than remove it.

Evaluation uses strictly proper scoring rules, so that reporting your true predictive distribution is the optimum: CRPS, which generalizes the Brier score to continuous outcomes, and the log score, which rewards putting density where the outcome lands. Alongside them, a calibration check on interval coverage.

One caution follows from the measurement above: an average score over all outcomes is dominated by the middle, where the families agree. A model can win on CRPS and still be the wrong one for a value-at-risk limit. If the decision lives in the tail, score the tail — coverage of the 1% quantile, or a weighted CRPS — rather than trusting an average that the tail barely moves.

Why it is worth the trouble

One forecast answers every question downstream. Inventory buffers, capacity limits and financial risk limits are all thresholds on the same distribution, and a point forecast cannot supply any of them without a separate uncertainty estimate bolted on afterward. Producing the distribution once and reading thresholds off it keeps those answers consistent with each other, which a collection of separately fitted numbers does not guarantee.

References

• SciPy documentation, scipy.stats.norm, scipy.stats.t and scipy.stats.laplace. docs.scipy.org