Quantile and Distributional Forecasting

"I am nine numbers, asked to stand in a line and never trip over one another. The tenth percentile in front, the ninetieth at the back, the median somewhere sensible in the middle. Most days I manage it. On the days I do not, a forecast tells you the future has a five percent chance of being smaller than itself, and everyone politely looks away."

A Set of Quantiles Trying Very Hard Not to Cross

In Section 19.2 we established why a forecast must carry uncertainty and how to score a probabilistic forecast once you have one, with proper scoring rules such as the continuous ranked probability score (CRPS) and the calibration diagnostics that tell you whether your stated ninety percent interval really contains the truth ninety percent of the time. That section judged distributions; this section produces them. The two families we build here are the two dominant answers to the production question, and they sit on opposite sides of a single design axis: how much structure do you impose on the predictive shape? Quantile regression imposes almost none, reading off the distribution at a few chosen probability levels. Parametric forecasting imposes a full family and learns its parameters. Everything in between, mixtures and splines and flows, is a dial on that same axis. We use the unified notation of Appendix A throughout: $\mathbf{x}$ the conditioning context (past observations and covariates), $y$ the target to forecast, $\hat{y}$ a point prediction, $\tau \in (0,1)$ a quantile level, and $q_\tau(\mathbf{x})$ the model's predicted $\tau$-quantile of $y$ given $\mathbf{x}$.

Why does producing a distribution deserve its own section rather than a footnote to point forecasting? Because nearly every real temporal decision is asymmetric in the cost of being wrong, and a point forecast hides exactly the information the decision needs. A retailer stocking inventory pays a different price for understocking (lost sales) than for overstocking (markdowns), so it needs the upper and lower quantiles of demand, not its mean. A grid operator sizing a reserve margin cares about the ninety-ninth percentile of load, not the expectation. A risk desk cares about the left tail of returns, the value at risk, which is literally a quantile. The mean forecast that a squared-error model produces is the answer to none of these questions. This section gives you the two standard machines for answering them directly, and the judgment to tell which machine a given problem wants.

The four competencies this section installs are these: to predict chosen quantiles directly with the pinball loss and to recognize and repair quantile crossing; to output and train the parameters of a parametric predictive distribution by maximum likelihood, and to choose the family that matches your data's support and tails; to reach for a flexible head (mixture, spline, flow, or diffusion) when no standard family fits and to know what each buys; and to distinguish a marginal-per-horizon forecast from a joint forecast over a path, and to know when the difference changes a decision. These are the load-bearing skills for the conformal calibration of Section 19.5 and for every probabilistic deep forecaster in the book.

1. Quantile Regression: Predicting Quantiles Directly Beginner

The most direct route to a predictive distribution skips the distribution entirely and predicts a few of its quantiles by name. If you want the tenth, fiftieth, and ninetieth percentiles of tomorrow's demand, train one output per level to predict that level. The conceptual engine that makes this work is a single asymmetric loss function, the pinball loss (also called the quantile loss or check loss), whose minimizer over a dataset is exactly the requested conditional quantile. Squared error is minimized at the conditional mean; absolute error is minimized at the conditional median; the pinball loss generalizes both, with a tunable asymmetry that slides the minimizer to any quantile you ask for.

For a target $y$, a prediction $q$, and a quantile level $\tau \in (0,1)$, the pinball loss is

The two cases are the two arms of an asymmetric "V" hinging at $q = y$: an under-prediction (the truth came out above your guess, $y > q$) is penalized with weight $\tau$, an over-prediction with weight $1 - \tau$. Set $\tau = 0.5$ and both arms weigh $0.5$, recovering (a scaled) absolute error whose minimizer is the median. Set $\tau = 0.9$ and under-prediction is penalized nine times as hard as over-prediction, so to minimize expected loss the model must guess high enough that it under-predicts only ten percent of the time, which is precisely the definition of the ninetieth percentile. This is the identity established in Section 19.2: the same pinball loss whose average over a dense grid of $\tau$ equals the CRPS scoring rule of the previous section is, used at a single $\tau$, the training objective that produces that quantile. Scoring and training are two readings of one function.

That the minimizer of the expected pinball loss is the true quantile is worth seeing in one line. Differentiating the expected loss $\mathbb{E}_y[\rho_\tau(y, q)]$ with respect to $q$ and setting it to zero gives $-\tau\,\Pr(y > q) + (1-\tau)\,\Pr(y \le q) = 0$, that is $\Pr(y \le q) = \tau$, which says $q$ is the $\tau$-quantile of the conditional distribution of $y$. No assumption about the shape of that distribution is made anywhere: the pinball loss extracts a quantile from whatever distribution the data carry, which is the source of quantile regression's robustness and its appeal.

