Under active development. A stable version is expected in November 2026. Feedback is welcome by email or via GitHub issues.

9 Distributional Neural Networks

9.1 Overview

In many economic and financial applications, the signal-to-noise ratio is low, and decisions depend on the entire predictive distribution (also called a distribution forecast), not just a single point estimate. A central bank may care about downside inflation risk, a financial economist may care about the left tail of returns, and a risk manager may care about the probability of a large loss. A neural network (NN) trained with mean squared error (MSE) can be useful for estimating a conditional mean, \mathbb{E}[Y \mid \mathbf{X}=\mathbf{x}], but it does not by itself provide a calibrated predictive distribution.

Distributional neural networks address this by reframing the prediction problem. Instead of predicting a single number, the network predicts the parameters of a full conditional distribution. This connects the neural-network chapters back to the earlier chapters on information theory and evaluating predictive distributions: likelihood-based training corresponds to the log score, and alternative proper scoring rules such as the continuous ranked probability score (CRPS) can also be used when we want different sensitivity to tails or outliers. The broader literature on probabilistic forecasting and its evaluation is reviewed in Gneiting and Katzfuss (2014).

9.2 Roadmap

We first contrast point-prediction networks with distributional networks.
We then connect negative log-likelihood training to cross-entropy, Kullback–Leibler (KL) divergence, and proper scoring rules.
Next, we study Gaussian distributional networks and the instability that arises when mean and variance are trained jointly.
We use Hemisphere Neural Networks (HNNs)—networks with a shared representation and separate mean and variance branches—as a case study for stabilizing distributional estimation in macro-financial settings.
We then turn to Mixture Density Networks (MDNs), which output the weights and component parameters of a mixture distribution, for richer conditional distributions.
We close with an empirical application—Gaussian density forecasts for equity realized variance, evaluated with proper scoring rules on a held-out test block—followed by key takeaways, common pitfalls, and exercises.

9.3 From Point Forecasts to Predictive Distributions

As discussed in the Feed-Forward Neural Networks chapter, a standard neural network trained with mean squared error produces a point forecast:

Traditional NN Output: \hat{y} = f_\theta(\mathbf{x})
Distributional NN Output: p_\phi(y \mid \mathbf{x}) = p(y \mid \boldsymbol{\eta}_\phi(\mathbf{x})), where \boldsymbol{\eta}_\phi(\mathbf{x}) are the distribution parameters learned by the network.

The distributional version replaces the scalar output by one or more output heads that map features into distributional parameters. For example, a Gaussian distributional network can produce a conditional mean and variance, while a Student-t network can additionally produce a degrees-of-freedom parameter. For the Gaussian case the parameter vector is \boldsymbol{\eta}_\phi(\mathbf{x})=(\mu(\mathbf{x}),\log\sigma^2(\mathbf{x})) (equivalently (\mu(\mathbf{x}),\sigma^2(\mathbf{x}))), so the scalar heads \mu_t and \sigma_t^2 that appear later are recognizably its components.

The model is typically trained by maximizing the log-likelihood of the data, which is equivalent to minimizing the negative log-likelihood (NLL) loss:

L(\phi) = -\frac{1}{N}\sum_{i=1}^N \log p(y_i \mid \mathbf{x}_i; \boldsymbol{\eta}_\phi(\mathbf{x}_i)).

Here, y_i is the observed outcome, \mathbf{x}_i is the feature vector, \boldsymbol{\eta}_\phi(\mathbf{x}_i) denotes the distributional parameters implied by the neural network, and \phi collects all trainable network weights. We write the objective as a sample average, matching the loss convention of the network chapters; the 1/N does not affect the minimizer.

9.4 Training Scores and Information Theory

Information-theoretic interpretation

This NLL objective is the cross-entropy loss encountered in the likelihood and KL-divergence discussion: maximum likelihood estimation (MLE) minimizes sample cross-entropy. At the population level, expected conditional NLL differs from the expected conditional Kullback–Leibler (KL) divergence from the true conditional distribution to the model by a term that does not depend on the model parameters. Provided the relevant expectations are finite, both criteria therefore have the same population minimizer: the KL projection of the true conditional distribution onto the represented model family. Convergence of the empirical estimator to that target additionally requires identification within the density family, dependence-appropriate control of model complexity, a uniform law of large numbers, and sufficiently accurate optimization.

The score and the model family play different roles. A realized logarithmic score (LogS) is local: it evaluates only the density assigned to the outcome that occurred. Its expected value is strictly proper and is uniquely minimized by the true density when that density belongs to the admissible class and the expectations are finite. A Gaussian network can adjust only its represented conditional mean and variance functions, so its population target under misspecification is the closest represented Gaussian in KL divergence. Student-t and mixture outputs enlarge the family of distributions to which the expected score can respond.

The same logic can be stated in the language of proper scoring rules. The negative log-likelihood corresponds to LogS: after observing y_i, the model is rewarded for assigning high density to the realized outcome. Minimizing average NLL is therefore the training analogue of minimizing the out-of-sample log score; see the discussion of LogS and CRPS in the Evaluating Predictive Distributions chapter.

Other proper scoring rules can also be used as training losses if they are computable and differentiable for the chosen predictive family. For example, the CRPS compares the full predictive cumulative distribution function (CDF) to the realized outcome. The choice of training score should match the forecasting object: log-likelihood is very sensitive to tail misspecification and extreme observations, while CRPS often behaves more like an integrated distributional error.

Question for Reflection

You fit a Gaussian distributional network by maximum likelihood for daily stock returns. The in-sample log-likelihood is high, but the realized probability integral transform (PIT) histogram, computed as in the predictive-distribution chapter, is U-shaped. What does the U-shape suggest, and why does it not establish that conclusion by itself? Name two model failures that are both compatible with it, and propose one additional diagnostic that helps separate them.

9.5 A Case Study: Hemisphere Neural Networks (HNN)

