DeepAR and Probabilistic RNN Forecasting

"They asked me for tomorrow's number and I refused to give just one. Instead I drew a hundred tomorrows: a calm one, a spiked one, three that quietly fail, a long tail of disasters nobody wants to fund inventory for. Then I shrugged and said, the truth was never a number. It was a shape, and here is its weight."

A DeepAR Decoder Sampling a Hundred Possible Tomorrows

Chapters 9 through 13 built the neural machinery of sequence modeling: the recurrent cell and backpropagation through time (Section 9.3), the gated RNNs of Chapter 10, the temporal convolutions of Chapter 11, attention in Chapter 12, and the state-space models of Chapter 13. Every one of those, as we trained it, produced a point: a single predicted value per step under a squared-error or absolute-error loss. This chapter turns that machinery toward forecasting as it is actually practiced, and the first and most consequential turn is the one this section makes. We stop predicting a number and start predicting a distribution. The recurrent backbone is unchanged; what changes is the head that sits on top of the hidden state and the loss that trains it. That single substitution, a likelihood head and a negative-log-likelihood objective in place of a linear head and a mean-squared error, converts every RNN of Part III into a probabilistic forecaster, and it is the conceptual seed of the entire uncertainty-quantification programme of Chapter 19, to which this section repeatedly looks forward. We use the unified notation of Appendix A: $z_t$ for the observed target at step $t$ (DeepAR's convention), $\mathbf{x}_t$ for covariates, $\mathbf{h}_t$ for the hidden state, $\theta_t$ for the per-step distribution parameters.

This section installs five competencies. First, to articulate why a decision needs a distribution and not a point, and what goes wrong when you hand a planner a mean. Second, to construct the DeepAR likelihood head: an RNN whose output layer emits the parameters of a chosen likelihood (Gaussian, Student-t, or negative binomial) and is trained by NLL. Third, to generate multi-step probabilistic forecasts by ancestral sampling, feeding each sampled value back into the recurrence. Fourth, to understand global training across many series, scale handling, and covariates, the three ideas that separate DeepAR from a per-series classical model. Fifth, to evaluate a probabilistic forecast honestly with quantile loss and the continuous ranked probability score. The worked example threads all five together in code.

1. From Point to Probabilistic Forecasting: Predict a Distribution, Not a Number Beginner

A point forecast collapses the future into a single number and, in doing so, throws away exactly the information a decision most needs: how uncertain that number is. Consider the difference between two forecasts that both report an expected demand of one hundred units. In the first, demand is almost certainly between ninety-five and one hundred five; in the second, it could plausibly be anywhere from twenty to one hundred eighty. A point forecast reports "one hundred" in both cases and is, in both cases, useless for the actual decision, which is how many units to stock. The first situation calls for ordering near one hundred; the second, if stockouts are costly, calls for ordering far above one hundred to buy a safety margin. The quantity that drives the order is not the mean of the future but a high quantile of it, and a point forecast does not carry quantiles.

This is the general shape of the argument, and it recurs across every domain in this book. Decisions under uncertainty are almost never functions of the expected outcome alone; they are functions of the whole distribution, usually through an asymmetric cost. Stockout cost differs from holding cost, so the optimal order is a quantile (the classic newsvendor solution, which we revisit in the decision chapters of Part VI). A capacity planner provisions for a peak, not an average, so they need an upper quantile of load. A trading desk sizes positions against a tail loss (value-at-risk is literally a quantile of the loss distribution), not against expected return. In each case the object that the decision consumes is a quantile, an interval, or an expectation of a nonlinear cost, and every one of those is a functional of the predictive distribution that a point forecast simply cannot provide. Predicting a distribution is therefore not a refinement of point forecasting for the statistically fastidious; it is the minimum the decision requires.

There is a second, quieter reason point forecasting under squared error is inadequate for many real series: it implicitly assumes a symmetric, constant-variance, Gaussian-shaped world, and real demand, traffic, and count data violate all three assumptions. Retail demand is a nonnegative integer that is often zero (intermittent demand), heavily right-skewed, and has a variance that grows with its mean. Squared-error regression on such data produces forecasts that can go negative, that ignore the skew, and that report a uniform uncertainty band when the true uncertainty scales with volume. The probabilistic framing fixes this at the root by letting you choose the likelihood to match the data: a negative binomial for skewed nonnegative counts, a Student-t for heavy-tailed continuous series, a Gaussian for well-behaved real values. The forecaster then learns the parameters of that correctly-shaped distribution rather than forcing a mismatched one. We forward-reference the full treatment of calibration, coverage, and conformal guarantees to Chapter 19; here the point is simply that the distribution is the deliverable.

2. DeepAR: An Autoregressive RNN That Outputs Likelihood Parameters Intermediate

DeepAR keeps the autoregressive recurrent backbone of Part III intact and changes only the top of the network. At each step the RNN consumes the previous observed value and the current covariates, updates its hidden state, and then, instead of emitting a single predicted value, emits the parameters of a probability distribution over the next value. Concretely, with an LSTM or GRU cell producing the hidden state $\mathbf{h}_t$, the model computes

where $[z_{t-1}, \mathbf{x}_t]$ concatenates the previous target with the current covariates (lagged-target autoregression), and $\theta(\cdot)$ is a small affine-plus-link readout that maps the hidden state to the parameters $\theta_t$ of the chosen likelihood $p(z \mid \theta_t)$. For a Gaussian likelihood, $\theta_t = (\mu_t, \sigma_t)$, produced by

where the softplus link $\mathrm{softplus}(x) = \log(1 + e^x)$ guarantees a positive standard deviation. The mean head is unconstrained; the scale head is passed through softplus precisely so the network can never emit a negative or zero variance. This positivity-link pattern is the general recipe: every distribution parameter that must be positive (the $\sigma$ of a Gaussian, the degrees of freedom and scale of a Student-t, the dispersion of a negative binomial) is produced by passing an unconstrained linear output through a link that maps the real line to the positive half-line. Figure 14.1.1 shows the head sitting atop the recurrent backbone.

Training maximizes the likelihood of the observed sequence under the model, equivalently minimizing the negative log-likelihood (NLL). Because the model is autoregressive and conditioned on the true past during training (teacher forcing), the joint likelihood factorizes into a product of per-step likelihoods, and the loss is the sum of per-step negative log-densities:

For the Gaussian head this is the familiar Gaussian negative log-likelihood, written out per step as

Read this loss carefully, because its structure is the entire reason probabilistic training behaves differently from point training. The third term, $(z_t - \mu_t)^2 / (2\sigma_t^2)$, is a squared error scaled by the predicted variance: the model is rewarded for putting its mean near the target, but the reward is divided by how much variance it claimed. The second term, $\log \sigma_t$, is a penalty for claiming large variance. Together they create a tension the network must resolve: it cannot simply inflate $\sigma_t$ to make the squared-error term small, because $\log \sigma_t$ punishes that; and it cannot shrink $\sigma_t$ toward zero to minimize $\log \sigma_t$, because then the squared-error term explodes whenever the mean misses. The minimum of this tension is reached when $\sigma_t$ equals the true conditional standard deviation of the residual. This is why the Gaussian NLL does not merely fit a mean; it fits a calibrated variance, learning to report large uncertainty where the data are noisy and small uncertainty where they are predictable. A plain mean-squared-error loss is exactly this NLL with $\sigma_t$ frozen at a constant, which is why MSE forecasters cannot express heteroscedastic uncertainty.

The likelihood is a modeling choice matched to the data, and DeepAR's three canonical choices cover the common cases. The Gaussian, with parameters $(\mu, \sigma)$, suits real-valued, roughly symmetric series. The Student-t, adding a degrees-of-freedom parameter $\nu$, gives heavier tails for series with occasional large shocks (financial returns, the Chapter 5 volatility series). The negative binomial, with a mean and a dispersion parameter, is the workhorse for nonnegative integer count data (retail unit sales, web traffic), where it captures the right skew and the mean-dependent variance that a Gaussian cannot. The negative-binomial mean and the Gaussian mean are read out the same way; only the second parameter and the density formula change. Selecting the likelihood is the single most important modeling decision in DeepAR, and getting it wrong (a Gaussian on intermittent counts) produces poorly calibrated intervals no amount of training can fix.

To turn the trained one-step likelihood into a multi-step forecast, DeepAR uses ancestral sampling (also called the autoregressive rollout). During training the RNN was fed the true previous value; at prediction time, beyond the end of the observed history, there is no true value to feed, so the model feeds back a sample. Starting from the hidden state at the end of the conditioning range, it draws $\hat{z}_{t+1} \sim p(z \mid \theta_{t+1})$ from the predicted distribution, then feeds that drawn value in as the input for step $t+2$, predicts $\theta_{t+2}$, draws $\hat{z}_{t+2}$, and continues to the forecast horizon. One such rollout produces a single sampled future trajectory. Repeating the rollout many times (typically a few hundred) produces an ensemble of trajectories, and the empirical distribution of those trajectories at each future step is the model's predictive distribution: its quantiles give prediction intervals, its mean gives a point forecast if one is wanted. The sampling is "ancestral" because each step's distribution is conditioned on the actual samples drawn at all earlier steps, so the trajectories correctly capture how uncertainty compounds the further out you forecast.

3. Global Training: One Model, Many Series, With Scale and Covariates Intermediate

The classical models of Part II were fit per series: one ARIMA for one product, with its own parameters estimated from that product's history (Chapter 5). For a catalog of a hundred thousand items this is a hundred thousand separate fits, each starved of data, each unable to forecast a freshly launched item with no history at all. DeepAR makes the opposite bet: train one model globally across all the series at once, so that a single set of RNN weights is shared by every series and learns patterns common to the whole population. This is the most important practical departure from classical forecasting, and it is what the "AR" in a deep model buys you over a per-series autoregression.

Global training delivers three concrete advantages. It lets the model borrow statistical strength across series, learning seasonality, promotion responses, and lifecycle shapes from the whole panel rather than re-estimating them from each thin history. It enables cold-start forecasting: a new product with zero observed sales can still be forecast from its covariates (category, price, launch week) because the model has seen thousands of comparable launches. And it amortizes into a single training and serving artifact rather than a hundred thousand models to fit, store, and monitor. The cost is that the model must be expressive enough to represent every series with one parameter set, which is exactly what a deep recurrent network, conditioned on per-series covariates, is built to do.

Global training across series of wildly different magnitudes raises an immediate problem: scale. One product sells two units a week, another sells fifty thousand; a single network with shared weights cannot natively span five orders of magnitude in its inputs and outputs without the large series dominating the loss and the small series being ignored. DeepAR handles this with a per-series scale factor, computed as a simple statistic of each series' history (in the original paper, $\nu_i = 1 + \frac{1}{t_0}\sum_{t=1}^{t_0} z_{i,t}$, roughly the mean level of series $i$). The model divides each series' inputs by its scale before the RNN sees them, predicts in the normalized space, and multiplies the predicted distribution's parameters back by the scale at the output. The network therefore only ever has to learn the shape of demand, not its absolute magnitude, and the same weights serve the two-unit item and the fifty-thousand-unit item. This scaling is the unglamorous engineering that makes global training actually work on real catalogs, and forgetting it is a common cause of a DeepAR that trains but forecasts the small series as flat zero.

The third ingredient is covariates: time-varying and static features concatenated into the RNN input alongside the lagged target. Time-varying covariates known into the future (day of week, month, holiday flags, scheduled promotions, price) let the model condition its forecast on what it knows about the future, which is how it forecasts a holiday spike it has never seen for this item but has seen for others. Static covariates (product category, region, an embedded item identifier) let one global model specialize its behavior per series without separate parameters. Because the future covariates are fed at every prediction step, the autoregressive rollout of subsection two carries them forward naturally: at each sampled step the model sees the real future calendar and price even as it samples the unknown demand. Covariates are what let a single global model be, in effect, a different forecaster for every series and every future date, all from one shared weight set.

4. Evaluating Probabilistic Forecasts: Quantile Loss and CRPS Advanced

A probabilistic forecast cannot be scored by mean-squared error, because MSE only sees the point and is blind to whether the claimed uncertainty was honest. Evaluating a distribution requires proper scoring rules, metrics minimized in expectation only by the true predictive distribution, so that a model cannot improve its score by lying about its uncertainty. Two such metrics dominate forecasting practice, and both reappear in full generality in Chapter 19; we introduce them here as the tools that make the worked example's evaluation meaningful.

The quantile loss (also called the pinball loss) scores a single predicted quantile against the realized value. For target quantile level $\rho \in (0,1)$, predicted quantile $\hat{q}_\rho$, and observation $z$, it is

The asymmetry is the point: for a high quantile like $\rho = 0.9$, under-predicting (the truth lands above your quantile) is penalized at weight $0.9$ while over-predicting is penalized at weight $0.1$, so the loss is minimized exactly when the predicted $\hat{q}_{0.9}$ sits at the true ninetieth percentile, with ten percent of observations above it. Averaging the quantile loss over a grid of levels (often the deciles $0.1, 0.2, \dots, 0.9$) gives a summary of how well the whole set of predicted quantiles is calibrated, and this averaged quantile loss is the standard reported metric on the retail benchmarks DeepAR targets.

The continuous ranked probability score (CRPS) scores the entire predictive distribution against the observation in one number. With predictive CDF $F$ and observation $z$, it is the integrated squared difference between the predicted CDF and the step function that jumps at the observation,

which rewards a forecast for concentrating probability mass near the truth while remaining proper. CRPS generalizes the absolute error to distributions: if the forecast is a point mass (a degenerate distribution at a single value), CRPS reduces exactly to the absolute error $|z - \hat{z}|$, so it is the natural distributional analogue of mean-absolute error and is reported in the same units as the data. In practice, when the forecast is represented by samples (as DeepAR's rollout produces), CRPS is estimated from the empirical CDF of the samples, or equivalently from the averaged quantile loss over a fine quantile grid, which is why the two metrics travel together. A model with a lower CRPS has, on average, a predictive distribution genuinely closer to the truth, sharpness and calibration both accounted for, and that is the single number to track when comparing probabilistic forecasters.

5. Worked Example: A Gaussian-Head Probabilistic RNN, From Scratch and From a Library Advanced

We now build the whole pipeline. The plan: implement a small DeepAR-style probabilistic RNN from scratch in PyTorch, a GRU backbone with a Gaussian likelihood head trained by NLL; train it on a synthetic seasonal series; generate a multi-step probabilistic forecast by ancestral sampling; plot the prediction intervals; and compute the NLL and quantile loss. Then we reproduce the same model in a handful of lines with the GluonTS and neuralforecast DeepAR and state the line-count reduction. Code 14.1.1 is the from-scratch model and its training loop.

import torch, torch.nn as nn, math

torch.manual_seed(0)

# --- synthetic data: a seasonal series with level + noise (one global panel of series) ---
def make_series(n=600, period=24):
    t = torch.arange(n, dtype=torch.float32)
    season = 10.0 * torch.sin(2 * math.pi * t / period)      # daily-style seasonality
    level  = 50.0 + 0.02 * t                                  # slow upward trend
    noise  = 3.0 * torch.randn(n)                             # observation noise
    return (level + season + noise)

z = make_series()
ctx, horizon = 480, 24                                        # conditioning range, forecast horizon

class DeepARGaussian(nn.Module):
    """Autoregressive GRU with a Gaussian likelihood head (mu, sigma)."""
    def __init__(self, hidden=32):
        super().__init__()
        self.rnn   = nn.GRU(input_size=1, hidden_size=hidden, batch_first=True)
        self.mu    = nn.Linear(hidden, 1)                     # unconstrained mean head
        self.presig = nn.Linear(hidden, 1)                    # pre-softplus scale head

    def forward(self, x, h=None):                            # x: (batch, time, 1) lagged target
        out, h = self.rnn(x, h)
        mu    = self.mu(out).squeeze(-1)                      # (batch, time)
        sigma = torch.nn.functional.softplus(self.presig(out)).squeeze(-1) + 1e-3  # > 0
        return mu, sigma, h

def gaussian_nll(z_true, mu, sigma):                         # per-step Gaussian NLL of subsection 2
    return (0.5 * math.log(2 * math.pi) + torch.log(sigma)
            + (z_true - mu) ** 2 / (2 * sigma ** 2)).mean()

model = DeepARGaussian()
opt = torch.optim.Adam(model.parameters(), lr=1e-2)

# teacher-forced training on the conditioning range: input is z shifted by one step
zc = z[:ctx]
inp = zc[:-1].view(1, -1, 1)                                  # z_{t-1} as input
tgt = zc[1:].view(1, -1)                                      # z_t as target
for epoch in range(400):
    opt.zero_grad()
    mu, sigma, _ = model(inp)
    loss = gaussian_nll(tgt, mu, sigma)                      # minimize NLL = maximize likelihood
    loss.backward(); opt.step()
print("final train NLL = %.4f" % loss.item())

final train NLL = 2.5631

With the model trained, we generate a probabilistic forecast by the ancestral sampling of subsection two: roll the GRU forward past the conditioning range, at each step sample a value from the predicted Gaussian, and feed that sample back as the next input. Code 14.1.2 runs many such rollouts and reads prediction intervals and the evaluation metrics from the sample ensemble.

import numpy as np

@torch.no_grad()
def sample_paths(model, history, horizon, n_samples=300):
    """Ancestral sampling: feed each drawn sample back in (subsection 2 rollout)."""
    # warm up the hidden state on the full conditioning history
    _, _, h = model(history.view(1, -1, 1))
    last = history[-1].view(1, 1, 1).repeat(n_samples, 1, 1)  # start from last observed value
    h = h.repeat(1, n_samples, 1)                             # one hidden state per sample path
    paths = []
    for _ in range(horizon):
        mu, sigma, h = model(last, h)                         # one-step predictive Gaussian
        draw = mu + sigma * torch.randn_like(mu)             # sample z_hat ~ N(mu, sigma)
        paths.append(draw.squeeze(-1))
        last = draw.view(n_samples, 1, 1)                     # feed the SAMPLE back in
    return torch.stack(paths, dim=1).numpy()                  # (n_samples, horizon)

paths = sample_paths(model, z[:ctx], horizon)                # 300 sampled future trajectories
truth = z[ctx:ctx + horizon].numpy()

q10, q50, q90 = np.percentile(paths, [10, 50, 90], axis=0)   # prediction-interval quantiles

def pinball(z_true, q_pred, rho):                            # quantile (pinball) loss, subsection 4
    e = z_true - q_pred
    return np.mean(np.maximum(rho * e, (rho - 1) * e))

# Gaussian NLL of the held-out truth under the sample mean/std (sharper test than train NLL)
mu_h, sd_h = paths.mean(0), paths.std(0) + 1e-3
nll_test = np.mean(0.5*np.log(2*np.pi) + np.log(sd_h) + (truth-mu_h)**2/(2*sd_h**2))
print("test NLL          = %.4f" % nll_test)
print("pinball q90       = %.4f" % pinball(truth, q90, 0.9))
print("90%% interval cov. = %.2f" % np.mean((truth >= q10) & (truth <= q90)))

test NLL          = 2.9214
pinball q90       = 1.1077
90% interval cov. = 0.83

Figure 14.1.2 shows what the sample ensemble looks like as a forecast: a median line wrapped in a shaded interval that fans out with the horizon, the visual signature of a probabilistic forecast that a point forecaster cannot produce.

The from-scratch model, its training loop, the sampling routine, and the metric helpers ran to roughly seventy lines. A production library collapses the architecture, the scale handling, the covariate plumbing, the multi-likelihood support, and the sampling-based evaluation to a few lines, which is the "Right Tool" point made concrete.

# Nixtla neuralforecast: global probabilistic DeepAR in ~8 lines
from neuralforecast import NeuralForecast
from neuralforecast.models import DeepAR
from neuralforecast.losses.pytorch import DistributionLoss

model = DeepAR(h=24, input_size=48,
               loss=DistributionLoss(distribution="Normal", level=[80, 90]),
               max_steps=400)
nf = NeuralForecast(models=[model], freq="H")
nf.fit(df=panel_df)                          # df: many series stacked (global training)
fc = nf.predict()                            # returns median + 80/90% prediction intervals

Read the three code blocks together as the section in miniature. Code 14.1.1 built the likelihood head and trained it by NLL, the substitution of subsection two. Code 14.1.2 turned the trained one-step model into a multi-step probabilistic forecast by ancestral sampling and scored it with the proper losses of subsection four. The library shortcut showed that everything we built by hand is one constructor argument away in production, and that switching likelihoods to match the data (subsection two) or training globally across a panel (subsection three) is a matter of passing a dataframe of stacked series rather than rewriting the model. The deliverable, in every version, is a distribution.

Exercises

The exercises split into conceptual, implementation, and open-ended. Attempt the conceptual ones before reading further; the implementation ones extend the code of subsection five; the open-ended one points toward research.

1. Conceptual. A colleague trains a point RNN under mean-squared error and proposes to "add uncertainty" afterward by reporting the historical residual standard deviation as a constant band around every forecast. Explain, using the Gaussian NLL of subsection two, why this constant-band forecast is exactly the special case of DeepAR with $\sigma_t$ frozen, and give one concrete situation (name the series property) where it produces badly miscalibrated intervals that a learned $\sigma_t$ would fix.
2. Implementation. Modify Code 14.1.1 to use a negative-binomial likelihood instead of Gaussian: replace the two-headed output with a mean head (softplus link, positive) and a dispersion head (softplus link), implement the negative-binomial NLL, and retrain on a synthetic count series (round the seasonal signal and clip at zero). Compare the eighty-percent interval coverage of the Gaussian model and the negative-binomial model on the count data and report which is better calibrated and why, connecting your answer to the likelihood-choice discussion of subsection two.
3. Open-ended. The ancestral sampling of Code 14.1.2 produces correlated trajectories whose uncertainty compounds with the horizon, but it samples each step's value independently given the history. Investigate one limitation of this scheme for multivariate forecasting (forecasting several correlated series jointly) and sketch how a research approach you find in the literature (for example a low-rank Gaussian copula head, as in the GluonTS DeepVAR, or a normalizing-flow output) would model the cross-series dependence that independent per-series DeepAR sampling misses. Connect your sketch forward to the probabilistic-forecasting treatment of Chapter 19.