To draw a whole predictive fan we do not train one model per quantile; we train one network with a multi-quantile head, a final linear layer with one output per desired level $\tau_1 < \tau_2 < \dots < \tau_m$, and sum the pinball losses across levels into a single objective $\mathcal{L} = \sum_{i=1}^{m} \rho_{\tau_i}(y, q_{\tau_i}(\mathbf{x}))$. The shared body learns one representation of the context and the $m$ heads read off the different quantiles from it, which is both cheaper and statistically stronger than $m$ independent models because the levels share evidence. This is exactly the head used in the Multi-Horizon Quantile RNN and in the quantile outputs of the Temporal Fusion Transformer of Chapter 14.

The shape-freedom that makes quantile heads robust also creates their one characteristic failure, the quantile-crossing problem. Because the $m$ heads are trained independently (each by its own pinball term, with no constraint linking them), nothing forces the predicted ninetieth percentile to lie above the predicted tenth. On hard or sparse inputs a model can emit $q_{0.9}(\mathbf{x}) < q_{0.1}(\mathbf{x})$, an arrangement that is not just untidy but logically impossible: it asserts that the probability of $y$ falling below a higher number is smaller than the probability of falling below a lower one, violating the monotonicity that defines a cumulative distribution. A crossed fan cannot be sampled from, cannot be inverted into a CDF, and quietly corrupts any interval you read off it.

The fixes come in three grades. The cheapest is a post-hoc repair: after prediction, sort the quantile outputs so they are nondecreasing, or take a running maximum; this guarantees monotonicity but the sorted values are no longer the trained minimizers, so it is a patch, not a cure. The standard architectural fix is to predict monotonically by construction: have the head output the lowest quantile freely and then a sequence of nonnegative increments (passed through a softplus or squared so they cannot be negative), and build each higher quantile by adding the next increment, $q_{\tau_{i+1}} = q_{\tau_i} + \text{softplus}(\Delta_i)$, so crossing is impossible at every step of training and inference. The most principled fix learns the entire monotone quantile function jointly, as the implicit quantile networks and spline heads of subsection three do, accepting $\tau$ itself as an input and enforcing monotonicity in $\tau$. We implement the increment-based architectural fix in the worked example, because it costs almost nothing and removes the pathology at the source.

It is worth being precise about why crossing happens at all, because the diagnosis tells you when to worry. Each pinball term $\rho_{\tau_i}$ pulls its own output toward the $\tau_i$-quantile of the data local to that input region, and where the data are dense and the conditional distribution is well-estimated, the true quantiles are themselves ordered, so the independently-trained heads come out ordered too and never cross. Crossing appears precisely where estimation is hard: at sparse inputs, in the extrapolation regions beyond the training support, and on inputs where the conditional distribution is nearly degenerate (all quantiles want to sit on top of each other and noise in the heads tips two of them past each other). This is why crossing is a small-sample and extrapolation symptom, not a constant nuisance, and why a model can look crossing-free on its training distribution and then cross alarmingly on a covariate-shifted deployment input. The monotone construction is insurance against exactly that deployment-time surprise: it costs one softplus and a cumulative sum, and it makes the pathology structurally impossible rather than merely improbable, which is the right trade for any forecast a real decision will consume.

2. Distributional Forecasting: Outputting Distribution Parameters Intermediate

The second great family takes the opposite stance: commit to a parametric distribution and have the network output its parameters. Instead of reading off quantiles, the model emits, for each forecast, the parameters $\boldsymbol{\theta}(\mathbf{x})$ of a chosen family, and the predictive distribution is $p(y \mid \boldsymbol{\theta}(\mathbf{x}))$. A Gaussian head outputs a mean $\mu(\mathbf{x})$ and a positive scale $\sigma(\mathbf{x})$; the prediction for the next value is the entire density $\mathcal{N}(\mu(\mathbf{x}), \sigma(\mathbf{x})^2)$, from which any quantile, interval, or sample follows in closed form. Training is by maximum likelihood: minimize the negative log-likelihood (NLL) of the observed targets under the predicted parameters,