A key challenge in distributional modeling is reliably estimating multiple distribution parameters at once. Consider the Gaussian case where we want to learn both the conditional mean and variance:

Y_{t+1}\mid\mathbf{X}_t=\mathbf{x}_t \sim \mathcal{N}\big(\mu_t, \sigma_t^2\big),\quad \mu_t = f_{\phi_\mu}(\mathbf{x}_t),\ \ \sigma_t^2 = g_{\phi_\sigma}(\mathbf{x}_t)>0,

where Y_{t+1} is the random outcome and \mathbf{X}_t the random predictor vector with realization \mathbf{x}_t, following the book’s case convention. Throughout this chapter, conditioning on \mathbf{x}_t abbreviates conditioning on the event \mathbf{X}_t=\mathbf{x}_t. The subscript t on \mu_t and \sigma_t^2 records the forecast origin: both are functions of \mathbf{x}_t and describe the conditional distribution of Y_{t+1}, consistent with the two-index convention \hat{y}_{t+1\mid t} used elsewhere in the book. The per-time-step negative log-likelihood loss is, up to an additive constant,

\ell_{t+1}(\phi_\mu,\phi_\sigma)= \frac{1}{2}\log g_{\phi_\sigma}(\mathbf{x}_t) + \frac{(y_{t+1}-f_{\phi_\mu}(\mathbf{x}_t))^2}{2\,g_{\phi_\sigma}(\mathbf{x}_t)},

the form analyzed in Exercise 9.1. Econometricians have met this objective before: it is exactly the Gaussian quasi-likelihood of a generalized autoregressive conditional heteroskedasticity (GARCH) model, with the parametric variance recursion replaced by the network map g_{\phi_\sigma}(\mathbf{x}_t). The quasi-maximum-likelihood theory of Bollerslev and Wooldridge (1992) also explains why this objective is a sensible target even when the true conditional distribution is not Gaussian: the criterion identifies the conditional mean and the conditional variance as long as those two moments are correctly specified.

In overparameterized models, jointly training the mean and variance heads by maximum likelihood can be unstable. The model can achieve a low loss by letting the mean network overfit, which drives residuals and therefore the fitted variance toward zero, or by letting the variance network absorb too much variation. This is a distributional overfitting problem: the model has to learn both the location and the uncertainty, and the two tasks can interfere with each other.

The HNN architecture, proposed by Goulet Coulombe, Frenette, and Klieber (2026), is designed to stabilize this process.

HNN architecture. The architecture splits the network after a shared common core into two specialized hemispheres. The core is a set of shared layers that learn a common representation of the inputs; because latent drivers can influence both moments through this shared representation, the two heads draw on common predictors—a representation-level link between mean and variance. Genuine autoregressive conditional heteroskedasticity (ARCH)-type dynamics, in which the variance is driven by past innovations, or volatility-in-mean effects, in which the conditional variance enters the mean equation, do not follow from shared layers alone: they arise only if the corresponding variables—lagged residuals, lagged variance, or an explicit cross-head connection—enter the model, as in the Temporal Dynamics extension below. The mean hemisphere is a dedicated stack of layers ending in a linear output head that predicts \mu_t. The variance hemisphere mirrors it with a softplus output head—the softplus transform \zeta(z)=\log(1+e^{z})>0 (introduced in the Feed-Forward Neural Networks chapter) maps any real pre-activation to a positive value—so that the predicted variance \sigma_t^2 is always positive.

The resulting architecture is illustrated in Figure 9.1.

Figure 9.1: Architecture of a Hemisphere Neural Network. The input vector \mathbf{x}_t (top) passes through a shared common core of L_c layers, after which the network splits into two dedicated branches: a mean hemisphere of L_m layers ending in a linear head that outputs \mu_t = f_{\phi_\mu}(\mathbf{x}_t) (lower left), and a variance hemisphere of L_v layers ending in a softplus head that outputs \sigma_t^2 = g_{\phi_\sigma}(\mathbf{x}_t) > 0 (lower right). Arrows show the flow of computation.

When \mathbf{x}_t contains relevant lagged shocks or variance measures, the network can learn reactive volatility in which variance rises after shocks, as in GARCH-type dynamics. With genuine leading indicators known at the forecast origin, it can instead learn proactive volatility in which predicted variance rises before the realization. The econometric intuition is that some variables may move the conditional mean, some may move conditional uncertainty, and some may affect both.

Stabilizing joint MLE in practice.

The HNN framework introduces three key ingredients to discipline the unstable MLE procedure:

Volatility-emphasis constraint. If Y_{t+1} has a finite second moment, the population law of total variance gives
$\operatorname{Var}(Y_{t+1}) = \mathbb{E}\bigl[\operatorname{Var}(Y_{t+1}\mid \mathbf{X}_t)\bigr] + \operatorname{Var}\bigl(\mathbb{E}[Y_{t+1}\mid \mathbf{X}_t]\bigr),$
where the outer expectation and variance average over the random predictor vector \mathbf{X}_t. Thus, the expected conditional variance is strictly less than the unconditional variance whenever the conditional mean varies with the predictors. Under stationarity and ergodicity, time averages can estimate these population moments. In finite samples, HNN uses the corresponding empirical relation as a regularization device:
$\frac{1}{T}\sum_{t=1}^{T} g_{\phi_\sigma}(\mathbf{x}_t) \approx \kappa\cdot \operatorname{Var}(Y_{t+1}), \qquad 0 < \kappa < 1.$
The multiplier \kappa is calibrated from out-of-bag residuals of a non-distributional baseline—residuals computed for each observation from fits whose training draws exclude it, the device used in the empirical application of the predictive-distribution chapter and treated fully in the random forests chapter. The constraint regularizes the variance hemisphere away from zero and prevents it from absorbing predictable variation that should be explained by the mean hemisphere.
Blocked subsampling (time-aware bagging). Blocked subsampling trains the model repeatedly on different time blocks (“bags”). It adapts bagging—short for bootstrap aggregation—to time-series data; the random forests chapter introduces bagging in full. Predictions for \mu_t and \sigma_t^2 are averaged only across fits for which time t was held out of the training sample. Averaging reduces variation across training samples and network initializations.
Blocked out-of-bag recalibration. After estimation, the fitted variance path is recalibrated against out-of-bag squared forecast errors—errors from fits that did not see the corresponding observations—so that predicted uncertainty is disciplined by prediction errors rather than by in-sample residuals (Goulet Coulombe, Frenette, and Klieber 2026). Mechanically, the log out-of-bag squared errors are regressed on the log fitted variance, so the recalibration adjusts both the level and, through the slope, the amplitude of the variance path, damping its variation when the raw variance head has overfit. This recalibration is distinct from the subsampling in ingredient 2 that generates the out-of-bag predictions.

