Under active development. A stable version is expected in November 2026. Feedback is welcome by email or via GitHub issues.
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 complete log-target pipeline with a complete level-target pipeline. In this sample the log transform removes most of the target’s right skewness, and the log-target pipeline performs better on the chapter’s main log-scale criterion. The comparison also changes the normalization of the lag inputs, however, so it does not isolate a causal effect of the target transformation alone. Second, it compares a linear benchmark—the heterogeneous autoregressive (HAR) model of daily, weekly, and monthly volatility components, estimated by ordinary least squares (OLS), abbreviated HAR-OLS—with four neural-network specifications: a parsimonious HAR-FNN (a feed-forward neural network, FNN, on the HAR inputs), a full FNN, a Simple RNN (recurrent neural network), and an LSTM—all evaluated under the same forecast-origin discipline. A central finding is that the linear benchmark is hard to beat: with the regressor set held fixed, the nonlinear HAR-FNN improves on HAR-OLS only marginally; among the full-history networks, only the full FNN edges past HAR-OLS—by well under one percent of RMSE—while the recurrent models do not beat it.
8.2 Roadmap
- We define the realized-variance target and the one-step-ahead forecast origin.
- We specify the local-data workflow and the time-ordered train-validation-test split.
- We define the inputs for HAR-OLS, HAR-FNN, the full FNN, the Simple RNN, and the LSTM.
- We tune the neural networks using only the validation block.
- We interpret the empirical results as evidence about target-and-normalization pipelines, benchmark choice, and architecture choice.
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. Realized variance RV_t is a canonical high-frequency volatility measure; Andersen et al. (2003) is a standard reference for modeling and forecasting it.1 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 log target, the one-step-ahead forecasting problem is
\widehat z_{t+1\mid t} = f_\theta(z_t,z_{t-1},\ldots,z_{t-21}),
where the input window contains the previous 22 trading days, roughly one trading month. The forecast target z_{t+1} is not part of the input window.
The full-history models all use this same forecast-origin history. The full FNN receives the 22 lags as a flattened feature vector. The Simple RNN and LSTM receive the same values as an ordered sequence. The HAR specifications are different by design: they receive a traditional volatility summary consisting of yesterday’s realized variance RV_t, the trailing 5-day average, and the trailing 21-day average. This three-component (daily, weekly, monthly) decomposition is the HAR model of Corsi (2009), designed to capture the long-memory structure of realized variance with a parsimonious linear specification. All three components are known at the forecast origin t.
The local replication code expects the daily SPY realized-measures file documented in the dataset appendix.
Figure 8.1 shows the forecast-origin convention. The training, validation, and test blocks are chronological. Scaling constants and hyperparameters are learned before the test block is touched.
The workflow below makes the econometric design explicit: it constructs forecast origins before splitting the sample, estimates scaling constants only on the training block, and keeps the test block untouched until the final comparison. Concretely, the 22-lag inputs and the target are transformed with the training-block mean and standard deviation of the target series—every lag coordinate is the same variable, so a single center and scale serves all of them—while the three HAR components are standardized column by column with their own training-block moments (cf. the Cross-Validation chapter on fitting preprocessing on training data only); forecasts are mapped back to the original scale. Replication uses the local SPY realized-variance file documented in the dataset appendix.
The first chunk loads the realized-variance series and constructs z_t=\log(RV_t).
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 optuna
import keras
from keras import layers
from sklearn.metrics import mean_absolute_error, root_mean_squared_error
from sklearn.linear_model import LinearRegression
optuna.logging.set_verbosity(optuna.logging.WARNING)
random.seed(0)
np.random.seed(0)
keras.utils.set_random_seed(0)
path = "data/taq_spy/SPY_daily_measures.csv"
rv_source = (
pd.read_csv(path, 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)
# Guard: no missing or nonpositive rv inside the usable sample, so that lag k
# below means k trading days back rather than k available observations back.
usable = rv_source["date"] <= rv["date"].max()
assert (rv_source.loc[usable, "rv"] > 0).all(), "gap inside the usable sample"
rv["z"] = np.log(rv["rv"])
print(f"observations: {len(rv)}, range: {rv['date'].min().date()} to {rv['date'].max().date()}")observations: 2516, range: 2015-01-02 to 2024-12-31The next helper constructs the supervised forecasting problem. For each valid origin t, the full-history models see (z_t,\ldots,z_{t-21}) and forecast z_{t+1}. The HAR helper uses the same origin t but replaces the 22 individual lags by three realized-variance summaries:
\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).
For the log-target models these three components enter in logs—the same regressors as HAR-OLS, matching the scale of the target—while the level-target models use them in levels.
The training block is the first 64% of valid forecast origins, the validation block is the next 16%, and the test block is the final 20%. The 64/16/20 allocation is a design choice, not an optimum: it reserves a contiguous validation block for tuning while keeping a sizable test block, and we do not investigate sensitivity to this split here. Optuna is a hyperparameter-search library: a trial is one sampled hyperparameter configuration, evaluated on the validation objective, and its default sampler is the tree-structured Parzen estimator (TPE), a Bayesian search strategy. See the chapter on hyperparameter optimization for the underlying selection-bias and search-strategy issues.
Show the code
window = 22
def make_windows(values, window):
X, y, target_index = [], [], []
for end in range(window - 1, len(values) - 1):
X.append(values[end - window + 1 : end + 1])
y.append(values[end + 1])
target_index.append(end + 1)
return np.asarray(X), np.asarray(y), np.asarray(target_index)
def make_har_features(rv_values, target_values, window):
X, y, target_index = [], [], []
for end in range(window - 1, len(rv_values) - 1):
X.append(
[
rv_values[end],
rv_values[end - 4 : end + 1].mean(),
rv_values[end - 20 : end + 1].mean(),
]
)
y.append(target_values[end + 1])
target_index.append(end + 1)
return np.asarray(X), np.asarray(y), np.asarray(target_index)
def make_log_har_ols_features(rv_values, window):
X, y, target_index = [], [], []
for end in range(window - 1, len(rv_values) - 1):
X.append(
np.log(
[
rv_values[end],
rv_values[end - 4 : end + 1].mean(),
rv_values[end - 20 : end + 1].mean(),
]
)
)
y.append(np.log(rv_values[end + 1]))
target_index.append(end + 1)
return np.asarray(X), np.asarray(y), np.asarray(target_index)
X_raw, y_raw, target_index = make_windows(rv["z"].to_numpy(), window)
X_har_ols_log_raw, y_har_ols_log_raw, _ = make_log_har_ols_features(
rv["rv"].to_numpy(),
window,
)
X_har_raw = X_har_ols_log_raw # log HAR components: identical regressors to HAR-OLS
X_har_rv_raw, y_har_rv_raw, _ = make_har_features(
rv["rv"].to_numpy(),
rv["rv"].to_numpy(),
window,
)
X_lag_rv_raw, y_lag_rv_raw, _ = make_windows(rv["rv"].to_numpy(), window)
dates = rv.loc[target_index, "date"].reset_index(drop=True)
n = len(y_raw)
train_end = int(0.64 * n)
valid_end = int(0.80 * n)
X_train_raw, y_train_raw = X_raw[:train_end], y_raw[:train_end]
X_valid_raw, y_valid_raw = X_raw[train_end:valid_end], y_raw[train_end:valid_end]
X_test_raw, y_test = X_raw[valid_end:], y_raw[valid_end:]
test_dates = dates.iloc[valid_end:].reset_index(drop=True)
X_har_train_raw = X_har_raw[:train_end]
X_har_valid_raw = X_har_raw[train_end:valid_end]
X_har_test_raw = X_har_raw[valid_end:]
X_har_ols_log_train_valid_raw = X_har_ols_log_raw[:valid_end]
X_har_ols_log_test_raw = X_har_ols_log_raw[valid_end:]
y_har_ols_log_train_valid = y_har_ols_log_raw[:valid_end]
X_har_rv_train_raw = X_har_rv_raw[:train_end]
X_har_rv_valid_raw = X_har_rv_raw[train_end:valid_end]
X_har_rv_test_raw = X_har_rv_raw[valid_end:]
y_rv_train_raw = y_har_rv_raw[:train_end]
y_rv_valid_raw = y_har_rv_raw[train_end:valid_end]
y_rv_test = y_har_rv_raw[valid_end:]
X_har_ols_rv_train_valid_raw = X_har_rv_raw[:valid_end]
X_har_ols_rv_test_raw = X_har_rv_raw[valid_end:]
y_har_ols_rv_train_valid = y_har_rv_raw[:valid_end]
X_lag_rv_train_raw = X_lag_rv_raw[:train_end]
X_lag_rv_valid_raw = X_lag_rv_raw[train_end:valid_end]
X_lag_rv_test_raw = X_lag_rv_raw[valid_end:]
center = y_train_raw.mean()
scale = y_train_raw.std()
x_har_center = X_har_train_raw.mean(axis=0)
x_har_scale = X_har_train_raw.std(axis=0)
x_har_scale = np.where(x_har_scale == 0, 1, x_har_scale)
x_har_rv_center = X_har_rv_train_raw.mean(axis=0)
x_har_rv_scale = X_har_rv_train_raw.std(axis=0)
x_har_rv_scale = np.where(x_har_rv_scale == 0, 1, x_har_rv_scale)
y_rv_center = y_rv_train_raw.mean()
y_rv_scale = y_rv_train_raw.std()
x_lag_rv_center = X_lag_rv_train_raw.mean(axis=0)
x_lag_rv_scale = X_lag_rv_train_raw.std(axis=0)
x_lag_rv_scale = np.where(x_lag_rv_scale == 0, 1, x_lag_rv_scale)
def standardize(a):
return (a - center) / scale
def unstandardize(a):
return a * scale + center
def standardize_har(a):
return (a - x_har_center) / x_har_scale
def standardize_har_rv(a):
return (a - x_har_rv_center) / x_har_rv_scale
def unstandardize_rv(a):
return a * y_rv_scale + y_rv_center
def standardize_lag_rv(a):
return (a - x_lag_rv_center) / x_lag_rv_scale
X_train = standardize(X_train_raw)
X_valid = standardize(X_valid_raw)
X_test = standardize(X_test_raw)
X_har_train = standardize_har(X_har_train_raw)
X_har_valid = standardize_har(X_har_valid_raw)
X_har_test = standardize_har(X_har_test_raw)
y_train = standardize(y_train_raw)
y_valid = standardize(y_valid_raw)
X_har_rv_train = standardize_har_rv(X_har_rv_train_raw)
X_har_rv_valid = standardize_har_rv(X_har_rv_valid_raw)
X_har_rv_test = standardize_har_rv(X_har_rv_test_raw)
y_rv_train = (y_rv_train_raw - y_rv_center) / y_rv_scale
y_rv_valid = (y_rv_valid_raw - y_rv_center) / y_rv_scale
X_lag_rv_train = standardize_lag_rv(X_lag_rv_train_raw)
X_lag_rv_valid = standardize_lag_rv(X_lag_rv_valid_raw)
X_lag_rv_test = standardize_lag_rv(X_lag_rv_test_raw)
print(f"forecast origins: {n}")
print(f"train: {train_end}, validation: {valid_end - train_end}, test: {n - valid_end}")forecast origins: 2494
train: 1596, validation: 399, test: 4998.4 Comparable Specifications
The specifications differ in how they represent the same forecasting problem.
- The HAR-OLS benchmark is linear. For the log target, it regresses \log(RV_{t+1}) on the logs of the daily, weekly, and monthly HAR components. For the level target, it regresses RV_{t+1} on the corresponding level components.
- The HAR-FNN is the most structured specification. It uses the three HAR volatility-component summaries—daily, weekly, and monthly averages of realized variance—from the Corsi (2009) model. For the log target it receives exactly the same three log-transformed components as HAR-OLS (standardized), so the two specifications share their target and predictor information while contrasting a linear OLS procedure with a tuned nonlinear neural-network procedure; for the level target, both receive the level components.
- The full FNN uses the full 22-day history, processing the fixed-position lag vector jointly through dense layers: each lag coordinate receives its own weights, with no recurrent parameter sharing and no sequential state update.
- The Simple RNN also uses the full 22-day history, but processes the lags sequentially through a recurrent state.
- The LSTM processes the same sequence with gated memory, allowing the network to regulate how much past information is retained.
Figure 8.2 summarizes the neural specifications. All models share the forecast origin, the one-step horizon, and roughly one month of history, but they do not receive identical information: the HAR specifications see three summaries of the past 21 days, while the full-history models see all 22 individual lags—and summarization discards information. Each model’s width is tuned separately, so parameter counts differ as well. The most controlled comparison is HAR-OLS versus HAR-FNN, which share the target and predictor information exactly; the other contrasts mix information design with architecture. Even this comparison is not a pure functional-form experiment: HAR-OLS uses direct least-squares fitting on raw HAR inputs, whereas HAR-FNN uses standardized inputs, stochastic optimization, and validation-based tuning. It therefore asks whether the complete tuned nonlinear neural-network procedure improves on the linear benchmark for the same information set.
All four neural networks are deliberately shallow. The HAR-FNN and full FNN have one hidden dense layer with a tuned width and a linear output layer. The Simple RNN and LSTM have one recurrent layer with a tuned number of units and a linear output layer. A deeper pyramidal FNN—for example, one with 22 inputs, hidden layers of 64 and 32 units, and one output—would be a reasonable robustness check. It is not the baseline here because it would add another capacity choice while the chapter’s main comparison is about target transformation and information representation.
Each architecture below uses the same skeleton—one hidden layer of tuned width followed by a linear output. The input shapes enforce the intended information sets: three HAR features for the HAR-FNN, 22 flattened lags for the full FNN, and a 22-step sequence for the recurrent models. The dense models use softplus hidden activations, as in the feed-forward chapter; the recurrent models use the \tanh cells whose gradient behavior is analyzed in Chapter 6 and Chapter 7.
Show the code
def make_fnn(n_neurons, learning_rate):
model = keras.Sequential(
[
keras.Input(shape=(window,)),
layers.Dense(n_neurons, activation="softplus"),
layers.Dense(1),
],
name="FNN",
)
model.compile(
loss=keras.losses.MeanSquaredError(),
optimizer=keras.optimizers.Adam(learning_rate=learning_rate),
)
return model
def make_har_fnn(n_neurons, learning_rate):
model = keras.Sequential(
[
keras.Input(shape=(3,)),
layers.Dense(n_neurons, activation="softplus"),
layers.Dense(1),
],
name="HARFNN",
)
model.compile(
loss=keras.losses.MeanSquaredError(),
optimizer=keras.optimizers.Adam(learning_rate=learning_rate),
)
return model
def make_simple_rnn(n_neurons, learning_rate):
model = keras.Sequential(
[
keras.Input(shape=(window, 1)),
layers.SimpleRNN(
n_neurons,
activation="tanh",
return_sequences=False,
stateful=False,
),
layers.Dense(1),
],
name="SimpleRNN",
)
model.compile(
loss=keras.losses.MeanSquaredError(),
optimizer=keras.optimizers.Adam(learning_rate=learning_rate),
)
return model
def make_lstm(n_neurons, learning_rate):
model = keras.Sequential(
[
keras.Input(shape=(window, 1)),
layers.LSTM(
n_neurons,
activation="tanh",
return_sequences=False,
stateful=False,
),
layers.Dense(1),
],
name="LSTM",
)
model.compile(
loss=keras.losses.MeanSquaredError(),
optimizer=keras.optimizers.Adam(learning_rate=learning_rate),
)
return modelThe code presents each overlapping 22-day window as a separate input sequence and sets stateful=False for both recurrent layers. Because no initial_state is supplied, Keras initializes the Simple RNN hidden state at h_0=0 and the LSTM hidden and cell states at h_0=C_0=0 for every window. The estimated weights persist across windows, but the final hidden and cell states do not: they are discarded rather than passed to the next window. Both the forward state horizon and the gradient horizon are therefore limited to 22 steps. This is the fixed-window training scheme defined in the RNN chapter.
8.5 Validation and Tuning
Hyperparameters are selected using the validation block, not the test block. This is not a nested evaluation design; it is a compact empirical illustration with a single holdout test period. Because the validation block is consulted across trials and epochs, the best attained validation loss is an adaptively selected criterion rather than an unbiased estimate of future performance. The final test block is used once after the selected hyperparameters have been fixed.
The load-bearing time-series restrictions are the chronological train-validation-test split and the look-ahead-free window construction—not the shuffle flag. Under the fixed-window scheme just described, each overlapping window is a self-contained input-target pair, so shuffling pre-built training windows between epochs (the Keras default for stateless models) neither leaks future information nor breaks the forecast-origin discipline. We set shuffle=False as a deliberate batch-composition choice—training batches are then contiguous, chronologically ordered blocks, whose gradient-variance implications under serial correlation are exactly those analyzed in Exercise 3.2 of the optimization chapter—not as a leakage protection. What must never be randomized is the split itself: this chapter does not use random K-fold cross-validation, because random folds would mix volatility regimes across training and validation sets and would no longer mimic the one-step-ahead forecasting problem.
Show the code
def as_sequence(X):
return X[..., np.newaxis]
def tune_model(model_factory, X_tr, X_va, y_tr, y_va, n_trials=20):
def objective(trial):
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),
"n_neurons": trial.suggest_int("n_neurons", 16, 96, step=4),
"epochs": 30,
}
model = model_factory(params["n_neurons"], params["learning_rate"])
history = model.fit(
X_tr,
y_tr,
validation_data=(X_va, y_va),
epochs=params["epochs"],
batch_size=params["batch_size"],
verbose=0,
shuffle=False,
)
val_losses = history.history["val_loss"]
trial.set_user_attr("best_epoch", int(np.argmin(val_losses)) + 1)
return min(val_losses)
study = optuna.create_study(
direction="minimize",
sampler=optuna.samplers.TPESampler(seed=1),
)
study.optimize(objective, n_trials=n_trials)
params = dict(study.best_params)
params["best_epoch"] = study.best_trial.user_attrs["best_epoch"]
return params
best_har_fnn = tune_model(make_har_fnn, X_har_train, X_har_valid, y_train, y_valid)
best_fnn = tune_model(make_fnn, X_train, X_valid, y_train, y_valid)
best_rnn = tune_model(make_simple_rnn, as_sequence(X_train), as_sequence(X_valid), y_train, y_valid)
best_lstm = tune_model(make_lstm, as_sequence(X_train), as_sequence(X_valid), y_train, y_valid)In the local run that produced this chapter’s output, chronological validation selected the following hyperparameters. Validation losses are comparable within a target-fit group but not across the log(rv) and rv target groups: both targets are standardized to unit training variance, but the two groups minimize mean squared error for different standardized targets, so the numbers measure performance on different objectives. The trainable-parameter count includes the hidden or recurrent layer and the final linear output layer; it makes clear that an equal number of tuning trials does not imply equal search difficulty across architectures.
HAR-OLS does not appear in this table because it has no neural tuning parameters. It is fit by ordinary least squares on the combined train-validation block before the final test evaluation.
| Target fit | Model | Tuned width | Trainable parameters | Learning rate | Batch size | Best epoch | Validation loss |
|---|---|---|---|---|---|---|---|
log(rv) |
HAR-FNN | 96 | 481 | 0.000972 | 48 | 30 | 0.2762 |
log(rv) |
FNN | 72 | 1,729 | 0.000346 | 36 | 28 | 0.2743 |
log(rv) |
Simple RNN | 20 | 461 | 0.000752 | 60 | 13 | 0.2717 |
log(rv) |
LSTM | 72 | 21,385 | 0.000715 | 16 | 30 | 0.2678 |
rv |
HAR-FNN | 52 | 261 | 0.000538 | 44 | 22 | 0.1099 |
rv |
FNN | 96 | 2,305 | 0.000972 | 48 | 21 | 0.1156 |
rv |
Simple RNN | 88 | 8,009 | 0.000109 | 24 | 2 | 0.1168 |
rv |
LSTM | 36 | 5,509 | 0.000984 | 64 | 28 | 0.1149 |
The log-target full-history models have similar validation losses, so the final test comparison should not be read as a decisive ranking of recurrent memory versus a flattened lag vector. The level-target validation losses also cluster tightly. Two entries deserve attention already here. The level-target Simple RNN’s validation-selected epoch count is just 2, paired with the smallest learning rate in the table. This is an observed outcome of the tuning rule, not by itself evidence that the model is undertrained; learning curves or repeated-seed fits would be needed for that diagnosis. Two of the log-target winners, the HAR-FNN and the LSTM, sit exactly at the 30-epoch cap, so for them the best-epoch rule did not stop training early; a longer cap might have selected differently. Validation loss alone does not settle the target-transformation question because each panel is fit and scaled on a different target.
After tuning, each selected model is refit on the combined train-validation block and evaluated once on the test block. Only the network weights are refit: the preprocessing moments remain fixed at their estimates from the original training block. This is leakage-safe because the test block remains untouched and preserves the transformation used during tuning, although it is not a full re-estimation of every data-dependent component on the combined pre-test sample. Each row in the final table is therefore the out-of-sample performance of one target choice and one architecture under this stated protocol.
Show the code
def qlike(rv_true, rv_pred):
ratio = rv_true / rv_pred
return np.mean(ratio - np.log(ratio) - 1)
def summarize_forecast(name, target, z_hat, rv_hat):
rv_true = np.exp(y_test)
n_nonpos = int(np.sum(rv_hat <= 0))
# Log-scale and QLIKE losses are defined only for strictly positive variance
# forecasts. If a model submits a nonpositive variance, we report them as
# undefined rather than rescuing the model with an arbitrary positive floor.
if n_nonpos:
rmse_log = mae_log = qlike_val = float("nan")
else:
z_for_log_loss = np.log(rv_hat)
rmse_log = root_mean_squared_error(y_test, z_for_log_loss)
mae_log = mean_absolute_error(y_test, z_for_log_loss)
qlike_val = qlike(rv_true, rv_hat)
return {
"target": target,
"model": name,
"rmse_log_rv": rmse_log,
"mae_log_rv": mae_log,
"rmse_rv": root_mean_squared_error(rv_true, rv_hat),
"mae_rv": mean_absolute_error(rv_true, rv_hat),
"qlike": qlike_val,
"nonpositive_rv_forecasts": n_nonpos,
"forecast": z_hat,
}
def fit_har_ols_log():
model = LinearRegression().fit(
X_har_ols_log_train_valid_raw,
y_har_ols_log_train_valid,
)
z_hat = model.predict(X_har_ols_log_test_raw)
return summarize_forecast("HAR-OLS", "log(rv)", z_hat, np.exp(z_hat))
def fit_har_ols_level():
model = LinearRegression().fit(
X_har_ols_rv_train_valid_raw,
y_har_ols_rv_train_valid,
)
rv_hat = model.predict(X_har_ols_rv_test_raw)
z_hat = np.log(np.maximum(rv_hat, 1e-12))
return summarize_forecast("HAR-OLS", "rv", z_hat, rv_hat)
def fit_and_score(name, model_factory, params, X_trva, y_trva, X_te):
model = model_factory(params["n_neurons"], params["learning_rate"])
model.fit(
X_trva,
y_trva,
epochs=params["best_epoch"],
batch_size=params["batch_size"],
verbose=0,
shuffle=False,
)
z_hat = unstandardize(model.predict(X_te).ravel())
rv_hat = np.exp(z_hat)
return summarize_forecast(name, "log(rv)", z_hat, rv_hat)
X_train_valid = np.vstack([X_train, X_valid])
X_har_train_valid = np.vstack([X_har_train, X_har_valid])
y_train_valid = np.concatenate([y_train, y_valid])
results = [
fit_har_ols_log(),
fit_and_score("HAR-FNN", make_har_fnn, best_har_fnn, X_har_train_valid, y_train_valid, X_har_test),
fit_and_score("FNN", make_fnn, best_fnn, X_train_valid, y_train_valid, X_test),
fit_and_score("Simple RNN", make_simple_rnn, best_rnn, as_sequence(X_train_valid), y_train_valid, as_sequence(X_test)),
fit_and_score("LSTM", make_lstm, best_lstm, as_sequence(X_train_valid), y_train_valid, as_sequence(X_test)),
]
performance = pd.DataFrame(results).drop(columns=["forecast"])
performanceThe level-target rows of the results below use the analogous pipeline on the rv-scale arrays constructed earlier—the changes are the target arrays, the rv-scale standardizers, and the back-transformation from standardized levels rather than logs. The 22 log-lag coordinates use one common center and scale, whereas the 22 level-lag coordinates are standardized coordinate by coordinate. Both transformations are estimated on the training block, but the log-versus-level comparison therefore changes input normalization as well as the target transformation and should not be read as isolating the target scale alone:
Show the code
def fit_and_score_level(name, model_factory, params, X_trva, y_trva, X_te):
model = model_factory(params["n_neurons"], params["learning_rate"])
model.fit(
X_trva,
y_trva,
epochs=params["best_epoch"],
batch_size=params["batch_size"],
verbose=0,
shuffle=False,
)
rv_hat = unstandardize_rv(model.predict(X_te).ravel())
z_hat = np.log(np.maximum(rv_hat, 1e-12))
return summarize_forecast(name, "rv", z_hat, rv_hat)
best_har_fnn_rv = tune_model(make_har_fnn, X_har_rv_train, X_har_rv_valid, y_rv_train, y_rv_valid)
best_fnn_rv = tune_model(make_fnn, X_lag_rv_train, X_lag_rv_valid, y_rv_train, y_rv_valid)
best_rnn_rv = tune_model(make_simple_rnn, as_sequence(X_lag_rv_train), as_sequence(X_lag_rv_valid), y_rv_train, y_rv_valid)
best_lstm_rv = tune_model(make_lstm, as_sequence(X_lag_rv_train), as_sequence(X_lag_rv_valid), y_rv_train, y_rv_valid)
X_har_rv_train_valid = np.vstack([X_har_rv_train, X_har_rv_valid])
X_lag_rv_train_valid = np.vstack([X_lag_rv_train, X_lag_rv_valid])
y_rv_train_valid = np.concatenate([y_rv_train, y_rv_valid])
results_level = [
fit_har_ols_level(),
fit_and_score_level("HAR-FNN", make_har_fnn, best_har_fnn_rv, X_har_rv_train_valid, y_rv_train_valid, X_har_rv_test),
fit_and_score_level("FNN", make_fnn, best_fnn_rv, X_lag_rv_train_valid, y_rv_train_valid, X_lag_rv_test),
fit_and_score_level("Simple RNN", make_simple_rnn, best_rnn_rv, as_sequence(X_lag_rv_train_valid), y_rv_train_valid, as_sequence(X_lag_rv_test)),
fit_and_score_level("LSTM", make_lstm, best_lstm_rv, as_sequence(X_lag_rv_train_valid), y_rv_train_valid, as_sequence(X_lag_rv_test)),
]The final table and figure are assembled directly from the two result lists:
Show the code
import matplotlib.pyplot as plt
all_results = results + results_level
performance = pd.DataFrame(all_results).drop(columns=["forecast"])
fig, (ax_l, ax_r) = plt.subplots(1, 2, figsize=(16, 6), width_ratios=[1.15, 1])
series_colors = {"HAR-OLS": "0.5", "HAR-FNN": "tab:brown", "FNN": "tab:blue",
"Simple RNN": "tab:orange", "LSTM": "tab:green"}
ax_l.plot(test_dates, y_test, color="black", linewidth=1.1, label="Realized log RV")
for r in results:
ax_l.plot(test_dates, r["forecast"], color=series_colors[r["model"]],
linewidth=0.8, label=r["model"])
ax_l.set_title("Log-target forecasts")
ax_l.set_xlabel("Date")
ax_l.set_ylabel("log realized variance")
ax_l.legend(frameon=False, ncol=2, fontsize=9)
models = ["HAR-OLS", "HAR-FNN", "FNN", "Simple RNN", "LSTM"]
log_rmse = [next(r for r in results if r["model"] == m)["rmse_log_rv"] for m in models]
rv_rmse = [next(r for r in results_level if r["model"] == m)["rmse_log_rv"] for m in models]
xpos = np.arange(len(models))
ax_r.bar(xpos - 0.19, log_rmse, 0.38, color="steelblue", label="Fit log(rv)")
ax_r.bar(xpos + 0.19, rv_rmse, 0.38, color="tab:orange", label="Fit rv")
ax_r.set_xticks(xpos)
ax_r.set_xticklabels(models)
ax_r.set_ylabel("test RMSE on log realized variance")
ax_r.set_title("Common-scale target comparison")
ax_r.legend(frameon=False)
ax_r.set_ylim(0, 0.85)
for bars in ax_r.containers:
for b in bars:
h = b.get_height()
if np.isnan(h):
ax_r.annotate("undefined\n(nonpositive\nforecasts)",
(b.get_x() + b.get_width() / 2, 0.02),
ha="center", va="bottom", fontsize=8, color="darkred")
else:
ax_r.annotate(f"{h:.3f}", (b.get_x() + b.get_width() / 2, h),
ha="center", va="bottom", fontsize=9)
ax_l.grid(alpha=0.3)
fig.text(0.01, 0.015,
"Derived from local SPY realized-variance data through 2024-12-31; see the dataset appendix.\n"
"HAR-OLS and HAR-FNN use the same components (logs for the log target, levels for the level target).\n"
"Recurrent models use tanh cells; validation-selected epoch counts are reused in the final refits.\n"
"Log-scale and QLIKE losses are undefined when a model submits nonpositive variance forecasts.",
fontsize=9, color="0.3")
plt.tight_layout(rect=[0, 0.10, 1, 1])
fig.savefig("figures/nn-example-performance.png", dpi=150)8.6 Results
The model-fitting and tuning code above is shown for exposition but is not executed during rendering; the results below come from a local run with access to the Wharton Research Data Services (WRDS) Trade and Quote (TAQ) database. The empirical output uses the 2015-01-02 through 2024-12-31 SPY realized-variance file described in the dataset appendix. The log-target rows come from the workflow shown above. The level-target rows repeat the same split, input construction, tuning budget, and evaluation protocol, but fit the models to RV_{t+1} rather than \log(RV_{t+1}), using the level pipeline shown above. The tuning run used 20 Optuna trials for each neural architecture. This is an equal nominal trial budget, not equal effective search coverage: architectures differ in parameter count, search geometry, and optimization variability. Each candidate is evaluated using one training run with a fixed random seed, so architecture rankings remain conditional on that run; a stronger empirical comparison would repeat the leading configurations across several seeds. Figure 8.3 is likewise a pre-generated figure from the same local run. The Python package versions are pinned in requirements.txt, but no versioned file of the test forecasts accompanies this frozen output; after a data, code, or environment change, the displayed table and figure must therefore be regenerated together rather than assumed current.
One protocol detail deserves explicit mention. During tuning, each configuration is scored by the lowest validation loss attained during its first 30 epochs—an implicit form of early stopping. The epoch count that achieved the selected validation loss is recorded for the winning configuration and reused in the final refit on the combined training and validation data, so the deployed training length matches the one that was selected. Validation selects this epoch count; it does not certify that the count is optimal for the refit sample. One approximation remains: the refit applies that epoch count to a larger sample, and with the validation block folded into training there is nothing left to monitor—an explicit early-stopping callback with a further holdout block would be the fully adaptive alternative.
The columns answer different questions. Here the root mean squared error, \text{RMSE}=\sqrt{\frac{1}{N_{\text{test}}}\sum_t(\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—for quasi-likelihood—is the volatility loss
\frac{1}{N_{\text{test}}}\sum_{t} \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 replication object records the count of such violations and RV-scale absolute error; the compact table below omits those two diagnostics.
| Target fit | Model | RMSE log RV | RMSE RV | QLIKE |
|---|---|---|---|---|
log(rv) |
HAR-OLS | 0.5771 | 0.3183 | 0.1957 |
log(rv) |
HAR-FNN | 0.5771 | 0.3166 | 0.1836 |
log(rv) |
FNN | 0.5735 | 0.3145 | 0.1860 |
log(rv) |
Simple RNN | 0.5873 | 0.3168 | 0.2046 |
log(rv) |
LSTM | 0.5795 | 0.3141 | 0.1917 |
rv |
HAR-OLS | 0.6641 | 0.3262 | 0.1861 |
rv |
HAR-FNN | 0.6783 | 0.3304 | 0.1903 |
rv |
FNN | 0.7130 | 0.3416 | 0.2106 |
rv |
Simple RNN | 0.7154 | 0.3507 | 0.2259 |
rv |
LSTM | 0.7005 | 0.3288 | 0.1988 |
The nonpositive_rv_forecasts column of performance records whether a level-target model produced nonpositive variance forecasts. In this run, no model did—every count is zero, so all ten rows above have well-defined entries. The check is not decorative, however. Both the log-scale losses and QLIKE require a strictly positive variance forecast: the first takes \log \widehat{RV}_t, the second the ratio RV_t/\widehat{RV}_t. A common expedient is to clip nonpositive forecasts to a small constant such as 10^{-12} and score them anyway, but the resulting number is then governed by the arbitrary constant rather than by the forecast—clipping at 10^{-12} turns a handful of inadmissible forecasts into a QLIKE in the billions, while clipping at 10^{-6} would produce a different figure entirely. We therefore treat these losses as undefined for any model that submits a nonpositive variance. RV-scale squared and absolute errors remain well defined and are retained in performance, along with the violation count; the displayed table reports RMSE but omits mean absolute error (MAE) and the all-zero count column. A model that cannot deliver an admissible variance forecast has failed a precondition of variance-scale evaluation; that is a more informative verdict than a large but arbitrary score.
The failure mode is not hypothetical. A level-target network with a linear output layer has nothing that keeps \widehat{RV}_t positive. In an earlier run of this design on a slightly different sample, the level-target Simple RNN returned eight nonpositive forecasts. In the current run it stays positive but ranks last among the level-target models on all three displayed losses. Its validation-selected two-epoch fit and small learning rate make insufficient optimization one plausible explanation, but the reported evidence does not distinguish that hypothesis from seed variation, model misspecification, or the tuning rule itself. A positivity-enforcing output such as the softplus head used in the distributional networks chapter makes nonpositive forecasts impossible by construction. It removes the inadmissibility count as a diagnostic, so practitioners should instead monitor the minimum forecast, the frequency of near-zero forecasts, individual QLIKE contributions, and variance calibration.
log(rv) (blue) versus fit on rv (orange); bars are annotated with their values; the vertical axis is linear and starts at zero, so visual differences are proportional to actual differences. Had a level-target model produced nonpositive variance forecasts, its bar would be replaced by an “undefined” marker; in this run all ten bars are defined.
Figure 8.3 makes the pipeline comparison visible. On the common log-realized-variance scale, the log-target full FNN is the best-performing model, but its RMSE advantage over HAR-OLS is 0.0036 log-RV units—about 0.6% of the baseline RMSE. The log-target HAR-FNN and HAR-OLS are identical to four decimals (0.5771 against 0.5771; the network is lower by 0.00006), the LSTM lands just behind them (0.5795), and the Simple RNN is last in the panel (0.5873). That is the main benchmark lesson: the neural networks are not automatically better just because they are more flexible. The level-target pipelines lose accuracy on the log-realized-variance scale throughout—by roughly 0.09 to 0.14 log-RV units for the same architecture—and the level-target Simple RNN is the weakest model in the figure. Because lag normalization also differs across target pipelines, these gaps do not identify the effect of the target transformation alone.
Two notes help interpret these rankings. First, the yardstick for whether a ranking is informative is the standard error of the mean loss differential, not the dispersion of the individual per-period differentials—a distinction worth being explicit about, because the two differ by a large factor here. Let \delta_t denote the per-period difference in squared forecast errors between two models. A Diebold and Mariano (1995) comparison contrasts \bar \delta = N_{\text{test}}^{-1}\sum_t \delta_t with \widehat{\text{se}}(\bar \delta) = (\hat\omega^2/N_{\text{test}})^{1/2}, where \hat\omega^2 is a heteroskedasticity- and autocorrelation-consistent (HAC) estimate of the long-run variance of \{\delta_t\}; the HAC correction is needed because loss differentials from volatility forecasts are serially correlated (Exercise 4.4 in the predictive-distributions chapter works through exactly this construction for score differentials). Up to that correction, the relevant band is narrower than the per-period standard deviation of \delta_t by a factor \sqrt{N_{\text{test}}} = \sqrt{499} \approx 22. Treating “smaller than one per-period standard deviation” as a tie would therefore declare essentially every comparison in the table a tie, which is the opposite mistake. We do not run the test here, so nothing in this chapter establishes which differences are statistically distinguishable or economically irrelevant. The 0.0036 log-RV-unit FNN-versus-HAR-OLS gap, and the 0.00006 HAR-FNN-versus-HAR-OLS gap, are numerically small relative to the baseline RMSE; without a loss-differential test or economic-value calculation, they indicate weak practical separation rather than established equivalence.
Second, the QLIKE column in the table above is the Patton (2011) quasi-likelihood loss for volatility, the standard comparison loss in the high-frequency volatility literature. It is robust in Patton’s sense that rankings by expected loss are preserved when the true variance is replaced by a conditionally unbiased proxy such as realized variance; sampling variation can still reverse the ordering of finite-sample average losses.
The two HAR rows hold the target and predictor information fixed while comparing a linear OLS procedure with a tuned nonlinear neural-network procedure. In the log-target panel, the nonlinear procedure leaves log-RV RMSE unchanged to four decimals, lowers RV-scale RMSE slightly, and produces a clearer QLIKE gain (0.1836 against 0.1957). In the level-target panel it performs worse: HAR-FNN loses to HAR-OLS on all three columns. In this sample and design, the selected HAR-FNN adds little incremental predictive accuracy beyond what the HAR components already deliver. Because the fitting, scaling, and selection procedures also differ, the comparison should not be interpreted as isolating nonlinearity alone.
QLIKE and log-RV RMSE also need not pick the same winner. QLIKE is scale-invariant, that is, homogeneous of degree zero: it depends on forecast and target only through the ratio RV_t/\widehat{RV}_t, so rescaling both by a common positive constant leaves the loss unchanged. It measures whether \widehat{RV}_t tracks RV_t accurately in ratio terms, and it penalizes underprediction more heavily than overprediction through the ratio RV_t/\widehat{RV}_t. A model that calibrates the relative size of volatility spikes well can win on QLIKE while losing on log-RV RMSE—which is why the HAR-FNN fit on \log(RV) attains the lowest QLIKE in the table while the full FNN fit on \log(RV) attains the lowest log-RV RMSE, and why the level-target HAR-OLS, mediocre on log-RV RMSE, sits within 0.003 of the best QLIKE.
With dependent financial data, an architecture comparison is interpretable only after the forecast target, lag construction, scaling, validation block, and tuning protocol have been fixed. The numerical ranking of LSTM, HAR-OLS, and the FNN is secondary to that design discipline.
8.7 Summary
Key Takeaways
- Empirical architecture comparisons must hold the forecast target, available predictor information, preprocessing and tuning budget fixed.
- The log-target pipeline performs best on log-scale accuracy and RV-scale RMSE here, while level-target HAR-OLS remains competitive under QLIKE.
- With common target and predictor information, the tuned nonlinear HAR-FNN changes accuracy only marginally relative to HAR-OLS in this design.
- These results concern predictive performance in one chronological evaluation and do not establish structural economic mechanisms.
Common Pitfalls
- Fit scalers on the training block and keep validation and test observations in chronological order.
- Construct HAR averages and input windows using only observations available before the target day.
- Check that level forecasts are positive before evaluating them with QLIKE.
- Do not attribute the log-pipeline comparison to target scaling alone, since its input normalization also changes.
8.8 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
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.↩︎