8 Empirical Illustration: Neural Networks for Realized Volatility

8.1 Overview

Neural-network architectures are often introduced as if the main empirical question were which architecture is most powerful. In time-series econometrics, that is usually the wrong starting point. Before comparing a feed-forward network, a recurrent network, or a Long Short-Term Memory (LSTM) network, we must define the forecast target, the information set at each forecast origin, the validation scheme, and the loss function used for evaluation. Throughout this chapter, volatility names the broad forecasting application; the measured target is realized variance, not its square root.

This chapter makes those choices explicit in a single empirical illustration. The object is one-step-ahead forecasting of the realized variance of the S&P 500 exchange-traded fund (ETF; ticker SPY), using the daily measures described in the dataset appendix.

The empirical comparison has two purposes. First, it compares a log-target pipeline with a level-target pipeline. The comparison changes both the target transformation and the normalization of the lag inputs, so it does not isolate the effect of the target transformation alone. Second, it compares a linear benchmark, the heterogeneous autoregressive (HAR) model estimated by ordinary least squares (OLS), with four neural specifications: a HAR feed-forward neural network (HAR-FNN) on the same three volatility summaries, a direct FNN on 22 individual daily observations, a Simple RNN (recurrent neural network), and an LSTM, with both recurrent models using 756 daily observations. The long input history lets us ask whether learning from roughly three years of daily observations improves on a parsimonious volatility summary. In this illustration, HAR remains competitive with the neural models, and the ranking depends on the evaluation criterion.

8.2 Roadmap

1. We define the realized-variance target and the information available at each forecast origin.
2. We separate the 756-day input history from the sliding estimation sample and specify annual refitting.
3. We compare HAR summaries, a dense lag vector, and recurrent sequence representations.
4. We select neural hyperparameters before the test period and apply the fixed refitting rule.
5. We compare forecast errors across models and loss functions and assess what the rankings can establish.

8.3 Data and Forecast Target

Let r_{t,j} denote the jth intraday return on trading day t. A standard realized-variance measure is

RV_t=\sum_{j=1}^{M} r_{t,j}^2,

where M is the number of intraday return intervals. The empirical series uses 5-minute log returns multiplied by 100, so RV and its root mean squared error (RMSE) on the level scale are measured in squared percentage points. Realized variance RV_t is a canonical high-frequency volatility measure; Andersen et al. (2003) is a standard reference for modeling and forecasting it.¹ It is positive, highly right-skewed, and persistent. For that reason, a common forecasting target is

z_t = \log(RV_t),

although the empirical comparison below also reports what happens when the networks are fit directly to RV_t in levels.

A note on case: this chapter is a worked empirical illustration, so its formulas operate on the observed series: lowercase z_t and r_{t,j} denote realized values, in line with the book’s convention. In the few analytic statements that treat the target as random, we write Z_t = \log RV_t for the random variable whose realization is z_t; the capitalization of the acronym RV_t does not itself mark randomness, and where RV_t appears inside a conditional expectation it denotes the random realized variance.

For the recurrent models with a log target, the one-step-ahead forecasting problem is

\widehat z_{t+1\mid t} = f_\theta(z_{t-755},\ldots,z_t),

where the input contains 756 daily observations, approximately three trading years at 252 observations per year. The sequence is ordered from oldest to newest. We forecast after observing day t; the target z_{t+1} is excluded from the input. The Simple RNN and LSTM process these observations sequentially. The direct FNN instead receives only the most recent 22 daily observations, (z_{t-21},\ldots,z_t), as a fixed-position vector. For the level-target versions, replace the log observations by the corresponding RV observations.

The HAR specifications instead receive three summaries known at origin t:

\left( RV_t,\quad \frac{1}{5}\sum_{j=0}^{4}RV_{t-j},\quad \frac{1}{21}\sum_{j=0}^{20}RV_{t-j} \right).

The daily, weekly, and monthly components are the parsimonious volatility representation of Corsi (2009). For the log target, we take the logarithm of each component; for the level target, we use each component in levels. In particular, the weekly and monthly log inputs are logs of average RV, not averages of log RV.

Input history and estimation sample. The input history describes one forecast example: 22 daily values for the direct FNN and 756 for the recurrent models. Separately, we fit each model on the most recent 756 training examples, each pairing a past input with its subsequently observed target. At an annual refit, older examples leave this sliding estimation sample as newer examples enter. The sample never expands. Each training example still needs its own lag history, so preparing 756 recurrent-model examples with 756 daily inputs requires 1,512 consecutive daily observations, roughly six trading years. The earlier observations supply predictors; they do not add training targets.

The data contain 2,516 trading-day observations from 2015-01-02 to 2024-12-31. After constructing the long recurrent input windows, 1,760 forecast examples remain; the other models use the same target dates for comparability. We use complete calendar years: 2022 for validation (251 trading-day predictions) and 2023–2024 for testing (502 predictions). Each annual fit uses the most recent 756 observed training outcomes available before that year starts. The validation year therefore follows the same fit-once, forecast-for-one-year rule as each test year.

Figure 8.1 distinguishes the information used for one forecast from the examples used to estimate its model. Within each validation or test year, the weights and scaling constants stay fixed, while the daily input window advances after every observed outcome.

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib.patches import Rectangle
from datetime import datetime