Together, these three ingredients provide the stabilization: dedicated output heads alone do not prevent joint-likelihood degeneracy.

Beyond Gaussian outputs. The hemisphere approach can be extended to a location-scale Student-t output by adding a third hemisphere for the degrees-of-freedom parameter, \nu. Write Y=\mu+sZ with Z\sim t_\nu. The density is defined for \nu>0, its mean is finite only for \nu>1, and its variance is finite only for \nu>2. Requiring \nu>2, for instance through \nu=2+\zeta(z), guarantees finite variance, but s^2 is not that variance:

\operatorname{Var}(Y)=s^2\frac{\nu}{\nu-2}.

If an existing variance head v(\mathbf{x}_t)—the role played by g_{\phi_\sigma}(\mathbf{x}_t) above—is meant to retain the interpretation of a conditional variance, the Student-t scale must satisfy s^2(\mathbf{x}_t)=v(\mathbf{x}_t)(\nu-2)/\nu.

Temporal dynamics. Lagged outcomes, innovations, or squared residuals can enter \mathbf{x}_t to capture GARCH-type effects directly. For more complex sequences, the common core could be a Long Short-Term Memory (LSTM) network.

9.6 Handling Complex Distributions: Mixture Density Networks (MDNs)

For distributions that are skewed, multimodal, or heavier-tailed over empirically relevant ranges than a single Gaussian, one parametric component can be too restrictive. MDNs, proposed by Bishop (1994), address this by modeling the conditional distribution as a mixture of simpler ones, such as Gaussians:

p(y \mid \mathbf{x})=\sum_{k=1}^K \pi_k(\mathbf{x})\,\mathcal{N}\!\big(y;\,\mu_k(\mathbf{x}),\,\sigma_k^2(\mathbf{x})\big).

A finite Gaussian mixture can be substantially more leptokurtic than one Gaussian over the range relevant for an application, but its asymptotic tails remain Gaussian and are dominated by the component with the largest variance. It is therefore not fat-tailed in the strict power-law sense.

The network has three output heads to predict the parameters for all K components:

Mixing Weights \pi_k(\mathbf{x}): A softmax layer over the K components, \pi_k(\mathbf{x})=\exp(a_k(\mathbf{x}))/\sum_{j=1}^K \exp(a_j(\mathbf{x})). Here, k = 1, \ldots, K indexes the mixture components and K is their number; the transform yields \pi_k>0 and \sum_{k=1}^K\pi_k=1.
Means \mu_k(\mathbf{x}): A linear layer.
Standard Deviations \sigma_k(\mathbf{x}): A softplus layer ensures positivity.

MDNs are particularly useful for modeling phenomena whose predictive distributions are multimodal. Examples include:

Multimodal inflation distributions (e.g., a low-inflation versus a high-inflation regime).
Asset returns with regime-switching behavior.
Labor market dynamics that might feature multiple equilibria.

MDN training has two transparent difficulties (Bishop 1994). First, component labels are not identified and a subset of components can receive almost all the weight. Second, the likelihood is unbounded if a sufficiently flexible component centers on an observation while its scale approaches zero and its mixing weight does not vanish quickly enough. Careful initialization, penalties on extreme mixing weights, and a positive scale floor are possible safeguards, although the latter two change the criterion or parameter space.

A different route to flexible distributional output is to train separate network heads with the quantile, or pinball, loss

\rho_\tau(y-q)=(y-q)\bigl(\tau-\mathbf{1}\{y<q\}\bigr), \qquad \tau\in(0,1),

at several quantile levels. The conditional expected loss is minimized by a conditional \tau-quantile q_\tau(\mathbf{x}). Unless the architecture or loss imposes an ordering restriction, the estimated quantiles can cross. A parametric density head avoids crossing because all quantiles come from one cumulative distribution function.

Econometric Interpretation

An MDN is a flexible way to approximate a conditional density with latent regimes. It should not automatically be interpreted as a structural regime-switching model: the mixture components are learned predictive components unless additional economic restrictions are imposed.

9.7 Empirical Application: Distributional Network for SPY Realized Variance

The realized-variance illustration of the previous chapter produced point forecasts of daily log realized variance, \log RV_{t+1}, for SPY, the S&P 500 exchange-traded fund, and evaluated them with squared and absolute losses. A Gaussian predictive distribution also requires a conditional variance. We therefore compare two such distributions using the same SPY data, log-transformed heterogeneous autoregressive (HAR) features, and the final 20% test block. Because all hyperparameters are fixed in advance, this chapter combines the previous chapter’s training and validation blocks into one estimation sample.¹

The HAR-OLS Gaussian baseline fits the HAR components by ordinary least squares (OLS) and uses the resulting linear conditional mean. Specifically, it regresses \log RV_{t+1} on \log RV_t, \log\big(\tfrac{1}{5}\sum_{j=0}^{4} RV_{t-j}\big), and \log\big(\tfrac{1}{21}\sum_{j=0}^{20} RV_{t-j}\big). Its conditional variance is the in-sample residual variance from the same fit, computed with the maximum-likelihood convention by dividing the sum of squared residuals by the number of observations. The variance stays constant over the test block, so the baseline is homoskedastic in log-RV space.