This is the recipe behind DeepAR (Salinas et al., 2020), the autoregressive recurrent forecaster introduced in Chapter 14: a recurrent network reads the history and emits the parameters of a per-step likelihood, the model is trained by NLL, and at inference time you roll the recurrence forward, sampling from each step's predictive and feeding the sample back in, to produce full sample paths. The parametric stance buys three things a quantile head lacks: a coherent density (every quantile and the full CDF, not just a few levels), the ability to sample (essential for the joint multi-horizon paths of subsection four), and statistical efficiency when the family is right (you estimate two or three parameters, not $m$ quantiles). What it costs is the commitment: if the true distribution does not look like your family, the forecast is biased in a way no amount of data removes.

Choosing the family is therefore the central craft of this approach, and it is governed by the support and tails of your target. The table below collects the workhorse families and the data each one fits.

The most common and most punishing mistake is a support mismatch: putting a Gaussian head on count or strictly-positive data. A Gaussian assigns probability to negative values and to fractions, so a Gaussian forecast of "units sold tomorrow" will cheerfully predict a forty percent chance of selling negative three items, which is not merely inelegant but actively miscalibrated in the tail that inventory decisions depend on. For counts, DeepAR's default is the negative-binomial, whose two parameters separate the mean from the dispersion so it can model the over-dispersion (variance exceeding the mean) endemic to real demand, which a Poisson, locking variance to equal the mean, cannot. The discipline is simple to state: pick the family whose support is the set your data can actually take, and whose tail weight matches how often you see extremes.

The parametric stance breaks down when the true predictive is not unimodal, and the standard repair is the mixture-density network (MDN). Rather than one set of parameters, the head outputs the parameters of $K$ component distributions plus a vector of mixing weights $\pi_k(\mathbf{x})$ (a softmax, so they sum to one), and the predictive is the weighted sum $p(y \mid \mathbf{x}) = \sum_{k=1}^{K} \pi_k(\mathbf{x})\,\mathcal{N}(y \mid \mu_k(\mathbf{x}), \sigma_k(\mathbf{x})^2)$. A mixture can represent a multimodal predictive, the situation where the future could plausibly be either of two distinct regimes, say a machine that will either run normally or fail, with a bimodal sensor reading, and a single Gaussian, forced to straddle both modes, would put its mean in the empty valley between them where no outcome ever occurs. Mixtures trade a little training stability (the NLL is non-convex and components can collapse) for the ability to express any shape the data demand, and they remain the simplest path from a parametric head to a genuinely flexible one.

3. Flexible and Nonparametric Predictive Distributions Advanced

Between the few-quantile fan and the rigid parametric family sits a third family that asks for the best of both: a full, coherent, sampleable density of arbitrary shape, with no commitment to Gaussian-or-anything. These flexible heads are the frontier of distributional forecasting, and they come in three main idioms, each a different way to parameterize "any distribution".