fig, axes = plt.subplots(2, 1, figsize=(11, 6.2), height_ratios=[1, 1.7], layout="constrained")
ax = axes[0]
ax.set(xlim=(0, 12), ylim=(0, 3.5))
for x, label in [(0.5, "$z_{t-755}$"), (2.3, "$z_{t-754}$"), (6, "$z_{t-1}$"), (7.8, "$z_t$")]:
    ax.add_patch(Rectangle((x, 1.2), 1.35, .8, facecolor="#d7e5f5", edgecolor="#4c78a8"))
    ax.text(x+.675, 1.6, label, ha="center", va="center", fontsize=12)
ax.text(4.8, 1.6, "...", ha="center", fontsize=18)
ax.add_patch(Rectangle((10.2, 1.2), 1.35, .8, facecolor="#f7d7c4", edgecolor="#d65f00"))
ax.text(10.875, 1.6, "$z_{t+1}$", ha="center", va="center", fontsize=12)
ax.annotate("", xy=(10.1, 1.6), xytext=(9.3, 1.6), arrowprops={"arrowstyle":"->"})
ax.text(4.8, 2.55, "RNN / LSTM: 756 inputs known at origin t", ha="center", fontsize=12)
ax.text(10.875, 2.55, "Target", ha="center", fontsize=12)
ax.axis("off")

ax = axes[1]
def day(s):
    return mdates.date2num(datetime.fromisoformat(s))
rows = [
    (2, "Validation", "2019-01-03", "2021-12-31", 2022, "#fff0bf"),
    (1, "Test", "2020-01-02", "2022-12-30", 2023, "#f7d7c4"),
    (0, "Test refit", "2020-12-29", "2023-12-29", 2024, "#f7d7c4"),
]
for y, label, start, end, year, color in rows:
    ax.barh(y, day(end)-day(start), left=day(start), height=.56, color="#cfe8cf", edgecolor="0.45")
    ax.text((day(start)+day(end))/2, y, "756 training outcomes", ha="center", va="center", fontsize=10)
    prediction_start, prediction_end = day(f"{year}-01-01"), day(f"{year+1}-01-01")
    ax.barh(y, prediction_end-prediction_start, left=prediction_start, height=.56, color=color, edgecolor="0.45")
    ax.text((prediction_start+prediction_end)/2, y, str(year), ha="center", va="center", fontsize=11)
ax.set_yticks([r[0] for r in rows], [r[1] for r in rows])
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
ax.set_xlim(day("2018-10-01"), day("2025-03-01"))
ax.set_xlabel("Green: training outcomes; yellow/orange: prediction year. Lag inputs start earlier.")
for spine in ["top", "right", "left"]:
    ax.spines[spine].set_visible(False)
plt.show()

Figure 8.1: Two different windows in the forecasting design. Top: a recurrent-model forecast uses the preceding 756 daily observations, ordered from oldest to newest. Bottom: each row is one annual fit, with its 756 training outcomes in green and its prediction year in yellow (validation) or orange (test). The training window slides forward each year. Validation predictions cover 2022; only predictions for 2023 and 2024 enter the reported test scores.

The table separates validation, whose prediction errors select the neural configuration, from testing, whose prediction errors assess the selected procedure. Every row uses a sliding estimation sample of 756 training outcomes. The middle column is the information cutoff for fitting; the final column is the calendar year being predicted, one trading day at a time.

Phase	Fit using outcomes observed through	Predict each trading day in
Validation: choose the configuration	December 2021	2022
Test: first evaluation year	December 2022	2023
Test: annual refit	December 2023	2024

The validation predictions run from 3 January through 30 December 2022. Hyperparameters are then fixed before the first test prediction on 3 January 2023. At the beginning of 2024, we refit the weights and scalers using the same chosen configuration. Only the 2023 and 2024 prediction errors enter the reported test scores.

A 756-example training sample is approximately three trading years; its start need not coincide with 1 January. The exact training ranges are 2019-01-03–2021-12-31 for validation, 2020-01-02–2022-12-30 for the first test year, and 2020-12-29–2023-12-29 for the second test year. The 2022 validation outcomes can enter the 2023 fit once observed, and the 2023 test outcomes can enter the 2024 fit. Neither update changes an earlier forecast or uses the test errors to choose hyperparameters.

Question for Reflection

The 2024 model is estimated using realized variances observed during 2023, even though 2023 belongs to the test period. Why does this not leak future information into the reported forecasts?

Suggested Answer

The annual refitting rule is fixed before testing begins. Each 2023 forecast uses only information available at its forecast origin, and the 2024 fit uses only outcomes already known at its forecast origin. Reusing an observed outcome for a later fit is compatible with sequential evaluation; using it to revise an earlier forecast or choose hyperparameters after inspecting test performance would change the protocol.

The following code constructs the examples and selects the sliding estimation sample. The code indexes rows by the target date, so a row with target at array position t uses positions t-756 through t-1 for the recurrent models and the last 22 of those positions for the direct FNN. This is the same timing as the formula above, with the target index shifted by one. The complete implementation is available in the replication module.

from pathlib import Path
import time
import numpy as np
import pandas as pd
from nn_example_pipeline import (
    INPUT_DAYS, FNN_INPUT_DAYS, ESTIMATION_EXAMPLES, VALIDATION_START, TEST_START,
    TEST_END, N_TRIALS, MAX_EPOCHS, SEED,
)