The Gaussian distributional HAR-FNN is a feed-forward neural network (FNN) using the same three log-HAR inputs. One hidden layer with 32 softplus units feeds a linear mean head \mu_\phi(\mathbf{x}_t) and an unconstrained variance pre-activation a_\phi(\mathbf{x}_t). The predictive variance is \sigma_\phi^2(\mathbf{x}_t)=c+\operatorname{softplus}(a_\phi(\mathbf{x}_t)) for a small floor c>0. This smooth parameterization avoids the zero-gradient region created by hard clipping with \max\{\exp(a),c\}. We write a_\phi for this scalar head to distinguish it from the full parameter vector \boldsymbol{\eta}_\phi of Section 9.3; here \boldsymbol{\eta}_\phi(\mathbf{x}_t) = (\mu_\phi(\mathbf{x}_t), a_\phi(\mathbf{x}_t)). The heads share the hidden representation, in the spirit of the HNN common core in Section 9.5, but the model does not use the volatility-emphasis constraint, blocked subsampling, or out-of-bag recalibration. It minimizes Gaussian negative log-likelihood and is evaluated with LogS, CRPS, and the probability integral transform (PIT) on the held-out test block, as in the empirical application of the predictive-distribution chapter.

Show the code

import os
os.environ["KERAS_BACKEND"] = "tensorflow"
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"

import random
import numpy as np
import pandas as pd
import keras
from keras import layers, ops
from sklearn.linear_model import LinearRegression
from scipy.stats import norm
import properscoring as ps
import matplotlib.pyplot as plt

random.seed(0)
np.random.seed(0)
keras.utils.set_random_seed(0)

rv_source = (
    pd.read_csv("data/taq_spy/SPY_daily_measures.csv", parse_dates=["date"])
    .sort_values("date")
    .loc[:, ["date", "rv"]]
)
rv = rv_source.dropna(subset=["rv"])
rv = rv.loc[rv["rv"] > 0].reset_index(drop=True)

usable = rv_source["date"] <= rv["date"].max()
assert (rv_source.loc[usable, "rv"] > 0).all(), "gap inside the usable sample"

window = 22

def make_log_har(rv_arr):
    X, y, idx = [], [], []
    for end in range(window - 1, len(rv_arr) - 1):
        X.append(np.log([
            rv_arr[end],
            rv_arr[end - 4 : end + 1].mean(),
            rv_arr[end - 20 : end + 1].mean(),
        ]))
        y.append(np.log(rv_arr[end + 1]))
        idx.append(end + 1)
    return np.asarray(X), np.asarray(y), np.asarray(idx)

X, y, target_idx = make_log_har(rv["rv"].to_numpy())
dates = rv.loc[target_idx, "date"].reset_index(drop=True)

n = len(y)
split_valid = int(0.80 * n)
X_trva, y_trva = X[:split_valid], y[:split_valid]
X_te, y_te = X[split_valid:], y[split_valid:]
dates_te = dates.iloc[split_valid:].reset_index(drop=True)

x_center = X_trva.mean(axis=0)
x_scale = X_trva.std(axis=0)
y_center = float(y_trva.mean())
y_scale = float(y_trva.std())

Xs_trva = (X_trva - x_center) / x_scale
ys_trva = (y_trva - y_center) / y_scale
Xs_te = (X_te - x_center) / x_scale

print(f"forecast origins: {n}, train+validation: {split_valid}, test: {n - split_valid}")

forecast origins: 2494, train+validation: 1995, test: 499

Show the code

print(f"test window: {dates_te.iloc[0].date()} to {dates_te.iloc[-1].date()}")

test window: 2023-01-06 to 2024-12-31

The make_log_har helper builds the same three HAR averages used in the Chapter 8 illustration but applies \log(\cdot) before standardization, so the linear baseline is a log-log HAR regression and the FNN inputs are on the same scale. We standardize features and the target on the train-plus-validation block, which avoids leaking test-block information into the scaling constants.

The Gaussian baseline is just OLS plus an in-sample residual variance.

Show the code

ols = LinearRegression().fit(X_trva, y_trva)
mu_ols_te = ols.predict(X_te)
resid_trva = y_trva - ols.predict(X_trva)
sigma_ols = float(np.std(resid_trva))  # ML convention: divide the SSE by n
sigma_ols_te = np.full_like(mu_ols_te, sigma_ols)

The distributional HAR-FNN has two heads on top of a shared hidden layer. The custom loss implements the per-observation Gaussian NLL of Section 9.5, with a smooth positive variance parameterization that prevents a degenerate density when residuals are near zero (the pathology analyzed in Exercise 9.1). Exercise 9.1 derives the optimum under a hard inequality constraint; the implementation uses a smooth map because applying a hard maximum to a network output would give zero derivative whenever the floor binds.

Two protocol notes distinguish this workflow from the previous chapter’s. First, we fix all hyperparameters a priori—32 hidden units, learning rate 10^{-3}, 200 epochs, and batch size 32—without a data-driven tuning step. We therefore use the first 80% of forecast origins as one estimation sample; the variable names X_trva and Xs_trva record that this block combines the earlier training and validation samples. Second, standardizing on this full block is appropriate because no validation sample needs to remain protected, while the test block still never enters the scaling constants. The Chapter 8 warning about scaling beyond the training block applies when validation data guide model selection. The results below describe one fixed configuration and a single training run with a fixed random seed, not a stable architecture ranking; repeated seeds or resamples would be needed to quantify optimization variability.

Show the code

VAR_FLOOR = 1e-3