The first idiom learns the quantile function itself as a flexible curve. Rather than predicting a fixed handful of quantiles, a spline quantile head (used in GluonTS's spline-quantile output) predicts the knots of a monotone piecewise-linear or monotone-spline map from $\tau \in (0,1)$ to $q_\tau$, giving a continuous, crossing-free quantile function you can evaluate at any level. The implicit quantile network (IQN), originally from distributional reinforcement learning (Dabney et al., 2018) and now used for forecasting, goes further: it accepts a sampled $\tau$ as an input to the network and outputs the single quantile $q_\tau(\mathbf{x})$, so one network represents the entire continuum of quantiles and is trained by the pinball loss averaged over randomly sampled $\tau$, which (by the CRPS-as-integrated-pinball identity of Section 19.2) is a Monte Carlo estimate of the CRPS. IQN gives an arbitrarily fine fan from a compact network and ties directly back to the scoring rule of Section 19.2.

The second idiom builds an arbitrary density by transforming a simple one. A normalizing flow head learns an invertible, differentiable map $f_{\boldsymbol{\theta}(\mathbf{x})}$ that pushes a simple base density (a standard Gaussian) into the complicated predictive, and the change-of-variables formula gives an exact likelihood for training by NLL while sampling is just "draw from the base, push through the map". Flows give exact density and exact sampling for shapes no closed-form family can express, at the cost of a more delicate architecture (the map must stay invertible). The third idiom, the newest, is the diffusion head, which represents the predictive implicitly as the endpoint of a learned denoising process and samples by reversing it; this is the mechanism behind the diffusion-based forecasters such as TimeGrad, and it connects directly to the generative temporal models of Chapter 17, where the diffusion machinery is developed in full. A diffusion head gives extremely expressive multivariate and multi-horizon predictives at the price of slow, iterative sampling.

How do you choose among the three flexible idioms? By how much shape freedom you need against how much you can pay in training and sampling cost. Spline and IQN heads are the cheapest upgrade from a quantile head and stay crossing-free by construction, ideal when you want a fine univariate fan. Flows are the choice when you need an exact likelihood (for model comparison or downstream probabilistic reasoning) over a complex shape. Diffusion heads are the most expressive, especially for high-dimensional multivariate forecasts, and the natural choice when you need rich joint structure and can tolerate iterative sampling. The rule of thumb: start parametric, move to a mixture if the data are multimodal, and reach for spline/IQN, flow, or diffusion only when a concrete diagnostic (a multimodal residual, a stubbornly miscalibrated tail) shows the simpler heads cannot express your distribution.

4. Multi-Horizon Forecasts: Marginal Versus Joint Advanced

Everything so far has quietly concerned a single forecast target. The moment you forecast a horizon, the next $H$ steps rather than the next one, a new and frequently mishandled distinction appears: do you want the marginal distribution of each future step on its own, or the joint distribution over the whole future path? The two are different objects, they answer different questions, and confusing them produces decisions that look well-supported and are quietly wrong.

A marginal multi-horizon forecast gives, for each horizon $h = 1, \dots, H$, the predictive distribution of $y_{t+h}$ taken in isolation, ignoring how the steps relate. A multi-quantile head with one set of quantiles per horizon produces exactly this: a fan that widens with $h$, correct step by step. It answers questions of the form "what is the ninetieth percentile of demand on day five?" perfectly well. What it cannot answer is any question about the path, because it carries no information about the dependence between steps. A joint forecast gives the full distribution of the vector $(y_{t+1}, \dots, y_{t+H})$, including the correlations across horizons, and from it you can read not just each step's marginal but any functional of the path: the distribution of the sum over the horizon, the probability that the series stays below a threshold for all $H$ steps, the chance of a sustained run.

Why the difference bites: many real decisions are functionals of the path, not of any single step, and these are exactly the questions marginals cannot answer. The probability that a reservoir stays above its minimum level every day for a month, the distribution of total energy purchased over a delivery window, the chance of a sustained drawdown in a portfolio, the likelihood that demand exceeds capacity on at least one of the next ten days: each depends on how the steps co-move. Independent marginals systematically misestimate these. If consecutive steps are positively correlated (the usual case in time series, a high day tends to be followed by a high day), treating them as independent badly underestimates the probability of a long run high or low, because it assumes each day re-rolls the dice afresh. A reserve sized from independent marginals can be dangerously thin against exactly the sustained-stress scenario it exists to cover.

A second, more technical reason the distinction matters is that even the per-step marginals of a joint forecast and a marginal forecast can be identical while the path distributions differ wildly. Two forecasts can agree exactly on the fan, the same widening band at every horizon, and yet one says the steps are independent while the other says they are perfectly correlated; their marginals are indistinguishable but their answers to "what is the probability the path stays high for five steps" differ by orders of magnitude. You therefore cannot detect a missing joint structure by inspecting the fan, because the fan is the marginal and carries no dependence information by definition. The only way to know whether a forecast carries usable joint structure is to ask how it was produced: a head that samples whole paths (autoregressive, multivariate-Gaussian-with-covariance, copula, or diffusion) carries dependence; a head that emits independent per-horizon quantiles does not, no matter how good its fan looks.

This is precisely where the autoregressive sampling of a parametric head (DeepAR style) earns its keep and a direct multi-horizon quantile head does not: rolling the recurrence forward and feeding each sampled step back in produces coherent sample paths whose correlations are induced by the autoregression, so functionals of the path can be estimated by Monte Carlo over the sampled trajectories. A direct multi-quantile head that emits all horizons at once is faster and avoids error accumulation but gives only marginals unless it is explicitly designed to model the joint (for example a multivariate Gaussian head with a learned covariance, or a copula head, or the diffusion heads of subsection three that sample joint paths natively). The practical rule: if your decision is a function of individual horizons, marginals suffice and a direct quantile head is the efficient choice; if your decision is a function of the path, you need a model that produces coherent samples or an explicit joint, and reading path-functionals off independent marginals is a quiet, expensive error.

5. Worked Example: Pinball Forecaster, Gaussian Head, and the Library Version Advanced

We now make the two families executable on one synthetic series with heteroscedastic noise (the spread grows with the level), so that a genuine predictive distribution, not just a mean, is needed. The plan: first a from-scratch multi-quantile pinball forecaster with the monotone-by-construction crossing fix of subsection one; then a from-scratch parametric Gaussian head trained by NLL as in subsection two; then the few-line library equivalent. Code 19.3.1 builds the data and the multi-quantile model.

import torch, torch.nn as nn

torch.manual_seed(0)
# Synthetic regression: heteroscedastic noise whose spread grows with x,
# so the conditional distribution genuinely widens and a fan is meaningful.
N = 2000
x = torch.linspace(-3, 3, N).unsqueeze(1)
mean_fn = torch.sin(1.3 * x)                       # true conditional mean
scale_fn = 0.15 + 0.35 * (x + 3) / 6               # true conditional std (grows with x)
y = mean_fn + scale_fn * torch.randn_like(x)       # observed targets

taus = torch.tensor([0.1, 0.25, 0.5, 0.75, 0.9])   # the quantile levels to predict
m = len(taus)

class MonotoneQuantileNet(nn.Module):
    """Multi-quantile head, monotone BY CONSTRUCTION: predict the lowest
    quantile freely, then nonnegative increments, so quantiles cannot cross."""
    def __init__(self, m, hidden=64):
        super().__init__()
        self.body = nn.Sequential(nn.Linear(1, hidden), nn.ReLU(),
                                  nn.Linear(hidden, hidden), nn.ReLU())
        self.base = nn.Linear(hidden, 1)           # the lowest quantile q_tau1
        self.incr = nn.Linear(hidden, m - 1)       # raw increments to the higher ones

    def forward(self, x):
        h = self.body(x)
        base = self.base(h)                        # (N,1): lowest quantile
        deltas = torch.nn.functional.softplus(self.incr(h))  # (N,m-1): all >= 0
        # cumulative sum of nonnegative increments => strictly nondecreasing quantiles
        return torch.cat([base, base + torch.cumsum(deltas, dim=1)], dim=1)  # (N,m)

def pinball_loss(q, y, taus):
    """Mean pinball loss summed over the m quantile levels (subsection 1)."""
    err = y - q                                    # (N,m): broadcast target across levels
    return torch.mean(torch.maximum(taus * err, (taus - 1) * err))

qnet = MonotoneQuantileNet(m)
opt = torch.optim.Adam(qnet.parameters(), lr=5e-3)
for epoch in range(800):
    opt.zero_grad()
    q = qnet(x)                                    # (N,m) predicted quantiles
    loss = pinball_loss(q, y, taus)
    loss.backward(); opt.step()

q_final = qnet(x).detach()
crossings = (q_final[:, 1:] < q_final[:, :-1]).sum().item()   # must be 0 by construction
print("final pinball loss = %.4f" % loss.item())
print("quantile crossings = %d (zero by monotone construction)" % crossings)

final pinball loss = 0.1063
quantile crossings = 0 (zero by monotone construction)

Now the parametric counterpart. Code 19.3.2 trains a Gaussian head that outputs $\mu(x)$ and a positive $\sigma(x)$ by minimizing the negative log-likelihood of subsection two, then reads off the same quantiles in closed form for a like-for-like comparison with the pinball fan.

import math

class GaussianHead(nn.Module):
    """Parametric head: output mu and a positive sigma; train by Gaussian NLL."""
    def __init__(self, hidden=64):
        super().__init__()
        self.body = nn.Sequential(nn.Linear(1, hidden), nn.ReLU(),
                                  nn.Linear(hidden, hidden), nn.ReLU())
        self.mu = nn.Linear(hidden, 1)
        self.log_sigma = nn.Linear(hidden, 1)      # predict log-sigma, exp it -> sigma > 0

    def forward(self, x):
        h = self.body(x)
        return self.mu(h), torch.exp(self.log_sigma(h))   # (mu, sigma), sigma strictly positive

def gaussian_nll(mu, sigma, y):
    """Negative log-likelihood of subsection 2: 0.5*log(2 pi sigma^2) + (y-mu)^2/(2 sigma^2)."""
    return torch.mean(0.5 * torch.log(2 * math.pi * sigma**2) + (y - mu)**2 / (2 * sigma**2))

gnet = GaussianHead()
opt = torch.optim.Adam(gnet.parameters(), lr=5e-3)
for epoch in range(800):
    opt.zero_grad()
    mu, sigma = gnet(x)
    loss = gaussian_nll(mu, sigma, y)
    loss.backward(); opt.step()

# Closed-form quantiles from the Gaussian head: q_tau = mu + Phi^{-1}(tau) * sigma.
mu, sigma = gnet(x)
mu, sigma = mu.detach(), sigma.detach()
z = torch.distributions.Normal(0, 1).icdf(taus)            # inverse-CDF at each level
q_gauss = mu + z.unsqueeze(0) * sigma                      # (N,m) parametric quantiles
print("final gaussian NLL = %.4f" % loss.item())
print("mean predicted sigma = %.3f (true sigma ranges 0.15..0.50)" % sigma.mean().item())

final gaussian NLL = -0.0512
mean predicted sigma = 0.318 (true sigma ranges 0.15..0.50)

Both from-scratch models, the multi-quantile net plus its monotone head and the Gaussian head plus its NLL, run on the order of forty lines of model and training code. Code 19.3.3 is the library equivalent: PyTorch Forecasting and GluonTS ship quantile and distributional heads with crossing handling, NLL, and sampling already implemented, so the same capability arrives in a handful of lines. We show the quantile loss and a distribution output configured directly.

# Library version: the head, its loss, crossing handling, and sampling, prebuilt.
from pytorch_forecasting.metrics import QuantileLoss, NormalDistributionLoss

# A multi-quantile head is one object: pass the levels, get the summed pinball loss
# plus built-in monotone decoding and prediction utilities.
quantile_loss = QuantileLoss(quantiles=[0.1, 0.25, 0.5, 0.75, 0.9])

# A parametric Gaussian head is likewise one object: it owns the NLL and sampling.
normal_loss = NormalDistributionLoss()
# In a full model (e.g. TemporalFusionTransformer or DeepAR) you simply pass
# loss=quantile_loss or loss=normal_loss; the network's output layer is sized and
# decoded automatically, and .to_prediction()/.to_quantiles() handle the rest.

Read together, the three blocks make the section concrete. Code 19.3.1 built the quantile family with its signature crossing fix; Code 19.3.2 built the parametric family with its closed-form coherence; Code 19.3.3 showed both are one line each in production. The predictive fan from either model widens with $x$ exactly as the heteroscedastic data demand, which is the whole point: a point forecaster would have drawn a single line through the middle and told you nothing about the spread that every downstream decision needs.

6. Exercises

These exercises build directly on the worked example; the implementation tasks reuse the data generator of Code 19.3.1. Solutions to selected exercises appear in Appendix G.

Conceptual

1. The asymmetry at the optimum. Starting from the expected pinball loss $\mathbb{E}_y[\rho_\tau(y, q)]$, differentiate with respect to $q$ and show its minimizer satisfies $\Pr(y \le q) = \tau$. Explain in one sentence why this means the $\tau = 0.5$ pinball loss recovers the median and squared error recovers the mean.
2. Family choice by support. For each target, name the parametric family from Figure 19.3.2 you would use and state in one phrase why: (a) daily log-returns of an equity index; (b) units of a product sold per day; (c) hours until the next machine failure; (d) rainfall (mm) per day in a desert region with many dry days. Identify the specific failure a Gaussian head would exhibit on (b) and (d).

Implementation

1. Break and fix crossing. Modify Code 19.3.1 to use a naive head that outputs all five quantiles with one unconstrained linear layer (no monotone construction). Train it and count crossings; then add post-hoc sorting at inference and re-count. Compare both against the monotone-by-construction version on (i) number of crossings and (ii) final pinball loss, and explain why sorting changes the loss.
2. Calibrate the fan. Using the trained Gaussian head of Code 19.3.2, compute the empirical coverage of the predicted 10-to-90 interval on a held-out split (fraction of targets falling between $q_{0.1}$ and $q_{0.9}$); it should be near 0.80. Repeat for the quantile head of Code 19.3.1 and compare which head is better calibrated on this heteroscedastic data, and why a misspecified family would show up here.

Open-ended

1. Marginal versus joint, measured. Build a tiny two-step autoregressive forecaster (predict $y_{t+1}$, feed the sample back, predict $y_{t+2}$) with a Gaussian head, and a direct two-output marginal head, on a series with strong positive autocorrelation. For both, estimate the probability that both steps exceed a threshold. Show numerically that the marginal model underestimates this path probability, quantify the gap, and connect the gap to the across-step correlation the autoregressive sampler captures (subsection four).