def load_data(path):
    data = pd.read_csv(path, parse_dates=["date"]).sort_values("date")
    data = data.loc[data.date <= TEST_END, ["date", "rv"]].reset_index(drop=True)
    if data.date.duplicated().any() or not data.date.is_monotonic_increasing:
        raise ValueError("Dates must be distinct and increasing.")
    if not np.isfinite(data.rv).all() or not (data.rv > 0).all():
        raise ValueError("Every observed trading-day RV must be finite and positive.")
    return data

def make_inputs(data, input_days=INPUT_DAYS):
    """One row per target date; all inputs end on its preceding trading day."""
    rv = data.rv.to_numpy(dtype=float)
    z = np.log(rv)
    targets = np.arange(input_days, len(data))
    lag_rv = np.stack([rv[t-input_days:t] for t in targets])
    har_rv = np.array([[rv[t-1], rv[t-5:t].mean(), rv[t-21:t].mean()]
                       for t in targets])
    return {
        "dates": data.date.iloc[targets].reset_index(drop=True),
        "target_positions": targets,
        "lag_log": np.log(lag_rv), "lag_rv": lag_rv,
        "har_log": np.log(har_rv), "har_rv": har_rv,
        "y_log": z[targets], "y_rv": rv[targets],
    }

def estimation_indices(dates, first_forecast_date, n_examples=ESTIMATION_EXAMPLES):
    """The last n_examples outcomes observed before this forecast target."""
    idx = np.flatnonzero(np.asarray(dates < first_forecast_date))
    if len(idx) < n_examples:
        raise ValueError("Not enough forecast-origin examples for the estimation window.")
    return idx[-n_examples:]

data = load_data(Path("data/taq_spy/SPY_daily_measures.csv"))
panel = make_inputs(data)
train_idx = estimation_indices(panel["dates"], VALIDATION_START)

8.4 Comparable Specifications

With the forecast dates fixed, the specifications differ in how they represent the past.

• HAR-OLS regresses the target on the three HAR components, in logs or levels according to the target.
• HAR-FNN uses exactly the same three components through one hidden dense layer and a linear output.
• The direct FNN jointly processes the most recent 22 individual daily inputs through one hidden dense layer and a linear output; each lag position receives its own weights.
• The Simple RNN processes the 756 inputs from oldest to newest through one recurrent layer and a linear output.
• The LSTM processes the same sequence through one recurrent layer with gated memory and a linear output.

The HAR-OLS and HAR-FNN comparison holds the target and predictor information fixed. It asks whether a tuned nonlinear procedure improves on linear OLS with the same summaries. The direct FNN receives 22 individual observations instead of the HAR summaries, while the recurrent models receive 756. Their comparison with HAR therefore changes both the representation and the available information. Comparing the direct FNN with the recurrent models also changes the history length, so the results do not isolate the effect of recurrence. The direct FNN is a neural network on individual lags, not a HAR model without averaging.

Figure 8.2 shows where aggregation and parameter sharing enter. Increasing the input history raises the dense FNN’s parameter count but does not raise the recurrent layers’ parameter counts at a fixed width. It does increase the number of recurrent updates and the length of the gradient path.

import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyArrowPatch

fig, axes = plt.subplots(2, 2, figsize=(11.5, 6.4))
axes = axes.ravel()
titles = ["HAR-FNN", "FNN", "Simple RNN", "LSTM"]
for ax, title in zip(axes, titles):
    ax.set_xlim(0, 10)
    ax.set_ylim(0, 6)
    ax.set_title(title, fontsize=13)
    ax.axis("off")

for i, label in enumerate(["$RV_t$", "avg 5", "avg 21"]):
    axes[0].add_patch(Rectangle((0.7, 4.0 - i * 1.0), 1.3, 0.55, facecolor="#d7e5f5", edgecolor="#4c78a8"))
    axes[0].text(1.35, 4.27 - i * 1.0, label, ha="center", va="center", fontsize=9)
axes[0].text(1.35, 4.85, "HAR inputs", ha="center", fontsize=10)
axes[0].add_patch(Rectangle((4.2, 2.35), 1.5, 1.3, facecolor="#e6d8f2", edgecolor="#7b4fa3"))
axes[0].text(4.95, 3.0, "1 hidden\nlayer", ha="center", va="center", fontsize=10)
axes[0].add_patch(Rectangle((7.6, 2.65), 1.25, 0.7, facecolor="#f7d7c4", edgecolor="#c44e52"))
axes[0].text(8.22, 3.0, "forecast", ha="center", va="center", fontsize=9)
axes[0].add_patch(FancyArrowPatch((2.2, 3.0), (4.2, 3.0), arrowstyle="->", mutation_scale=13))
axes[0].add_patch(FancyArrowPatch((5.7, 3.0), (7.6, 3.0), arrowstyle="->", mutation_scale=13))

for i in range(5):
    axes[1].add_patch(Rectangle((0.8 + i * 0.5, 2.55), 0.32, 0.9, facecolor="#d7e5f5", edgecolor="#4c78a8"))
axes[1].text(1.8, 3.8, "22 lags flattened", ha="center", fontsize=10)
axes[1].add_patch(Rectangle((4.2, 2.35), 1.5, 1.3, facecolor="#e6d8f2", edgecolor="#7b4fa3"))
axes[1].text(4.95, 3.0, "1 hidden\nlayer", ha="center", va="center", fontsize=10)
axes[1].add_patch(Rectangle((7.6, 2.65), 1.25, 0.7, facecolor="#f7d7c4", edgecolor="#c44e52"))
axes[1].text(8.22, 3.0, "forecast", ha="center", va="center", fontsize=9)
axes[1].add_patch(FancyArrowPatch((3.6, 3.0), (4.2, 3.0), arrowstyle="->", mutation_scale=13))
axes[1].add_patch(FancyArrowPatch((5.7, 3.0), (7.6, 3.0), arrowstyle="->", mutation_scale=13))