def make_dist_har_fnn(n_neurons=32, learning_rate=1e-3):
    inputs = keras.Input(shape=(3,))
    h = layers.Dense(n_neurons, activation="softplus")(inputs)
    mu_head = layers.Dense(1)(h)
    raw_var_head = layers.Dense(1)(h)
    outputs = layers.Concatenate()([mu_head, raw_var_head])
    model = keras.Model(inputs, outputs, name="DistHARFNN")

    def gaussian_nll(y_true, y_pred):
        mu = y_pred[..., 0:1]
        raw_var = y_pred[..., 1:2]
        var = VAR_FLOOR + ops.softplus(raw_var)
        return 0.5 * (ops.log(2.0 * np.pi * var) + ops.square(y_true - mu) / var)

    model.compile(loss=gaussian_nll, optimizer=keras.optimizers.Adam(learning_rate))
    return model

dist_model = make_dist_har_fnn(n_neurons=32, learning_rate=1e-3)
dist_model.fit(
    Xs_trva, ys_trva.reshape(-1, 1),
    epochs=200, batch_size=32, verbose=0, shuffle=False,
)

<keras.src.callbacks.history.History object at 0x3231e4440>

Show the code


pred_std = dist_model.predict(Xs_te, verbose=0)
mu_std, raw_var_std = pred_std[:, 0], pred_std[:, 1]
mu_dist_te = mu_std * y_scale + y_center
var_dist_te = (VAR_FLOOR + np.logaddexp(0.0, raw_var_std)) * (y_scale ** 2)
sigma_dist_te = np.sqrt(var_dist_te)

Because the network is trained on the standardized target \tilde y_{t+1}=(y_{t+1}-\bar y)/s_y paired with \mathbf{x}_t, the map c+\operatorname{softplus}(a_\phi(\mathbf{x}_t)) represents the conditional variance of the standardized outcome \tilde Y_{t+1}:=(Y_{t+1}-\bar y)/s_y given \mathbf{x}_t, with \bar y and s_y treated as fixed training constants. Linear standardization is variance-equivariant:

\operatorname{Var}(Y_{t+1}\mid \mathbf{x}_t) =s_y^2\operatorname{Var}(\tilde Y_{t+1}\mid \mathbf{x}_t).

Thus, the back-transformed predictive variance is [c+\operatorname{softplus}(a_\phi(\mathbf{x}_t))]s_y^2, exactly as implemented. This back-transformation reverses only the linear standardization and leaves the forecast on the log-RV scale. Obtaining a predictive distribution for RV instead requires the nonlinear transformation RV=\exp(\log RV): a Gaussian predictive density for \log RV implies a lognormal density for RV, whose level-scale mean is \mathbb{E}[RV]=\exp(\mu+\sigma^2/2). The present application evaluates forecasts and outcomes on the log-RV scale and does not perform this level-scale transformation. Under the scores considered here, consistent linear rescaling preserves model rankings (Kleen 2024), although absolute numerical score levels generally change. Forecasts and outcomes must therefore be placed on one common scale before their scores are compared.

Show the code

logs_ols = -norm.logpdf(y_te, mu_ols_te, sigma_ols_te)
logs_dist = -norm.logpdf(y_te, mu_dist_te, sigma_dist_te)
crps_ols = np.array([ps.crps_gaussian(y, mu=m, sig=s)
                     for y, m, s in zip(y_te, mu_ols_te, sigma_ols_te)])
crps_dist = np.array([ps.crps_gaussian(y, mu=m, sig=s)
                      for y, m, s in zip(y_te, mu_dist_te, sigma_dist_te)])
pit_ols = norm.cdf(y_te, mu_ols_te, sigma_ols_te)
pit_dist = norm.cdf(y_te, mu_dist_te, sigma_dist_te)

scores = pd.DataFrame({
    "model": ["HAR-OLS Gaussian", "Distributional HAR-FNN"],
    "mean_LogS": [logs_ols.mean(), logs_dist.mean()],
    "mean_CRPS": [crps_ols.mean(), crps_dist.mean()],
})
scores.round(4)

                    model  mean_LogS  mean_CRPS
0        HAR-OLS Gaussian     0.8767     0.3218
1  Distributional HAR-FNN     0.8966     0.3246

Figure 9.2 places the two test-block forecasts on a common scale. Compare the constant-width baseline band in the top panel with the feature-dependent band below it.

Show the code

fig, axes = plt.subplots(2, 1, figsize=(11, 6.4), sharex=True)

for ax, mu, sig, label, color in zip(
    axes,
    [mu_ols_te, mu_dist_te],
    [sigma_ols_te, sigma_dist_te],
    ["HAR-OLS Gaussian", "Distributional HAR-FNN"],
    ["C1", "C0"],
):
    ax.plot(dates_te, y_te, color="0.55", linewidth=0.9, label="realized log RV")
    ax.plot(dates_te, mu, color=color, linewidth=1.1, label=f"{label} mean")
    ax.fill_between(dates_te, mu - 1.96 * sig, mu + 1.96 * sig,
                    color=color, alpha=0.20, label="95% band")
    ax.set_ylabel("log RV")
    ax.set_title(label)
    ax.legend(frameon=False, loc="upper right", ncol=3, fontsize=9)
    ax.grid(True, alpha=0.25)

axes[1].set_xlabel("Date")
plt.tight_layout()
plt.show()

Figure 9.2: Realized \log RV_{t+1} on the test block (gray) with one-step-ahead Gaussian density forecasts. In each panel, the colored line is the predictive mean and the shaded band is the central 95% predictive interval; the top panel shows the HAR-OLS Gaussian baseline, the bottom panel the distributional HAR-FNN.

The two predictive bands give the first part of the story directly. The HAR-OLS Gaussian band is a parallel pair of lines: a single training-block residual variance is replayed at every test-block date, so the forecast distribution is just a translated copy of the same bell curve. The distributional HAR-FNN band tightens in quiet stretches where the HAR averages are small and widens in stress periods—the predictive standard deviation varies by roughly a factor of 1.5 across the test block, while the baseline holds it constant. The variance head therefore produces a feature-dependent scale; whether that variation improves the density forecast is precisely what the score and calibration diagnostics below must judge.

Show the code