for ax, kind in [(axes[2], "rnn"), (axes[3], "lstm")]:
    ax.text(5, 5.15, "756 steps; final four shown", ha="center", fontsize=10)
    ax.text(0.35, 2.02, "...", ha="center", fontsize=11)
    xs = [0.8, 2.7, 4.6, 6.5]
    for j, x in enumerate(xs):
        lag = len(xs) - 1 - j
        suffix = "t" if lag == 0 else f"t-{lag}"
        lag_label = f"$z_{{{suffix}}}$"
        cell_label = f"$h_{{{suffix}}}$" if kind == "rnn" else f"$h_{{{suffix}}}$\n$C_{{{suffix}}}$"
        ax.add_patch(Rectangle((x, 3.35), 0.9, 0.75, facecolor="#d7e5f5", edgecolor="#4c78a8"))
        ax.text(x + 0.45, 3.73, lag_label, ha="center", va="center", fontsize=8)
        ax.add_patch(Rectangle((x, 1.65), 0.9, 0.75, facecolor="#e6d8f2", edgecolor="#7b4fa3"))
        ax.text(x + 0.45, 2.02, cell_label, ha="center", va="center", fontsize=7)
        ax.add_patch(FancyArrowPatch((x + 0.45, 3.35), (x + 0.45, 2.4), arrowstyle="->", mutation_scale=11))
        if j < len(xs) - 1:
            ax.add_patch(FancyArrowPatch((x + 0.9, 2.02), (xs[j+1], 2.02), arrowstyle="->", mutation_scale=11))
    ax.add_patch(Rectangle((8.35, 1.72), 1.1, 0.6, facecolor="#f7d7c4", edgecolor="#c44e52"))
    ax.text(8.9, 2.02, "forecast", ha="center", va="center", fontsize=8)
    ax.add_patch(FancyArrowPatch((7.4, 2.02), (8.35, 2.02), arrowstyle="->", mutation_scale=11))

plt.tight_layout()
plt.show()

Figure 8.2: Four volatility-forecasting specifications. The HAR-FNN uses the same three HAR components as HAR-OLS: yesterday’s realized variance, the trailing 5-day average, and the trailing 21-day average (in logs for the log target, in levels for the level target). The direct FNN flattens the 22-day window. The Simple RNN and LSTM process a 756-day ordered sequence; the last four steps are shown.

The dense hidden layers use softplus activations, as in the feed-forward chapter. The recurrent layers use \tanh updates, as in the RNN chapter and LSTM chapter. The input shapes are (3,) for HAR-FNN, (22,) for the direct FNN, and (756, 1) for the recurrent models.

Both recurrent models use stateful=False and receive no initial state. The hidden state starts at zero for each window, as does the LSTM cell state. Estimated weights persist between daily forecasts within a year; states do not pass from one input window to the next. The forward state and gradient paths are therefore limited to 756 steps. A long input window makes distant observations available without ensuring that their information survives the recurrent updates.

All neural models minimize mean squared error (MSE) with Adam. We apply global gradient-norm clipping at 1.0 to the stacked raw gradient before Adam’s update, following the RNN chapter’s clipping discussion. Clipping limits large raw gradients; it neither restores vanished gradients nor directly bounds the final Adam parameter step. No weight decay or dropout is used in this illustration.

def make_network(model_name, input_shape, width, learning_rate):
    import keras
    from keras import layers
    keras.backend.clear_session()
    keras.utils.set_random_seed(SEED)
    if model_name in ("FNN", "HAR-FNN"):
        hidden = layers.Dense(width, activation="softplus")
    elif model_name == "Simple RNN":
        hidden = layers.SimpleRNN(width, activation="tanh", stateful=False)
    elif model_name == "LSTM":
        hidden = layers.LSTM(width, activation="tanh", stateful=False)
    else:
        raise ValueError(model_name)
    model = keras.Sequential([keras.Input(shape=input_shape), hidden, layers.Dense(1)])
    # Clip the stacked raw gradient, following the recurrent-network chapter.
    model.compile(loss="mse", optimizer=keras.optimizers.Adam(
        learning_rate=learning_rate, global_clipnorm=1.0))
    return model

Training-only normalization. For raw log inputs, every lag and the target use the current estimation sample’s target mean and standard deviation. For raw level inputs, each lag coordinate is standardized using its own moments across the current training examples. The HAR components are standardized column by column in either pipeline, and the neural target uses its own training mean and standard deviation. At an annual refit, we reestimate these constants from the new sliding sample, then hold them fixed during that year’s forecasting. HAR-OLS is fit directly on the unstandardized HAR inputs; with an intercept, nonsingular affine standardization would leave its fitted predictions unchanged.

The normalization rule differs between the log and level pipelines. Thus their empirical comparison changes the complete training procedure, rather than identifying the effect of the target transform alone. The cross-validation chapter explains why fitting preprocessing on future observations would invalidate the forecast exercise.