fig, (ax_logs, ax_crps, ax_pit_ols, ax_pit_dist) = plt.subplots(1, 4, figsize=(15, 4))

labels = ["HAR-OLS\nGaussian", "Distributional\nHAR-FNN"]
xpos = np.arange(len(labels))
for ax, vals, title in [
    (ax_logs, [logs_ols.mean(), logs_dist.mean()], "Average LogS"),
    (ax_crps, [crps_ols.mean(), crps_dist.mean()], "Average CRPS"),
]:
    ax.bar(xpos, vals, 0.55, color=["C1", "C0"], alpha=0.85)
    ax.set_xticks(xpos)
    ax.set_xticklabels(labels)
    ax.set_ylabel("Average score (lower = better)")
    ax.set_title(title)
    ax.grid(True, axis="y", alpha=0.3)

bins = np.linspace(0, 1, 16)
for ax, pit_vals, title, color in zip(
    [ax_pit_ols, ax_pit_dist],
    [pit_ols, pit_dist],
    ["PIT: HAR-OLS Gaussian", "PIT: Distributional HAR-FNN"],
    ["C1", "C0"],
):
    ax.hist(pit_vals, bins=bins, density=True, color=color,
            edgecolor="white", alpha=0.85)
    ax.axhline(1.0, color="red", linestyle="--", linewidth=1, label="Uniform")
    ax.set_xlabel("PIT value")
    ax.set_ylabel("Density")
    ax.set_title(title)
    ax.set_ylim(0, 3.0)
    ax.legend(frameon=False)
    ax.grid(True, alpha=0.25)

plt.tight_layout()
plt.show()

Figure 9.3: Average LogS (first panel) and average CRPS (second panel) on the test block—each score in its own panel, because the two are measured in different units and are comparable only across models within a panel—and PIT histograms for the two forecasts (third and fourth panels). A well-calibrated forecast yields a flat PIT histogram close to the uniform density (red dashed line).

The averages and PITs in Figure 9.3 give the second part: the time-varying variance does not improve either the scores or the PIT diagnostic. The HAR-OLS Gaussian baseline has the lower average LogS (0.8767 against 0.8966) and the lower average CRPS (0.3218 against 0.3246). As in the previous chapter, we do not test these differences formally; a Diebold–Mariano-type comparison of the score differentials with heteroskedasticity- and autocorrelation-consistent (HAC) standard errors would be the appropriate formal comparison. The PIT histograms are centrally concentrated rather than flat: both forecasts produce too few observations in the outer bins and too many near the middle, descriptive evidence consistent with predictive distributions that are too wide. The central concentration is more pronounced for the distributional HAR-FNN. Because the PIT sample is finite and the PIT sequence may be serially dependent, the histograms are diagnostics rather than formal calibration tests.

This comparison ranks two complete forecasting procedures; it does not isolate the effect of adding a variance head. The procedures differ in their mean specification, variance specification, and optimization problem. Their score and PIT differences could therefore reflect the nonlinear mean, feature-dependent variance, joint optimization, or interactions among them. A sharper diagnostic comparison would hold the fitted mean fixed and vary only the variance estimator—for example, a constant or two-step variance against a jointly trained variance head. Applying the full HNN stabilization procedure would be a natural next experiment, but the present comparison cannot attribute its result to any single ingredient.

The methodological lesson is that no single diagnostic suffices. A practitioner who reads only the predictive bands might conclude that the distributional model is better because its uncertainty changes over time. The score comparison and PIT histograms show why that conclusion would be premature: in this test block, the extra variation in the predicted scale coincides with worse average scores and more pronounced overdispersion. Together, Figure 9.2 and Figure 9.3 distinguish what the variance head changes from whether that change improves the density forecast.

9.8 Summary

Key Takeaways

Distributional networks predict parameters of a conditional distribution rather than only a point forecast.
Negative log-likelihood minimizes LogS; CRPS provides an alternative when it is differentiable for the chosen model.
Hemisphere networks stabilize joint learning through a constraint emphasizing volatility, blocked subsampling and blocked out-of-bag recalibration, in addition to separate heads.
Mixture density networks combine component distributions to represent shapes such as skewness and multiple modes.

Common Pitfalls

Producing a valid density does not establish calibration, which must be checked on held-out outcomes.
Enforce positive scale outputs without overlooking zero-gradient regions created by hard floors.
Mixture components are not identified economic regimes without additional restrictions.
Validate added distributional flexibility against simpler forecasts when sample size is limited.

9.9 Exercises

Exercise 9.1: Why Joint Gaussian NLL Training Can Be Ill-Posed

Consider a Gaussian distributional model for one-step-ahead forecasting. Maximizing the likelihood is equivalent to maximizing the log-likelihood, which in turn is equivalent to minimizing the negative log-likelihood. We therefore work with the per-observation Gaussian negative log-likelihood, up to an additive constant,

\ell_{t+1}(e_{t+1},\sigma_t^2)=\frac{1}{2}\log(\sigma_t^2)+\frac{e_{t+1}^2}{2\sigma_t^2}, \qquad e_{t+1}=y_{t+1}-\mu_t,\qquad \sigma_t^2>0.

For fixed e_{t+1} \neq 0, derive the value of \sigma_t^2 that minimizes \ell_{t+1}(e_{t+1},\sigma_t^2). Verify that this critical point is a minimum.
Show that if e_{t+1}=0, then the infimum of \ell_{t+1}(e_{t+1},\sigma_t^2) over \sigma_t^2>0 equals -\infty and is not attained at any finite variance value. Conclude that if an overparameterized mean network interpolates the training sample exactly, so that e_{t+1}=0 for all t=1,\dots,T, and the variance head can choose each \sigma_t^2 freely, the sample Gaussian NLL has no finite minimizer.
Suppose instead that the model imposes a variance floor \sigma_t^2 \ge c for some constant c>0. Solve the constrained minimization problem for one observation and show that the optimal variance is \sigma_t^{2\ast}=\max\{e_{t+1}^2,c\}. Explain briefly why this removes the pathology from Part 2, but does not by itself guarantee good out-of-sample calibration.
Assume Y_{t+1} and \mathbb{E}[Y_{t+1}\mid\mathbf{X}_t] have finite second moments and that the process is stationary and ergodic. Use the law of total variance to show that \mathbb{E}[\operatorname{Var}(Y_{t+1}\mid\mathbf{X}_t)]\leq\operatorname{Var}(Y_{t+1}), with strict inequality when the conditional mean is nonconstant.
Suppose a flexible mean interpolates every training outcome. State its in-sample average squared residual, explain why those residuals cannot calibrate the variance scale, and explain why blocked out-of-bag forecast errors are preferable when overlapping time-series windows share observations.