def arrays_for(panel, model_name, target):
    prefix = "har" if model_name.startswith("HAR") else "lag"
    suffix = "log" if target == "log(rv)" else "rv"
    X = panel[f"{prefix}_{suffix}"]
    if model_name == "FNN":
        X = X[:, -FNN_INPUT_DAYS:]
    return X, panel[f"y_{suffix}"]

def fit_scaler(X, y, model_name, target):
    """Fit only on the current estimation examples, including their inputs."""
    y_center, y_scale = float(y.mean()), float(y.std())
    if y_scale <= 0:
        raise ValueError("The estimation targets have zero variance.")
    if not model_name.startswith("HAR") and target == "log(rv)":
        x_center, x_scale = y_center, y_scale
    else:
        x_center, x_scale = X.mean(axis=0), X.std(axis=0)
        x_scale = np.where(x_scale > 0, x_scale, 1.0)
    return {"x_center": x_center, "x_scale": x_scale,
            "y_center": y_center, "y_scale": y_scale}

def transform(X, y, scaler, model_name):
    X_scaled = ((X-scaler["x_center"])/scaler["x_scale"]).astype("float32")
    y_scaled = ((y-scaler["y_center"])/scaler["y_scale"]).astype("float32")
    if model_name in ("Simple RNN", "LSTM"):
        X_scaled = X_scaled[..., np.newaxis]
    return X_scaled, y_scaled

8.5 Validation and Tuning

The same input and scaling rules now apply to every candidate configuration. We tune each neural architecture separately for each target using 20 Optuna trials. Optuna is a hyperparameter-search library; its tree-structured Parzen estimator (TPE) selects configurations using earlier validation results. The hyperparameter-optimization chapter develops this search strategy. Here the search ranges are width 16–96 in steps of 4, batch size 16–64 in steps of 4, and learning rate 10^{-4}–10^{-3} on a logarithmic scale.

For each candidate, we fit on the 756 training examples available at the start of 2022, run 30 training epochs, and evaluate each epoch’s fitted weights on all 251 validation targets in 2022. Each epoch is evaluated with one set of weights held fixed across that year; the daily inputs still update using observations available at each forecast origin. We select the configuration and epoch count with the lowest validation MSE, then reuse that epoch count at the start-of-year refits for 2023 and 2024. Validation and testing thus share the same annual fitting schedule. One validation year remains a limited basis for tuning: it does not measure sensitivity to choosing a different validation year, and the selected epoch count need not be optimal for each later training sample.

We use a fixed neural initialization seed of 0 and an Optuna sampler seed of 1. Every model has the same nominal trial and epoch budget, but the search spaces differ in parameter count and optimization difficulty. Each configuration is assessed with one initialization. The comparison therefore does not measure variability across seeds or establish that every architecture has been equally well optimized.

The time-series restrictions concern when outcomes become available and which examples are used to fit, tune, and evaluate the model. For these self-contained, stateless windows, shuffling the already constructed training examples would not introduce future observations. We set shuffle=False to use chronologically ordered mini-batches, a batch-composition choice. The optimization chapter describes mini-batch sampling and reshuffling. We never randomize the train-validation-test split.

def tune_network(panel, model_name, target, trials=N_TRIALS, max_epochs=MAX_EPOCHS):
    import optuna
    optuna.logging.set_verbosity(optuna.logging.WARNING)
    dates = panel["dates"]
    tr = estimation_indices(dates, VALIDATION_START)
    va = np.flatnonzero(np.asarray((dates >= VALIDATION_START) & (dates < TEST_START)))
    X, y = arrays_for(panel, model_name, target)
    scaler = fit_scaler(X[tr], y[tr], model_name, target)
    X_tr, y_tr = transform(X[tr], y[tr], scaler, model_name)
    X_va, y_va = transform(X[va], y[va], scaler, model_name)
    trial_rows = []

    def objective(trial):
        start = time.monotonic()
        params = {
            "learning_rate": trial.suggest_float("learning_rate", 1e-4, 1e-3, log=True),
            "batch_size": trial.suggest_int("batch_size", 16, 64, step=4),
            "width": trial.suggest_int("width", 16, 96, step=4),
        }
        model = make_network(model_name, X_tr.shape[1:], params["width"], params["learning_rate"])
        history = model.fit(X_tr, y_tr, validation_data=(X_va, y_va),
                            epochs=max_epochs, batch_size=params["batch_size"],
                            shuffle=False, verbose=0)
        losses = np.asarray(history.history["val_loss"], dtype=float)
        if not np.isfinite(losses).all():
            raise FloatingPointError("Nonfinite validation loss; investigate before using results.")
        best_epoch = int(np.argmin(losses))+1
        best_loss = float(losses[best_epoch-1])
        trial.set_user_attr("best_epoch", best_epoch)
        trial.set_user_attr("parameters", model.count_params())
        row = dict(trial=trial.number, **params, best_epoch=best_epoch,
                   validation_loss=best_loss, parameters=model.count_params(),
                   seconds=round(time.monotonic()-start, 2))
        trial_rows.append(row)
        print(json.dumps({"stage": "trial", "model": model_name, "target": target, **row}), flush=True)
        return best_loss

    study = optuna.create_study(direction="minimize", sampler=optuna.samplers.TPESampler(seed=1))
    study.optimize(objective, n_trials=trials)
    best = dict(study.best_params, **study.best_trial.user_attrs,
                validation_loss=float(study.best_value))
    return best, trial_rows