Exam level. Parts 2–5 test whether students understand the instability behind joint mean-variance training and the motivation for out-of-bag calibration rather than only the Gaussian likelihood formula.

Hint for Part 3

First solve the unconstrained problem from Part 1. Then check whether that solution satisfies the inequality constraint \sigma_t^2 \ge c.

Solution

Part 1

Differentiate with respect to \sigma_t^2:

\frac{\partial \ell_{t+1}}{\partial \sigma_t^2} = \frac{1}{2\sigma_t^2}-\frac{e_{t+1}^2}{2(\sigma_t^2)^2}.

Setting this equal to zero gives

\frac{1}{2\sigma_t^2}=\frac{e_{t+1}^2}{2(\sigma_t^2)^2} \quad\Longrightarrow\quad \sigma_t^2=e_{t+1}^2.

To verify that this is a minimum, compute the second derivative:

\frac{\partial^2 \ell_{t+1}}{\partial (\sigma_t^2)^2} = -\frac{1}{2(\sigma_t^2)^2}+\frac{e_{t+1}^2}{(\sigma_t^2)^3}.

Evaluating at \sigma_t^2=e_{t+1}^2 yields

\frac{\partial^2 \ell_{t+1}}{\partial (\sigma_t^2)^2}\Big|_{\sigma_t^2=e_{t+1}^2} = \frac{1}{2e_{t+1}^4}>0.

Hence the minimizer for fixed e_{t+1}\neq 0 is \sigma_t^2=e_{t+1}^2.

Part 2

If e_{t+1}=0, then

\ell_{t+1}(0,\sigma_t^2)=\frac{1}{2}\log(\sigma_t^2).

As \sigma_t^2 \downarrow 0,

\frac{1}{2}\log(\sigma_t^2)\to -\infty.

Therefore the infimum of the objective is -\infty. However, no finite positive value of \sigma_t^2 attains this infimum, so the problem has no minimizer.

For the sample statement, if the mean network interpolates the training sample exactly, then e_{t+1}=0 for all t. The sample objective is

\sum_{t=1}^{T}\ell_{t+1}(0,\sigma_t^2)=\frac{1}{2}\sum_{t=1}^{T}\log(\sigma_t^2).

By sending each \sigma_t^2 arbitrarily close to zero, this sum can be made arbitrarily negative. Hence the training Gaussian NLL has no finite minimizer. This is the basic ill-posedness behind unstable joint training of very flexible mean and variance networks.

Part 3

From Part 1, the unconstrained minimizer is \sigma_t^2=e_{t+1}^2. Under the constraint \sigma_t^2\ge c, there are two cases:

If e_{t+1}^2\ge c, the unconstrained optimum is feasible, so \sigma_t^{2\ast}=e_{t+1}^2.
If e_{t+1}^2<c, the unconstrained optimum violates the constraint, so the constrained optimum lies at the boundary, \sigma_t^{2\ast}=c.

Therefore

\sigma_t^{2\ast}=\max\{e_{t+1}^2,c\}.

The positive variance floor removes the pathology from Part 2 because even if e_{t+1}=0, the smallest admissible variance is now c, so the objective cannot be driven to -\infty by shrinking the variance to zero. But this does not guarantee good out-of-sample calibration: a variance floor prevents a mathematical degeneracy, yet the learned variance process can still be badly misspecified for future data.

Part 4

The law of total variance gives

\operatorname{Var}(Y_{t+1}) =\mathbb{E}\!\left[\operatorname{Var}(Y_{t+1}\mid\mathbf{X}_t)\right] +\operatorname{Var}\!\left(\mathbb{E}[Y_{t+1}\mid\mathbf{X}_t]\right).

The second term is nonnegative, which proves the weak inequality. It is positive when the conditional mean is nonconstant, giving strict inequality. Stationarity and ergodicity justify estimating these population moments with time averages.

Part 5

If the fitted mean interpolates every training outcome, then e_{t+1}=0 for every training origin and

\frac{1}{T}\sum_{t=1}^{T}e_{t+1}^2=0.

These in-sample residuals contain no information about the variance of genuine forecast errors and therefore cannot calibrate \kappa or the overall variance scale. An out-of-bag error at time t+1 comes from a fit that did not train on that target. Holding out contiguous blocks also reduces the direct information sharing created when adjacent, overlapping windows contain many of the same observations. Blocked out-of-bag errors are therefore more credible proxies for time-series forecast errors than either interpolating residuals or errors from randomly held-out individual dates.

Exercise 9.2: Moment-Matched Gaussian Forecasts versus MDNs

Consider the symmetric two-component conditional density

f_M(y\mid \mathbf{x}_t) = \frac{1}{2}\,\varphi(y;-m_t,s^2)+\frac{1}{2}\,\varphi(y;m_t,s^2), \qquad m_t>0,\ s>0,

where \varphi(y;\mu,\sigma^2) denotes the Gaussian density with mean \mu and variance \sigma^2.

Now consider the Gaussian forecast

f_G(y\mid \mathbf{x}_t)=\varphi(y;0,m_t^2+s^2),

which matches the first two moments of f_M.