The validation-selected configurations are reported below. Validation losses are comparable within a target group; between groups, the standardized targets differ. Parameter counts include the hidden or recurrent layer and the final linear output. HAR-OLS has no neural hyperparameters and is fit by OLS at each annual refit.

Target fit	Model	Width	Parameters	Learning rate	Batch size	Epochs	Validation loss
log(rv)	HAR-FNN	96	481	0.000987	52	5	0.2380
log(rv)	FNN	72	1,729	0.000930	32	11	0.2345
log(rv)	Simple RNN	96	9,505	0.000461	48	4	0.2336
log(rv)	LSTM	80	26,321	0.000653	56	18	0.2260
rv	HAR-FNN	96	481	0.000671	44	15	0.0913
rv	FNN	56	1,345	0.000513	44	26	0.0905
rv	Simple RNN	72	5,401	0.000930	32	3	0.0873
rv	LSTM	72	21,385	0.000930	32	18	0.0936

8.6 Annual Forecasting and Replication

After validation, we freeze the selected hyperparameters and apply the annual sliding-window rule. At the first test origin in each year, we select the most recent 756 observed training targets, reestimate the scalers, and fit each selected network from a fresh initialization. We then forecast daily with those weights until the next refit. The first test fit occurs before the target on 2023-01-03; the next occurs before the target on 2024-01-02. No test loss is used to alter the hyperparameters or refitting rule.

The implementation below reconstructs the sequential forecasts retrospectively. It can predict a whole year in one batch because each row uses only its own past observations and the weights stay fixed within that year. It stores a separate forecast for each target date. Its scoring helper also counts nonpositive variance forecasts; log-scale losses and the volatility loss introduced in the next section require positive predictions and are returned as undefined if that condition fails. Squared and absolute errors on the RV scale remain well defined.

def score_forecasts(rv_true, rv_hat):
    rv_true, rv_hat = np.asarray(rv_true), np.asarray(rv_hat)
    if not np.isfinite(rv_hat).all():
        raise FloatingPointError("Nonfinite forecast.")
    nonpositive = int(np.count_nonzero(rv_hat <= 0))
    if nonpositive:
        rmse_log = mae_log = qlike = None
    else:
        log_error = np.log(rv_hat)-np.log(rv_true)
        ratio = rv_true/rv_hat
        rmse_log = float(np.sqrt(np.mean(log_error**2)))
        mae_log = float(np.mean(np.abs(log_error)))
        qlike = float(np.mean(ratio-np.log(ratio)-1))
    return {"rmse_log_rv": rmse_log, "mae_log_rv": mae_log,
            "rmse_rv": float(np.sqrt(np.mean((rv_hat-rv_true)**2))),
            "mae_rv": float(np.mean(np.abs(rv_hat-rv_true))), "qlike": qlike,
            "nonpositive_rv_forecasts": nonpositive}

def annual_forecasts(panel, model_name, target, params):
    """Refit once per test year on a sliding window; freeze within the year."""
    from sklearn.linear_model import LinearRegression
    dates = panel["dates"]
    X, y = arrays_for(panel, model_name, target)
    predictions, records = [], []
    for year in sorted(dates.loc[dates >= TEST_START].dt.year.unique()):
        te = np.flatnonzero(np.asarray((dates >= TEST_START) & (dates.dt.year == year)))
        first_date = dates.iloc[te[0]]
        tr = estimation_indices(dates, first_date)
        if model_name == "HAR-OLS":
            model = LinearRegression().fit(X[tr], y[tr])
            target_hat = model.predict(X[te])
        else:
            scaler = fit_scaler(X[tr], y[tr], model_name, target)
            X_tr, y_tr = transform(X[tr], y[tr], scaler, model_name)
            X_te, _ = transform(X[te], y[te], scaler, model_name)
            model = make_network(model_name, X_tr.shape[1:], params["width"], params["learning_rate"])
            model.fit(X_tr, y_tr, epochs=params["best_epoch"], batch_size=params["batch_size"],
                      shuffle=False, verbose=0)
            target_hat = model.predict(X_te, verbose=0).ravel()*scaler["y_scale"]+scaler["y_center"]
        rv_hat = np.exp(target_hat) if target == "log(rv)" else target_hat
        with np.errstate(invalid="ignore", divide="ignore"):
            z_hat = np.where(rv_hat > 0, np.log(rv_hat), np.nan)
        for j, index in enumerate(te):
            predictions.append({"date": dates.iloc[index].strftime("%Y-%m-%d"),
                                "rv": float(panel["y_rv"][index]),
                                "forecast_rv": float(rv_hat[j]),
                                "forecast_log_rv": float(z_hat[j]) if np.isfinite(z_hat[j]) else None})
        records.append({"year": int(year), "estimation_examples": len(tr),
                        "estimation_start": dates.iloc[tr[0]].strftime("%Y-%m-%d"),
                        "estimation_end": dates.iloc[tr[-1]].strftime("%Y-%m-%d"),
                        "first_test_target": first_date.strftime("%Y-%m-%d"),
                        "last_test_target": dates.iloc[te[-1]].strftime("%Y-%m-%d"),
                        "forecasts": len(te)})
    frame = pd.DataFrame(predictions)
    return score_forecasts(frame.rv, frame.forecast_rv), predictions, records

All ten model-target combinations are reproduced with the commands below, using the local SPY file documented in the dataset appendix. The Python dependencies are pinned in requirements.txt. Model fitting is not repeated during chapter rendering. The first command saves the selected configurations, trial histories, annual estimation dates, metrics, and daily forecasts in the private data directory. The second verifies the data and pipeline hashes, recomputes metrics from those forecasts, and updates the public aggregate results and figure together.

#| eval: false

python nn_example_pipeline.py --all --jobs 4
python nn_example_outputs.py
quarto render nn_example.qmd

The aggregate results file records the settings, package versions, selected configurations, annual estimation dates, and metrics. The output helper builds both displayed tables from this file and the performance figure from the private forecasts. The data directory remains private; the aggregate file and figure accompany the chapter. The fitting command can resume completed model runs, so use a fresh output directory after a data, code, or environment change rather than mixing old and new results.

8.7 Results

The following results use the common 502 test targets and the annual estimation schedule above. The log-target and level-target rows differ in their target and normalization rules, while preserving the calendar boundaries, input construction, tuning budget, and annual refit schedule.

The columns answer different questions. Write \mathcal T_{\mathrm{test}} for the set of N_{\mathrm{test}}=502 test target dates and \mathcal F_{t-1} for the information available after observing day t-1. Here \text{RMSE}=\sqrt{\frac{1}{N_{\text{test}}}\sum_{t\in\mathcal T_{\mathrm{test}}}(\hat z_t-z_t)^2} aggregates the squared forecast errors over the N_{\text{test}} test origins; \hat z_t denotes the forecast of z_t formed at t-1; we drop the two-index subscript of \hat z_{t+1\mid t} in this section to keep the formulas readable. RMSE on \log(RV) evaluates forecasts on the scale used by the log-target models and is the main comparison criterion here. For level-target models, the log-scale metrics require the level forecast \widehat{RV}_t to be strictly positive so that its logarithm can be taken; they are reported as undefined for any model that violates this. RMSE on RV evaluates the implied variance forecasts after exponentiating log-target forecasts. This back-transformation carries a subtlety: if \hat z_t targets the conditional mean \mu_t:=\mathbb{E}[Z_t\mid\mathcal F_{t-1}] of the random log target Z_t=\log RV_t, then \exp(\hat z_t) is a conditional geometric-mean forecast. It equals the conditional median only under an additional condition that makes the conditional mean and median of Z_t coincide, such as conditional symmetry. Writing Z_t=\mu_t+\varepsilon_t with \mathbb{E}[\varepsilon_t\mid\mathcal F_{t-1}]=0, the conditional mean on the level scale is

\mathbb{E}[RV_t\mid\mathcal F_{t-1}] = \exp(\mu_t)\,\mathbb{E}[\exp(\varepsilon_t)\mid\mathcal F_{t-1}].

The familiar correction \exp(\mu_t+\sigma_t^2/2) follows only when the log-scale error is conditionally Gaussian with variance \sigma_t^2; in general the smearing factor depends on its full conditional distribution. Transforming a complete predictive density is a separate change-of-variables step, discussed in the distributional networks chapter. Because RV-scale squared-error and QLIKE losses are minimized by conditional-mean forecasts, the uncorrected exponentiated log forecasts do not generally target the same level-scale functional as the level-target models. The log-scale RMSE column, the main criterion here, is unaffected. QLIKE (quasi-likelihood) is the volatility loss

\frac{1}{N_{\text{test}}}\sum_{t\in\mathcal T_{\mathrm{test}}} \left( \frac{RV_t}{\widehat{RV}_t} - \log\frac{RV_t}{\widehat{RV}_t} -1 \right),

which is defined only for positive variance forecasts; with the same indexing convention as above, \widehat{RV}_t is the forecast of RV_t formed at t-1. All three measures are losses, so smaller values are better. If a model produces a nonpositive variance forecast, the two losses that require positivity, log-RV RMSE and QLIKE, are reported as undefined rather than rescued by clipping. The table includes the number of nonpositive forecasts; the aggregate results file also records mean absolute error (MAE) on both scales.

Target fit	Model	RMSE log RV	RMSE RV	QLIKE	Nonpositive
log(rv)	HAR-OLS	0.5784	0.3191	0.1880	0
log(rv)	HAR-FNN	0.5836	0.3203	0.1918	0
log(rv)	FNN	0.6257	0.3212	0.1828	0
log(rv)	Simple RNN	0.5866	0.3257	0.1968	0
log(rv)	LSTM	0.5772	0.3216	0.1956	0
rv	HAR-OLS	0.7129	0.3384	0.2048	0
rv	HAR-FNN	0.6201	0.3211	0.1777	0
rv	FNN	0.6223	0.3217	0.1843	0
rv	Simple RNN	0.9063	0.4231	0.2938	0
rv	LSTM	0.6598	0.3257	0.1964	0

The two HAR specifications show that adding nonlinear layers does not uniformly improve the benchmark. With a log target, HAR-OLS has lower loss than HAR-FNN on all three displayed criteria (log-RV RMSE 0.5784 versus 0.5836). With a level target, HAR-FNN improves on HAR-OLS on all three criteria. These comparisons concern the selected complete procedures, not an isolated effect of nonlinearity.

All ten specifications produce positive variance forecasts in this run, so every displayed loss is defined. The linear output layer does not enforce positivity in the level-target models, however. The nonpositive count makes that restriction visible: a positive count renders log-scale RMSE and QLIKE undefined for that row. Clipping to a small arbitrary constant would change both the forecast and its score. A positive output transformation is an alternative specification, which should be chosen before test evaluation.