Derive the conditional mean and conditional variance of f_M(y\mid \mathbf{x}_t) and verify that f_G indeed matches these two moments.
Show that the point forecast that minimizes conditional MSE under both f_M and f_G is zero. Explain why an MSE comparison cannot distinguish these two forecast distributions.
Evaluate both densities at the realized outcome y_{t+1}=m_t: f_M(m_t\mid \mathbf{x}_t) \qquad\text{and}\qquad f_G(m_t\mid \mathbf{x}_t). Show that as m_t \to \infty with s fixed, the mixture density at y_{t+1}=m_t stays bounded away from zero, while the Gaussian density goes to zero. Conclude that LogS can strongly prefer the MDN-style forecast even when MSE cannot distinguish the two models.
Show formally that swapping the labels of the two mixture components leaves f_M(y\mid \mathbf{x}_t) unchanged. Explain why this means that mixture components are not structurally identified without additional restrictions.

Exam level. The point of this exercise is that matching means and variances is not enough to match the predictive distribution, and that distributional scores can detect this while MSE cannot.

Hint for Part 1

Use the symmetry of the mixture. For the variance, compute \mathbb{E}[Y_{t+1}^2\mid \mathbf{x}_t] first.

Solution

Part 1

By symmetry,

\mathbb{E}[Y_{t+1}\mid \mathbf{x}_t] = \frac{1}{2}(-m_t)+\frac{1}{2}(m_t)=0.

For the second moment,

\mathbb{E}[Y_{t+1}^2\mid \mathbf{x}_t] = \frac{1}{2}(m_t^2+s^2)+\frac{1}{2}(m_t^2+s^2) = m_t^2+s^2.

Hence

\operatorname{Var}(Y_{t+1}\mid \mathbf{x}_t) = \mathbb{E}[Y_{t+1}^2\mid \mathbf{x}_t] - \big(\mathbb{E}[Y_{t+1}\mid \mathbf{x}_t]\big)^2 = m_t^2+s^2.

The Gaussian forecast f_G(y\mid \mathbf{x}_t)=\varphi(y;0,m_t^2+s^2) therefore matches the conditional mean and variance of the mixture forecast.

Part 2

Under squared-error loss, the optimal point forecast is the conditional mean. Since both forecast distributions have conditional mean zero, the MSE-optimal point forecast is

\hat{y}_{t+1\mid t}=0.

Therefore an MSE comparison cannot distinguish the two forecast distributions: both imply exactly the same optimal point forecast, even though one forecast is bimodal and the other is unimodal.

Part 3

For the mixture density at y_{t+1}=m_t,

f_M(m_t\mid \mathbf{x}_t) = \frac{1}{2}\varphi(m_t;-m_t,s^2)+\frac{1}{2}\varphi(m_t;m_t,s^2).

Using the Gaussian density formula,

f_M(m_t\mid \mathbf{x}_t) = \frac{1}{2}\frac{1}{\sqrt{2\pi}s}e^{-2m_t^2/s^2} + \frac{1}{2}\frac{1}{\sqrt{2\pi}s} = \frac{1}{2\sqrt{2\pi}s}\Big(1+e^{-2m_t^2/s^2}\Big).

For the Gaussian forecast,

f_G(m_t\mid \mathbf{x}_t) = \frac{1}{\sqrt{2\pi(m_t^2+s^2)}}\exp\!\left(-\frac{m_t^2}{2(m_t^2+s^2)}\right).

As m_t \to \infty with s fixed, we have

f_M(m_t\mid \mathbf{x}_t)\to \frac{1}{2\sqrt{2\pi}s},

which is strictly positive, whereas

f_G(m_t\mid \mathbf{x}_t) \sim \frac{e^{-1/2}}{\sqrt{2\pi}\,m_t}\to 0.

So when the realization lands near one of the two modes and the modes are well separated, the mixture forecast assigns much higher density to the realized value. Therefore the mixture forecast achieves a substantially lower (better) LogS, even though MSE cannot distinguish it from the moment-matched Gaussian.

Part 4

If we swap the component labels, the predictive density becomes

\frac{1}{2}\varphi(y;m_t,s^2)+\frac{1}{2}\varphi(y;-m_t,s^2),

which is exactly the same expression as the original density, just with the two terms reversed. Hence the predictive density is unchanged by relabeling the components.

This is the label-switching problem. It implies that the two mixture components are not structurally identified from predictive fit alone. To interpret them as economic regimes, one would need additional restrictions, theory, or identifying assumptions beyond the MDN likelihood.

9.10 References

Bishop, Christopher M. 1994. “Mixture Density Networks.” Technical Report NCRG/94/004. Birmingham, U.K.: Neural Computing Research Group, Aston University.

Bollerslev, Tim, and Jeffrey M. Wooldridge. 1992. “Quasi-Maximum Likelihood Estimation and Inference in Dynamic Models with Time-Varying Covariances.” Econometric Reviews 11 (2): 143–72. https://doi.org/10.1080/07474939208800229.

Gneiting, Tilmann, and Matthias Katzfuss. 2014. “Probabilistic Forecasting.” Annual Review of Statistics and Its Application 1: 125–51. https://doi.org/10.1146/annurev-statistics-062713-085831.

Goulet Coulombe, Philippe, Mikael Frenette, and Karin Klieber. 2026. “From Reactive to Proactive Volatility Modeling With Hemisphere Neural Networks.” Journal of Applied Econometrics 41 (3): 265–79. https://doi.org/10.1002/jae.70042.

Kleen, Onno. 2024. “Scaling and Measurement Error Sensitivity of Scoring Rules for Distribution Forecasts.” Journal of Applied Econometrics 39 (5): 833–49. https://doi.org/10.1002/jae.3056.

Footnotes

The data file data/taq_spy/SPY_daily_measures.csv covers 2015-01-02 through 2024-12-31 and is documented in the dataset appendix. The longest HAR component averages 21 trading days. The pipeline nevertheless begins forecasting only after 22 lagged observations so that its forecast origins and final 20% test block match those in the previous chapter.↩︎