Figure 8.3 displays the forecasts and compares the two pipelines on the chapter’s main log-scale criterion.

Figure 8.3: Empirical output for the annually refitted SPY experiment, with 22 daily inputs for the direct FNN and 756 for the recurrent models. Top: observed log realized variance and daily forecasts from all five log-target models over the common test period; the dotted line marks the 2024 refit. Bottom: test RMSE on log RV, with blue bars for log-target fits and orange bars for level-target fits. The vertical axis starts at zero; an undefined marker replaces a score whenever a model produces a nonpositive variance forecast.

Among the individual log-target models, LSTM has the lowest test log-RV RMSE (0.5772), with HAR-OLS close behind (0.5784). Their unrounded gap is about 0.0011 log-RV units, or 0.19% of the HAR-OLS RMSE. HAR-FNN follows at 0.5836, the Simple RNN at 0.5866, and the direct FNN at 0.6257. The parsimonious HAR-OLS benchmark therefore remains competitive with the long recurrent input history in this run.

The log-target pipeline has lower log-RV RMSE for both HAR specifications, the Simple RNN, and LSTM. The direct FNN instead slightly favors the level-target pipeline (0.6223 versus 0.6257). Among these individual specifications, secondary criteria give different rankings: log-target HAR-OLS has the lowest RV-scale RMSE (0.3191), while level-target HAR-FNN has the lowest QLIKE (0.1777). Because normalization also changes between target pipelines, these results do not isolate the effect of taking logs of the target.

The common test period makes the forecast errors comparable, but the information sets still differ. The direct FNN sees 22 daily inputs, while the recurrent models see 756. The figure therefore cannot identify the benefit of recurrence holding the input history fixed. Nor can a result from one initialization establish that an architecture is generally better optimized or more accurate.

Numerical rankings do not establish statistically significant performance differences. Let \mathcal T_{\mathrm{test}} be the set of 502 test target dates and \delta_t the difference between two models’ squared forecast errors at date t. A Diebold and Mariano (1995) comparison relates the mean differential \bar\delta=N_{\mathrm{test}}^{-1}\sum_{t\in\mathcal T_{\mathrm{test}}}\delta_t to its standard error, (\hat\omega^2/N_{\mathrm{test}})^{1/2}, where \hat\omega^2 is a heteroskedasticity- and autocorrelation-consistent (HAC) estimate of the differential’s long-run variance. The relevant uncertainty is the standard error of the mean, not the dispersion of individual daily differences. Exercise 4.3 in the predictive-distributions chapter develops the corresponding construction for score differentials. We do not run this test here, and the results do not establish equivalence between close performers.

QLIKE is the Patton (2011) quasi-likelihood loss for volatility. Its robustness result concerns rankings by expected loss: under the required conditional-unbiasedness condition for the variance proxy, replacing latent variance by the proxy preserves the expected-loss ranking. Sampling variation can still reverse finite-sample rankings. QLIKE is scale-invariant because it depends only on the ratio RV_t/\widehat{RV}_t; it penalizes underprediction more heavily than overprediction. Thus QLIKE and log-RV RMSE can favor different forecasts.

HAR remains a strong benchmark in this comparison. A longer raw history, nonlinear mapping, or gated recurrent state does not by itself guarantee improved forecasts, and the ranking cannot be assigned to memory or nonlinearity alone.

8.8 Summary

Key Takeaways

1. Log-target LSTM has the lowest test log-RV RMSE among the individual models, only slightly below HAR-OLS.
2. Input history differs across models: 22 daily observations for the direct FNN, 756 for the recurrent models, and three HAR summaries.
3. Annual refits keep 756 training examples, while the neural hyperparameters stay fixed.
4. Empirical rankings depend on the loss and complete forecasting procedure; numerical gaps alone do not establish statistical significance.

Common Pitfalls

• Allowing the annual estimation sample to expand changes the intended forecasting rule.
• Attributing performance differences to recurrence alone ignores the different input histories across models.
• Clipping negative variance forecasts to an arbitrary positive constant changes the submitted forecast and its score.
• Choosing hyperparameters after inspecting test performance turns the test period into additional validation data.

8.9 References

Andersen, Torben G., Tim Bollerslev, Francis X. Diebold, and Paul Labys. 2003. “Modeling and Forecasting Realized Volatility.” Econometrica 71 (2): 579–625. https://doi.org/10.1111/1468-0262.00418.

Corsi, Fulvio. 2009. “A Simple Approximate Long-Memory Model of Realized Volatility.” Journal of Financial Econometrics 7 (2): 174–96. https://doi.org/10.1093/jjfinec/nbp001.

Diebold, Francis X, and Robert S Mariano. 1995. “Comparing predictive accuracy.” Journal of Business and Economic Statistics 13 (3): 253–63. https://doi.org/10.1080/07350015.1995.10524599.

Patton, Andrew J. 2011. “Volatility Forecast Comparison Using Imperfect Volatility Proxies.” Journal of Econometrics 160 (1): 246–56. https://doi.org/10.1016/j.jeconom.2010.03.034.

Footnotes

1. Under standard semimartingale conditions, realized variance consistently estimates quadratic variation. In a continuous-path model without price jumps, quadratic variation equals integrated variance; with jumps, it also contains the sum of squared jump sizes.