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

2 Cross-Validation

2.1 Overview

When building a machine learning model, we face several decisions that affect performance: which algorithm to use, and what values to choose for the tuning constants that govern a model’s flexibility—for instance, the degree of a polynomial regression. Here it is worth fixing terminology: ordinary parameters (such as regression coefficients \beta) minimize the in-sample objective for a given set of tuning constants, whereas hyperparameters (such as the polynomial degree d of this chapter’s running example or the architecture of a neural network) determine the objective or fitting procedure and are chosen by a separate criterion. That criterion is often out of sample, as in this chapter, although in-sample penalized criteria such as the Akaike or Bayesian information criterion can also be used. The shared challenge here is that we cannot evaluate predictive performance on the same data used for fitting—doing so produces optimistic estimates that do not generalize to new observations.

This chapter introduces cross-validation as the principled answer to that challenge. Cross-validation simulates the experience of predicting on unseen data by systematically holding out portions of the training set. The idea dates to Stone (1974) and Geisser (1975), who proposed leave-one-out (holding out one observation at a time) and predictive-sample-reuse procedures for model assessment well before “machine learning” was a separate field. For econometricians, the issues are familiar: overfitting, specification search, and the bias introduced by repeated use of the same data arise in classical model selection just as in machine learning. The econometric analogue of the problem is pretest and specification-search bias; the machine-learning community’s answer to it is nested cross-validation—carrying out model selection on one layer of held-out data while keeping a separate outer layer untouched for the final performance assessment (Varma and Simon 2006; Cawley and Talbot 2010)—an idea we develop at the end of this chapter.

2.2 Roadmap

We illustrate overfitting with a polynomial regression example to motivate why the training mean squared error (MSE) is not a reliable guide to predictive performance.
We introduce an honest workflow—separating training, validation, and test sets—and explain when each set should be touched.
We define K-fold cross-validation—partitioning the sample into K groups and rotating which group is held out—discuss common choices of K, and explain why letting information from held-out observations slip into the training step typically makes estimated predictive performance look better than it is.
We extend to time-series cross-validation—designs in which the training sample either grows over time or moves forward at a fixed length—and explain why standard K-fold is often invalid for dependent data, and when it can survive.
We close with a practical example of hyperparameter tuning, in which cross-validation selects among candidate model configurations, and then develop nested cross-validation as a two-loop procedure for assessing the entire tuning-and-fitting workflow.

2.3 A Simple Illustration of Overfitting

Our running example throughout this chapter is polynomial regression. The example is deliberately familiar: it isolates overfitting and model selection without introducing new estimation machinery. Richer model classes in later chapters, such as neural networks and tree-based methods, add more choices but do not change the validation logic.

Model and flexibility. Let X_i denote the scalar predictor and Y_i the outcome of unit i, and write x_i and y_i for their observed realizations, following the convention fixed in the information theory chapter. Consider a regression model whose errors are independent and identically distributed (i.i.d.) Gaussian conditional on the predictor:

Y_i = f(X_i) + \varepsilon_i, \qquad \varepsilon_i \mid X_i \stackrel{\text{iid}}{\sim} \mathcal{N}(0,\sigma^2),

so that f(x) = \mathbb{E}[Y_i \mid X_i = x] is the conditional mean function at a generic predictor value x. We approximate f by fitting a polynomial of degree d. At predictor value x, the resulting prediction is

\hat y = \hat f(x) = \hat\beta_0 + \hat\beta_1 x + \cdots + \hat\beta_d x^d.

We estimate the coefficients by ordinary least squares (OLS) on an intercept and the constructed regressors x_i,x_i^2,\dots,x_i^d. Each x_i^j is the jth power of the same observed predictor for unit i, not a separate variable. The degree d is a hyperparameter: increasing it moves the model from rigid to increasingly flexible.

A word on vocabulary. What econometrics calls estimating a model, the machine-learning literature calls training or fitting it: choosing parameter values by optimizing an objective on a given sample. The three verbs describe the same operation, and we use them interchangeably throughout the book; the sample used for this step is the training set (or estimation sample), and the objective optimized on it is the training objective introduced next. Where the distinction between the estimator, the rule that maps a sample to \hat\beta, and the numerical algorithm that computes it matters, as in the neural-network and optimization chapters, the text says so explicitly. What differs across the two fields is emphasis rather than mechanics: an econometrician estimating \beta usually wants to interpret the estimate and attach a standard error to it, whereas training is judged by how well the fitted function \hat f predicts, which is the perspective of this chapter.

Training objective. Fix a loss function L(y,\hat y) that measures the cost of predicting \hat y when the outcome is y. In this example, we use squared loss, L(y,\hat y)=(y-\hat y)^2. For a fitted model \hat f and a sample (x_i,y_i)_{i=1}^N, the empirical risk is the average loss

\frac{1}{N}\sum_{i=1}^{N}L\big(y_i,\hat{f}(x_i)\big),

which under squared loss is the MSE on that sample. Recall from the information theory chapter that maximizing a likelihood is equivalent to minimizing the sample average negative log-likelihood. Under Gaussian errors with fixed variance, that objective is, up to an additive constant and a positive scale factor, the training MSE. Thus least squares and maximum likelihood select the same polynomial coefficients.

Evaluating the empirical risk on the data used for fitting gives the training MSE. Because \hat f is chosen to make this quantity small, a sufficiently flexible polynomial can reduce it by following sample noise rather than the conditional mean.

Generalization target. The generalization (prediction) risk is the expected loss on a fresh draw (X,Y) from the same distribution, independent of the data used to fit \hat f:

R(\hat{f})=\mathbb{E}_{(X,Y)}\big[L\big(Y,\hat{f}(X)\big)\,\big|\,\hat{f}\big].

The conditioning bar holds the fitted rule \hat f fixed; the expectation runs only over the new observation. The training MSE is optimistic for this risk because the same observations helped determine \hat f. By contrast, the test MSE evaluates the fitted rule on observations excluded from estimation and therefore provides an out-of-sample estimate of R(\hat f), with accuracy that depends on the size of the test sample. Cross-validation pursues the same out-of-sample target without committing a large part of the available sample to one validation split. As we explain later, however, K-fold cross-validation targets the risk of a fitting procedure trained on the fold-induced sample size rather than exactly R(\hat f) for the final full-sample fit.

Simulation design. For Figure 2.1, we simulate N=25 training observations and, separately, 200 test observations from the same data-generating process. The test observations are independent fresh draws, not a subset removed from the training sample. We draw the predictor on [-1,1] and interpret y as an economic outcome whose conditional mean varies smoothly but nonlinearly with x, as in an Engel curve or a response of output growth to financial conditions. The left panel compares a rigid linear fit with a flexible degree-12 fit; the right panel plots training and test MSE across polynomial degrees.

Show the code

import numpy as np
import matplotlib.pyplot as plt

# Set random seed for reproducibility
np.random.seed(42)

# Generate data from a nonlinear signal plus Gaussian noise
n_train = 25
n_test = 200
x_train = np.sort(np.random.uniform(-1, 1, n_train))
x_test = np.sort(np.random.uniform(-1, 1, n_test))

def f(x):
    return np.sin(np.pi * x)

y_train = f(x_train) + np.random.normal(0, 0.25, n_train)
y_test = f(x_test) + np.random.normal(0, 0.25, n_test)
x_grid = np.linspace(-1, 1, 400)
y_true = f(x_grid)

degrees = range(1, 16)
train_mse = []
test_mse = []

for d in degrees:
    coefs = np.polyfit(x_train, y_train, d)
    yhat_train = np.polyval(coefs, x_train)
    yhat_test = np.polyval(coefs, x_test)
    train_mse.append(np.mean((y_train - yhat_train) ** 2))
    test_mse.append(np.mean((y_test - yhat_test) ** 2))

# Two representative fits
coefs_low = np.polyfit(x_train, y_train, 1)
coefs_high = np.polyfit(x_train, y_train, 12)
yhat_low = np.polyval(coefs_low, x_grid)
yhat_high = np.polyval(coefs_high, x_grid)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))

# Left panel: fitted curves
axes[0].scatter(x_train, y_train, color='black', s=35, alpha=0.8, label='Training data')
axes[0].plot(x_grid, y_true, color='gray', linestyle='--', linewidth=2, label='True mean')
axes[0].plot(x_grid, yhat_low, color='C0', linewidth=2, label='Degree 1 fit')
axes[0].plot(x_grid, yhat_high, color='C3', linewidth=2, label='Degree 12 fit')
axes[0].set_title('Fits on the Training Sample')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].legend(frameon=False)
axes[0].grid(True, alpha=0.3)

# Right panel: training vs test MSE
axes[1].plot(list(degrees), train_mse, 'o-', color='C0', label='Training MSE')
axes[1].plot(list(degrees), test_mse, 'o-', color='C3', label='Test MSE')
axes[1].axvline(12, color='gray', linestyle='--', alpha=0.7)
axes[1].set_title('Error vs Model Complexity')
axes[1].set_xlabel('Polynomial degree')
axes[1].set_ylabel('Mean squared error')
axes[1].legend(frameon=False)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Figure 2.1: Overfitting in polynomial regression. Left panel: the N=25 training observations (black dots), the true conditional mean (dashed gray), and two least-squares fits—degree 1 (blue) and degree 12 (red). Right panel: training MSE (blue) and test MSE (red) as functions of the polynomial degree, with the dashed vertical line marking degree 12.

The left panel of Figure 2.1 shows the fitted mean functions. The degree-12 polynomial follows local noise in the training sample, while the linear model cannot reproduce the nonlinear signal. The right panel shows the resulting overfitting pattern: as the degree rises, the training MSE falls monotonically, but the test MSE eventually rises. The independent test sample reveals this divergence, but repeatedly using the test MSE to choose d would turn the test set into a model-selection device rather than a final evaluation sample. This tension motivates the train-test workflow in the next section.

2.4 The Problem with Simple Train-Test Splits

A common first approach is to split the data into a training set and one held-out set:

Train different models/configurations on the training set
Evaluate their performance on the held-out set
Choose the best-performing option

Because the held-out set is used to choose among candidates, it functions as a validation set, regardless of whether we initially called it a test set. It supports model selection but no longer provides an untouched final assessment. By contrast, if one fitting procedure had been fully specified in advance and the held-out set were used exactly once after all decisions, it would be a legitimate test set.

Specifically, relying on a single split has three major drawbacks:

High variance: The estimated performance varies with the random assignment of observations to the two sets, so a single split gives a noisy ranking.
Test set contamination: If we use the test set to make multiple decisions (choose hyperparameters, select the model specification), we are effectively “training” on the test set, leading to overly optimistic performance estimates.
Limited data usage: We are not using all available data for training.

The K-fold design introduced below addresses the first and third drawbacks: by averaging over K splits, it reduces the dependence of the estimate on any one arbitrary split, and it lets every observation serve in validation exactly once. Cross-validation does not by itself prevent contamination from repeated model comparison.

Reliable Workflow: Validation and Test Sets

The clean workflow is:

Use the training data to fit candidate models.
Compare those candidates and tune hyperparameters on a validation set: data held out from fitting and used only to choose among models. Reuse should nevertheless be limited, because extensive adaptive comparison can overfit the validation set through the selection process even when no model is fitted directly on its observations motivating nested cross-validation later in this chapter
Obtain a final performance assessment either from a test set kept untouched until the very end or from the held-out losses of an outer layer of nested cross-validation, with all model selection and preprocessing repeated inside each outer training sample.

2.5 K-Fold Cross-Validation

We now need a systematic way to rotate which observations are held out. Cross-validation serves two related purposes:

Two Uses of Cross-Validation

Performance estimation: Estimate how well a model will perform on data not used in training.
Model selection: Choose between different hyperparameters, model specifications, or algorithms based on their cross-validated performance.

Typical choices made this way include the polynomial degree of this chapter’s running example and the many structural and estimation settings of neural networks and tree-based methods in later chapters.

The most common type of cross-validation is K-fold cross-validation. In the form given here it is designed for settings in which the observations can be treated as exchangeable—in particular for independently sampled cross-sectional data without cluster or spatial dependence; Section 2.6 develops the variants appropriate for time series. The procedure is as follows:

Shuffle the dataset randomly.
Split the dataset into K equal-sized groups (or “folds”).
For each fold:
1. Use the fold as a validation set.
2. Use the remaining K-1 folds as a training set.
3. Train the model on the training set and evaluate it on the validation set.
Average the held-out losses from the K folds to obtain a single loss estimate.

Formally, write the dataset as \{(x_i,y_i)\}_{i=1}^N, where i indexes observations, x_i denotes all predictor information supplied to the model for observation i, and y_i is the corresponding observed outcome. Thus x_i is the scalar predictor in the polynomial example of Figure 2.1, but it may be a vector of predictors more generally. Let D_k be the index set of the observations in fold k, let |D_k| be its size, and let \hat{f}^{-k} denote the model refit from scratch on all data except fold k—that is, trained only on the K-1 remaining folds, including any data-dependent preprocessing re-estimated on those folds alone. Consequently, \hat f^{-k}(x_i) is the prediction for held-out observation i\in D_k, obtained without using any observation in D_k during fitting or preprocessing. The pooled K-fold cross-validation (CV) estimate of the risk is

\widehat{\mathrm{CV}}_K =\frac{1}{N}\sum_{k=1}^{K}\sum_{i\in D_k}L\big(y_i,\hat{f}^{-k}(x_i)\big).

This form weights every held-out observation equally. When the folds have equal size, it is equivalent to the familiar average of the fold means,

\widehat{\mathrm{CV}}_K =\frac{1}{K}\sum_{k=1}^{K}\frac{1}{|D_k|}\sum_{i\in D_k}L\big(y_i,\hat{f}^{-k}(x_i)\big).

When K does not divide N, the folds differ slightly in size and an unweighted average of their means is not exactly the pooled per-observation average; the pooled form above uses weights |D_k|/N on the fold means.

What does \widehat{\mathrm{CV}}_K estimate? With equal folds, each fold model is trained on N(K-1)/K observations. Under the independent sampling assumed here, \widehat{\mathrm{CV}}_K estimates the average out-of-sample risk of the fitting procedure at that training-sample size, where the average is over hypothetical training samples and fresh test observations. It does not estimate the realized risk R(\hat f) of the particular model refitted on the observed full sample because the full sample is larger than any fold model’s training sample.

Optional: A special case for OLS

For OLS in a homoskedastic Gaussian linear model under squared loss, Bates, Hastie, and Tibshirani (2024) show that, conditional on the design matrix, the CV error estimate is independent of the prediction risk of the full-sample fit. This special-case result illustrates the distinction between a fitting procedure’s average performance and the risk of one fitted model.

Common Choices for K

A common choice for K is 5 or 10; the extreme case K=N, in which each validation fold contains a single observation, is called leave-one-out cross-validation. A higher K means each fold model is trained on a sample closer in size to the full sample. When prediction risk decreases with training-sample size, this reduces the training-size bias, but it also raises the computational cost. How the variance of \widehat{\mathrm{CV}}_K moves with K has no universal answer: it depends on the stability of the fitting algorithm and on the correlation between fold estimates, and either direction can occur (Arlot and Celisse 2010).

The number of folds K is a design choice, not a hyperparameter of the fitted model: K determines the validation criterion but is not selected by minimizing that same criterion.

Leakage Inside Cross-Validation

Cross-validation only works if each fold is treated like genuinely unseen data. The general principle is:

In each CV split, every preprocessing or transformation rule used to build predictors or outcomes must be a function of that split’s training observations only. Apply the resulting rule unchanged to the validation observations.

Concretely, the following steps must be re-estimated using each split’s training observations rather than once on the full sample:

scaling and standardization
imputation of missing values
feature selection or dimensionality reduction
outcome transformations chosen from the data

Full-sample preprocessing lets validation observations influence a step of the fitting procedure. Whether this changes predictions or the reported loss depends on the fitting rule: refitted unpenalized OLS with an intercept is invariant to centering and positive rescaling, whereas adding a fixed multiple of the squared standardized slope to the fitting loss can make the scale matter. When predictions change, the resulting change in squared validation loss can have either sign.

Question for Reflection

Suppose lagged excess returns used as predictors are winsorized before K-fold cross-validation. Analyst A clips these predictor values at thresholds fixed ex ante at -5\% and +5\%. Analyst B clips them at the sample’s 1st and 99th percentiles, computed once from the full dataset. Which procedure leaks validation information, and how should the data-dependent procedure be implemented instead?

2.6 Time Series Cross-Validation

Before discussing time-series validation designs, we need one concept that recurs throughout the book.

Definition: Forecast-Origin Information Set

For a forecast made at time t, the forecast-origin information set \mathcal{F}_t is the \sigma-algebra generated by all variables—outcomes, predictors, metadata—in the form in which they were actually available at real-world time t. For data subject to revision, this means the vintage published by time t: a provisional release available at t belongs to \mathcal{F}_t, while its later revised value does not. A validation design is time-aware if, at each forecast origin t, every predictor entering the forecast and every observation used for estimation is measurable with respect to \mathcal{F}_t. The validation outcome is subsequently realized and is not \mathcal{F}_t-measurable, since it is precisely what the forecast is evaluated against. In a time-series setting, leakage violates this measurability condition: an object on the right-hand side of the forecasting regression depends on information that was not available at t. This is the temporal version of the general leakage principle above. In both cases, a training-side rule depends on held-out information—validation observations in the K-fold setting and post-origin data here. The same concept reappears in the chapters on predictive-distribution evaluation, recurrent networks, hyperparameter optimization, and conformal prediction.

Standard K-fold cross-validation is generally not suitable for time-series data: the argument that each held-out fold mimics genuinely unseen data relies on the observations being exchangeable—as they are under i.i.d. sampling—and temporal dependence breaks exchangeability. Randomly shuffling and splitting time-series data gives some fitted models access to observations after the validation date, so the split no longer reproduces a real-time estimation sample. That mismatch is not automatically biasing: for a stationary, purely autoregressive model that nests the true lag structure, so its errors are serially uncorrelated, random cross-validation can remain asymptotically valid (Bergmeir, Hyndman, and Koo 2018). It becomes optimistic when a misspecified or flexible model can exploit dependence between the validation block and post-validation observations used for training, and feature timing or overlapping outcomes can create leakage even when the model itself is correctly specified. Time-aware designs are therefore the prudent default.

Time-series applications therefore generally need validation designs tailored to the prediction task. Expanding- and sliding-window designs are natural for strict real-time forecasting; random or two-sided designs require assumptions that justify training on observations after a validation date. More generally, the validation split should mirror the dependence structure of the predicted outcome; Roberts et al. (2017) give a useful taxonomy for temporal, spatial, hierarchical, and related dependence settings.

Econometric Warning

In forecasting applications, leakage can arise even without explicit shuffling. Examples include:

predictors that are revised ex post
features published with delay
overlapping forecast outcomes
preprocessing performed using information from the full sample

So “time-aware” cross-validation must reflect the actual information set available at the forecast origin.

Serial dependence near fold boundaries is a related but distinct problem. Autocorrelation by itself is not leakage in the sense defined above: it puts no unavailable information into the forecast, and a validation block that lies strictly after the training window remains genuinely unseen even when it is dependent on the training history. What dependence changes is what the validation loss measures. Observations just after the training window are partly predictable from recently shared shocks, so a design that always validates immediately next to the training data estimates short-horizon performance, and the fold-average losses are correlated draws rather than independent ones. The validation design must therefore match the intended deployment horizon. Training on observations after the validation block can make the estimate genuinely optimistic when a misspecified or flexible model can exploit dependence between that block and the post-validation observations used for training; the purely autoregressive exception described above remains important.

A separate source of genuine optimism is overlap in labels, the machine-learning term for the outcomes being predicted. If an outcome is constructed from a window of future periods, training observations whose label windows reach into the validation block use information from it. In financial applications one common response is purged cross-validation (Lopez de Prado 2018), which removes training observations whose labels overlap with the validation period, sometimes combined with an embargo period that drops a short buffer after the validation block before training resumes. The exact buffer length is an econometric design choice: it should reflect the forecast horizon, label construction, and persistence in the data.

Expanding Window (Recursive Estimation)

In this approach, the training set grows with each split, and the validation set is always a block of data that comes after the training data.

Start with a small training set.
The next block of data is the validation set.
In the next iteration, the previous validation set is added to the training set, and the subsequent block becomes the new validation set.

This retraining scheme mimics a real-world forecasting workflow in which the model is periodically re-estimated as new data become available.

Sliding Window (Rolling Estimation)

This method uses a fixed-size training set that “slides” through time.

Select a block of data for training.
The next block of data is the validation set.
In the next iteration, both the training and validation windows slide forward by a specified step.

A sliding window can be preferred when the data-generating process (DGP) exhibits slow non-stationarity or structural change, so that older observations are less informative about the current conditional distribution—provided the relevance gain outweighs the usual loss of estimation precision from using fewer observations.

In the econometric forecast-evaluation literature, these are pseudo-out-of-sample (pseudo-OOS) schemes: an expanding training window uses recursive estimation, whereas a fixed-length sliding window uses rolling estimation (West 1996; Clark and McCracken 2013). In both schemes the forecast origin advances through time and the model is re-estimated using only the observations that the chosen window makes available at that origin.

Visualizing Time Series Cross-Validation

The first two panels of Figure 2.3 illustrate expanding- and sliding-window validation. The third panel switches to a two-sided K-fold layout to show where purge and embargo buffers bind when labels overlap.

Show the code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from sklearn.model_selection import TimeSeriesSplit

def plot_cv_indices(splits, n_samples, ax, k, show_legend=True):
    """Visualizes cross-validation splits with three colors."""
    colors = ['tab:blue', 'tab:red', 'lightgray']  # training, validation, unused
    
    for i, (train, test) in enumerate(splits):
        color_indices = np.full(n_samples, 2)  # Default to unused
        color_indices[train] = 0  # Training = blue
        color_indices[test] = 1   # Validation = red
        
        ax.scatter(range(n_samples), [i + 0.5] * n_samples,
                   c=[colors[j] for j in color_indices], marker='_', lw=10)

    ax.set_yticks(np.arange(k) + 0.5)
    ax.set_yticklabels([f"Split {i+1}" for i in range(k)])
    ax.set_ylim(k, -0.2)
    ax.set_xlim(0, n_samples)
    ax.set_xlabel("Data Index")
    ax.set_ylabel("CV Iteration")
    
    if show_legend:
        legend_elements = [
            Rectangle((0, 0), 1, 1, color='tab:blue', label='Training Set'),
            Rectangle((0, 0), 1, 1, color='tab:red', label='Validation Set'),
            Rectangle((0, 0), 1, 1, color='lightgray', label='Unused')
        ]
        ax.legend(handles=legend_elements, loc=(1.02, 0.5), title="Data Role")
    return ax

n_samples = 100
n_splits = 5

fig, axs = plt.subplots(3, 1, figsize=(8, 8.5), sharex=True)

# Expanding Window
tscv_expand = TimeSeriesSplit(n_splits=n_splits, test_size=10)
splits_expand = list(tscv_expand.split(np.arange(n_samples)))
plot_cv_indices(splits_expand, n_samples, axs[0], n_splits, show_legend=False)
axs[0].set_title("Recursive CV: Expanding Window")

# Sliding Window
def sliding_window_split(n_samples, n_splits=5, train_size=30, test_size=10, step=15):
    splits = []
    for i in range(n_splits):
        train_start = i * step
        train_end = train_start + train_size
        test_end = train_end + test_size
        if test_end > n_samples: break
        splits.append((np.arange(train_start, train_end), np.arange(train_end, test_end)))
    return splits

splits_slide = sliding_window_split(n_samples, n_splits=n_splits)
plot_cv_indices(splits_slide, n_samples, axs[1], len(splits_slide), show_legend=False)
axs[1].set_title("Rolling CV: Sliding Window")

# Two-sided K-fold splits with binding purge and embargo buffers
def purged_embargoed_kfold(n_samples, n_splits=5,
                           purge_size=5, embargo_size=5):
    splits = []
    all_indices = np.arange(n_samples)
    for test in np.array_split(all_indices, n_splits):
        test_start, test_end = test[0], test[-1] + 1
        purge = np.arange(max(0, test_start - purge_size), test_start)
        embargo = np.arange(test_end, min(n_samples, test_end + embargo_size))
        excluded = np.concatenate((purge, test, embargo))
        train = np.setdiff1d(all_indices, excluded)
        splits.append((train, test, purge, embargo))
    return splits

def plot_purged_embargoed_indices(splits, n_samples, ax):
    colors = {
        'training': 'tab:blue',
        'validation': 'tab:red',
        'unused': 'lightgray',
        'purged': 'purple',
        'embargoed': 'dimgray'
    }
    for i, (train, test, purge, embargo) in enumerate(splits):
        roles = np.full(n_samples, 'unused', dtype=object)
        roles[train] = 'training'
        roles[test] = 'validation'
        roles[purge] = 'purged'
        roles[embargo] = 'embargoed'
        ax.scatter(range(n_samples), [i + 0.5] * n_samples,
                   c=[colors[role] for role in roles], marker='_', lw=10)

    ax.set_yticks(np.arange(len(splits)) + 0.5)
    ax.set_yticklabels([f"Split {i+1}" for i in range(len(splits))])
    ax.set_ylim(len(splits), -0.2)
    ax.set_xlim(0, n_samples)
    ax.set_xlabel("Data Index")
    ax.set_ylabel("CV Iteration")
    legend_elements = [
        Rectangle((0, 0), 1, 1, color=colors['training'], label='Training Set'),
        Rectangle((0, 0), 1, 1, color=colors['validation'], label='Validation Set'),
        Rectangle((0, 0), 1, 1, color=colors['purged'], label='Purged'),
        Rectangle((0, 0), 1, 1, color=colors['embargoed'], label='Embargoed'),
        Rectangle((0, 0), 1, 1, color=colors['unused'], label='Unused')
    ]
    ax.legend(handles=legend_elements, loc=(1.02, 0.28), title="Data Role")
    return ax

splits_purged = purged_embargoed_kfold(n_samples, n_splits=n_splits)
plot_purged_embargoed_indices(splits_purged, n_samples, axs[2])
axs[2].set_title("Two-Sided Purged K-Fold with Embargo")

fig.tight_layout(rect=[0, 0, 0.85, 1]) # Adjust layout to make space for legend
plt.show()

Figure 2.3: Time-series cross-validation designs over 100 ordered observations. Each row represents one split; blue marks training data, red the validation block, and light gray observations unused in that split. Top: recursive pseudo-out-of-sample evaluation with an expanding training window. Middle: rolling pseudo-out-of-sample evaluation with a fixed-length training window of 30 observations. Bottom: two-sided purged K-fold for outcomes that look forward five periods, so overlap-based purging becomes a five-observation buffer (purple) before validation; a five-observation embargo (dark gray) follows validation. Existing buffers are excluded from training, but at the sample boundaries only one can bind.

The contrast between the first two panels of Figure 2.3 is the window-length choice: the expanding window never discards early observations, while the sliding window drops old observations as it moves. Sliding can be preferable when the gain from emphasizing recent conditions outweighs the usual loss of estimation precision from using fewer observations; that loss is not universal under heteroskedasticity or structural change. In both designs the validation block lies strictly after the training data, so no split trains on the future. The third panel addresses a different issue in a two-sided K-fold design. Purple observations are excluded because their multi-period outcomes would overlap the validation block; dark-gray observations are excluded by the post-validation embargo before training resumes. The returned training set contains the blue observations on both sides but excludes every buffer that exists for that split; the first split has no preceding purge and the last has no following embargo. Such a two-sided design is appropriate only when its maintained assumptions justify training on post-validation observations; it does not reproduce a strictly real-time forecast exercise. Purging and embargoing therefore modify a base validation design rather than replace the choice between recursive and rolling estimation.

2.7 Practical Example: Tuning a Neural Network with Cross-Validation

We now show how to use cross-validation to select hyperparameters for a simple neural network. No knowledge of neural networks is needed here: treat the network as a flexible nonlinear regression model whose fitting procedure comes with two tuning constants—the hidden-layer architecture (which controls how rich the fitted function class is) and the learning rate (a step-size constant of the estimation algorithm). We study these objects properly in the chapter on feed-forward neural networks; here the network simply provides a realistic tuning problem, and we tune both hyperparameters using scikit-learn. As data we simulate a single-year cross-section of 1{,}000 firms: annual excess returns, in percent, generated from ten firm characteristics measured on deliberately heterogeneous scales, with an interaction and a quadratic term that a purely linear model would miss. Because the observations are independently simulated—with no cluster, spatial, or time-series dependence—the K-fold design of Section 2.5 is appropriate here. A real firm-year panel would instead require folds that respect firm-level clustering and time.

The architecture labels below report how many neurons each hidden layer contains. Thus, 50 denotes one hidden layer with 50 units, while 100-50 denotes two hidden layers with 100 and 50 units, in that order. These labels identify candidate architectures; width and depth both affect model capacity, so the labels do not place the candidates on a single ordered flexibility scale.

Standardizing the firm characteristics can make fitting easier because they are measured on very different scales; the neural-network chapter explains why. Scikit-learn’s Pipeline object combines these two ordered steps—standardization followed by model fitting—and treats them as one fitting procedure. Within each CV split, the pipeline estimates the means and standard deviations only from that split’s training observations, fits the network on the standardized training predictors, and applies those same scaling constants to the validation predictors. The validation observations therefore do not influence standardization.

We fix the fitting procedure explicitly at no more than 1,000 training epochs per CV split, where one epoch is one full pass through the training fold. In this example, scikit-learn’s max_iter counts epochs, while each weight update uses a mini-batch—a subset of the training fold. With 800 training observations per fold and the default mini-batch size of 200, a fit that reaches the epoch limit makes 4,000 gradient updates. The optimization chapter explains this distinction between epochs and updates.

Instead of suppressing convergence warnings, the code records how many of the five training runs reach 1,000 epochs for each configuration; these runs may stop before the algorithm’s convergence criterion is met. Within the manual grid-search loop, scikit-learn’s cross_validate function evaluates one fixed pipeline configuration: it refits the pipeline on the training observations in each of the five splits and returns the corresponding held-out fold results. The surrounding ParameterGrid loop repeats this evaluation for every configuration, after which the code selects the configuration with the smallest mean cross-validated MSE. Figure 2.4 summarizes the search: each point is one configuration in the grid, plotted against the learning rate in the left panel and against the architecture in the right panel, with the red star marking the selected configuration. Scikit-learn follows a higher-is-better scoring convention, so neg_mean_squared_error returns the negative of each fold’s MSE; the minus sign in the code converts the average back to ordinary MSE before comparison and plotting.

Show the code

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import cross_validate, ParameterGrid, KFold
from sklearn.exceptions import ConvergenceWarning
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import pandas as pd
import warnings
from IPython.display import display

# Set random seeds for reproducibility
np.random.seed(42)

# Simulated single-year cross-section: annual excess returns (%) on
# 10 firm characteristics with heterogeneous units
rng = np.random.default_rng(42)
n_obs, n_char = 1000, 10
X_standard = rng.standard_normal((n_obs, n_char))
char_scales = np.array([1.0, 50.0, 0.02, 10.0, 0.5,
                        100.0, 2.0, 0.1, 25.0, 5.0])
char_locations = np.array([0.0, 100.0, 0.05, 20.0, 1.0,
                           500.0, 0.0, 0.5, 50.0, 10.0])
X = char_locations + X_standard * char_scales
beta = np.array([1.5, -1.0, 0.8, 0.5, -0.5, 0.3, 0.0, 0.0, 0.0, 0.0])
y = 10.0 * (
    X_standard @ beta
    + 0.7 * X_standard[:, 0] * X_standard[:, 1]
    + 0.5 * X_standard[:, 2] ** 2
    + rng.standard_normal(n_obs) * 2.0
)

# Define concise labels: each number gives the units in one hidden layer
architecture_labels = {
    (50,): '50',
    (100,): '100',
    (50, 50): '50-50',
    (100, 50): '100-50',
}

# Define hyperparameter grid
param_grid = {
    'hidden_layer_sizes': list(architecture_labels),
    'learning_rate_init': [0.001, 0.01, 0.1],
}

# Create pipeline with scaling
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    # Switch off the coefficient penalty: this example tunes only the
    # architecture and learning rate.
    ('mlp', MLPRegressor(max_iter=1000, alpha=0.0, random_state=42))
])
cv = KFold(n_splits=5, shuffle=True, random_state=0)

# Perform grid search with cross-validation
results = []
for params in ParameterGrid(param_grid):
    # Update pipeline parameters
    pipeline_params = {f'mlp__{key}': value for key, value in params.items()}
    pipeline.set_params(**pipeline_params)
    
    # Perform the same reproducible 5-fold split for every configuration.
    # Capture convergence warnings and expose them through diagnostics below.
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always", ConvergenceWarning)
        cv_result = cross_validate(
            pipeline,
            X,
            y,
            cv=cv,
            scoring='neg_mean_squared_error',
            return_estimator=True,
        )

    neg_mse_values = cv_result['test_score']
    n_fits_at_limit = sum(
        estimator.named_steps['mlp'].n_iter_ >= 1000
        for estimator in cv_result['estimator']
    )
    n_convergence_warnings = sum(
        issubclass(warning.category, ConvergenceWarning)
        for warning in caught
    )
    
    results.append({
        'architecture': architecture_labels[params['hidden_layer_sizes']],
        'learning_rate': params['learning_rate_init'],
        'mean_cv_mse': -neg_mse_values.mean(),
        'fold_mse_sd': (-neg_mse_values).std(ddof=1),
        'n_fits_at_limit': n_fits_at_limit,
        'n_convergence_warnings': n_convergence_warnings,
    })

# Convert to DataFrame for easier analysis
results_df = pd.DataFrame(results)

# Find best configuration
best_idx = results_df['mean_cv_mse'].idxmin()
best_config = results_df.iloc[best_idx]

print("Best Configuration:")
print(f"Hidden units by layer: {best_config['architecture']}")
print(f"Learning Rate: {best_config['learning_rate']}")
print(f"Mean CV MSE: {best_config['mean_cv_mse']:.4f}")
print(f"Across-fold MSE standard deviation (ddof=1): {best_config['fold_mse_sd']:.4f}")
print(f"Fits reaching 1,000 epochs: {int(best_config['n_fits_at_limit'])}/5")
print(f"Captured convergence warnings: {int(best_config['n_convergence_warnings'])}")

# Visualize results
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Plot 1: Learning rate vs performance (all configurations)
for arch_val in results_df['architecture'].unique():
    subset = results_df[results_df['architecture'] == arch_val]
    ax1.scatter(subset['learning_rate'], subset['mean_cv_mse'], 
               label=arch_val, alpha=0.7, s=50)

ax1.set_xlabel('Learning Rate')
ax1.set_ylabel('CV MSE')
ax1.set_title('Performance vs Learning Rate')
ax1.set_xscale('log')
ax1.grid(True, alpha=0.3)
ax1.legend(title='Hidden units\nby layer')

# Plot 2: Architecture vs performance (all configurations)
# Create a mapping for architectures to x-positions
arch_names = results_df['architecture'].unique()
arch_mapping = {arch: i for i, arch in enumerate(arch_names)}
results_df['arch_pos'] = results_df['architecture'].map(arch_mapping)

# Add small random jitter to x-position for better visibility
np.random.seed(42)
jitter = np.random.normal(0, 0.05, len(results_df))

for lr_val in results_df['learning_rate'].unique():
    subset = results_df[results_df['learning_rate'] == lr_val]
    ax2.scatter(subset['arch_pos'] + jitter[subset.index], subset['mean_cv_mse'], 
               label=f'LR={lr_val}', alpha=0.7, s=50)

ax2.set_xlabel('Hidden units by layer')
ax2.set_ylabel('CV MSE')
ax2.set_title('Performance vs Hidden-Layer Architecture')
ax2.set_xticks(range(len(arch_names)))
ax2.set_xticklabels(arch_names)
ax2.grid(True, alpha=0.3)
ax2.legend(title='Learning Rate')

# Highlight best configuration
best_lr = best_config['learning_rate']
best_arch_pos = arch_mapping[best_config['architecture']]
best_mse = best_config['mean_cv_mse']

ax1.scatter([best_lr], [best_mse], color='red', s=100, marker='*',
           edgecolors='black', linewidth=1, label='Best Config', zorder=5)
ax2.scatter([best_arch_pos], [best_mse], color='red', s=100, marker='*',
           edgecolors='black', linewidth=1, label='Best Config', zorder=5)

# Update legends to include best config
ax1.legend(title='Hidden units\nby layer')
ax2.legend(title='Learning Rate')

plt.tight_layout()
plt.show()

Best Configuration:
Hidden units by layer: 50
Learning Rate: 0.001
Mean CV MSE: 448.4086
Across-fold MSE standard deviation (ddof=1): 38.5992
Fits reaching 1,000 epochs: 5/5
Captured convergence warnings: 5

Figure 2.4: Cross-validated MSE for every fixed-budget hyperparameter configuration in the grid search. Left panel: mean CV MSE against the learning rate (log scale), with marker color indicating the number of hidden units by layer. Right panel: the same configurations against those architecture labels, with color indicating the learning rate and horizontal jitter added for visibility. For example, 100-50 means two hidden layers with 100 and 50 units. In both panels the red star marks the configuration with the lowest mean CV MSE; the printed diagnostics report how many fold fits reached the 1,000-epoch cap.

The five configurations with the lowest mean CV MSE are reported separately in Table 2.1. Its first column uses exactly the same hidden-units-by-layer labels as the right panel of Figure 2.4. The Fits at epoch limit column counts how many of the five training runs reached 1,000 epochs; a value of 3 means that three runs used the full training budget.

Show the code

top_five = results_df.nsmallest(5, 'mean_cv_mse')[[
    'architecture', 'learning_rate', 'mean_cv_mse',
    'fold_mse_sd', 'n_fits_at_limit'
]].rename(columns={
    'architecture': 'Hidden units by layer',
    'learning_rate': 'Learning rate',
    'mean_cv_mse': 'Mean CV MSE',
    'fold_mse_sd': 'Across-fold MSE standard deviation',
    'n_fits_at_limit': 'Fits at epoch limit',
})
display(
    top_five.style
    .hide(axis='index')
    .format({
        'Learning rate': '{:.3g}',
        'Mean CV MSE': '{:.2f}',
        'Across-fold MSE standard deviation': '{:.2f}',
        'Fits at epoch limit': '{:.0f}',
    })
)

Table 2.1: Five configurations with the lowest mean cross-validated MSE. Hidden units by layer uses the same architecture labels as Figure 2.4; Fits at epoch limit counts, out of five, the training runs that reached 1,000 epochs.

Hidden units by layer	Learning rate	Mean CV MSE	Across-fold MSE standard deviation	Fits at epoch limit
50	0.001	448.41	38.60	5
100	0.001	454.40	36.12	5
100	0.1	612.38	52.91	0
50	0.1	637.94	72.81	0
50	0.01	643.34	97.46	4

How to Read the Diagnostics

The mean CV MSE ranks the candidate configurations. The Across-fold MSE standard deviation column—computed in the code as fold_mse_sd with ddof=1—describes how much the K validation-fold MSE values fluctuate, but it is not a standard error for the cross-validated risk: the fold-average losses are correlated because their training sets overlap (Bates, Hastie, and Tibshirani 2024).

Runs that reach the epoch limit may not have converged, so differences in CV MSE can reflect incomplete optimization as well as architecture or learning rate. Small differences in mean CV MSE should therefore not be read as precise rankings, especially when across-fold dispersion is large or different numbers of runs reach the limit.

The grid search above uses a single level of cross-validation. For each fixed configuration, cross_validate produces five returned values that the code converts to validation-fold MSEs. The ParameterGrid loop repeats that call for every configuration, and the subsequent idxmin step selects the configuration with the smallest mean CV MSE. This entire selection step plays the role of the inner tuning loop in the nested procedure developed next. Using the same splits for every configuration makes their loss estimates directly comparable, but the ranking remains noisy, as the dispersion and convergence diagnostics emphasize. Moreover, the mean CV MSE of the selected configuration is optimistically low, for exactly the reason formalized in Exercise 2.1 below: taking the smallest of many noisy validation-loss estimates favors lucky draws. Section 2.8 develops the separate outer layer needed to assess the tuning procedure without reusing its inner loss estimates.

Exhaustive grid search becomes costly as the number of hyperparameters grows. Random and sequential search methods provide alternatives, which we study in the chapter on advanced hyperparameter optimization. Regardless of how the candidates are searched, using the same validation-loss estimates for selection and assessment remains a separate problem.

2.8 Nested Cross-Validation

The practical example used cross-validation to answer a selection question: which candidate architecture and learning rate should we choose? Nested cross-validation adds an assessment layer for a second question: how well will the entire selection procedure predict on new data? Reporting the smallest inner-CV loss estimate as the answer reuses the same validation-loss estimates for both tasks. Because the minimum favors configurations that received unusually favorable validation draws, the reported estimate is generally optimistic even if each candidate’s CV loss estimate is unbiased before selection (Varma and Simon 2006; Cawley and Talbot 2010).

Figure 2.5 visualizes how nested cross-validation prevents that reuse. The left panel shows the outer loop that produces the final assessment losses. The right panel zooms into the first outer split and shows the separate inner loop used for tuning.

Show the code

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

outer_folds = 5
inner_folds = outer_folds - 1
row_height = 0.68

training_color = 'tab:blue'
outer_validation_color = 'tab:red'
inner_validation_color = 'tab:orange'

fig, (ax_outer, ax_inner) = plt.subplots(
    1, 2, figsize=(12, 5.2),
    gridspec_kw={'width_ratios': [1.0, 1.15]}
)

# Outer loop: rotate which fold is held out for final assessment.
for split in range(outer_folds):
    for fold in range(outer_folds):
        color = (
            outer_validation_color if fold == split else training_color
        )
        ax_outer.add_patch(Rectangle(
            (fold, split - row_height / 2), 0.96, row_height,
            facecolor=color, edgecolor='white'
        ))

# Outline the outer split expanded in the right panel.
ax_outer.add_patch(Rectangle(
    (-0.08, -row_height / 2 - 0.08),
    outer_folds + 0.04, row_height + 0.16,
    fill=False, edgecolor='black', linewidth=1.4
))
ax_outer.set_xlim(-0.15, outer_folds + 0.05)
ax_outer.set_ylim(outer_folds - 0.45, -0.65)
ax_outer.set_xticks(
    [fold + 0.48 for fold in range(outer_folds)],
    [fr'$D_{{{fold + 1}}}^{{\mathrm{{out}}}}$'
     for fold in range(outer_folds)]
)
ax_outer.set_yticks(
    range(outer_folds),
    [f'Outer split {split + 1}' for split in range(outer_folds)]
)
ax_outer.set_title('(a) Outer loop: assessment')
ax_outer.tick_params(length=0)

# Inner loop for outer split 1. The outer validation fold is displayed
# separately as a reminder that it never enters tuning.
outer_box_x = 0.0
inner_start_x = 1.35
for split in range(inner_folds):
    ax_inner.add_patch(Rectangle(
        (outer_box_x, split - row_height / 2), 0.82, row_height,
        facecolor=outer_validation_color, edgecolor='white', hatch='//'
    ))
    for fold in range(inner_folds):
        color = (
            inner_validation_color if fold == split else training_color
        )
        ax_inner.add_patch(Rectangle(
            (inner_start_x + fold, split - row_height / 2),
            0.96, row_height, facecolor=color, edgecolor='white'
        ))

ax_inner.axvline(1.08, color='0.5', linewidth=1, linestyle='--')
ax_inner.text(
    outer_box_x + 0.41, -0.57, 'untouched',
    ha='center', va='bottom', fontsize=9
)
ax_inner.set_xlim(-0.12, inner_start_x + inner_folds + 0.05)
ax_inner.set_ylim(inner_folds - 0.45, -0.65)
ax_inner.set_xticks(
    [outer_box_x + 0.41]
    + [inner_start_x + fold + 0.48 for fold in range(inner_folds)],
    [r'$D_1^{\mathrm{out}}$']
    + [fr'$D_{{{fold + 1}}}^{{\mathrm{{in}}}}$'
       for fold in range(inner_folds)]
)
ax_inner.set_yticks(
    range(inner_folds),
    [f'Inner split {split + 1}' for split in range(inner_folds)]
)
ax_inner.set_title('(b) Inner loop within outer split 1: tuning')
ax_inner.tick_params(length=0)

for ax in (ax_outer, ax_inner):
    for spine in ax.spines.values():
        spine.set_visible(False)

legend_handles = [
    Patch(facecolor=training_color, label='Training in current loop'),
    Patch(facecolor=inner_validation_color, label='Inner validation'),
    Patch(facecolor=outer_validation_color, hatch='//',
          label='Outer validation'),
]
fig.legend(
    handles=legend_handles, loc='lower center', ncol=3,
    frameon=False, bbox_to_anchor=(0.5, -0.01)
)
fig.tight_layout(rect=[0, 0.08, 1, 1])
plt.show()

Figure 2.5: Nested cross-validation with five outer folds and, within outer split 1, four inner folds. Left panel: each outer fold serves once as the outer validation fold (red), while the other folds form the outer training sample (blue). Right panel: the observations in the outer training sample are repartitioned into four inner folds, which rotate between inner training (blue) and inner validation (orange). The outer validation fold remains red and never enters tuning. After the inner loop selects a configuration, the model is refitted on the full outer training sample and evaluated once on the outer validation fold.

The right panel of Figure 2.5 makes the separation operational. The four inner folds rotate between training and validation, but D_1^{\mathrm{out}} never changes role. After the inner loss estimates select a configuration, denoted \hat\lambda_1 below, the procedure refits on the complete outer training sample and evaluates once on D_1^{\mathrm{out}}. Repeating the same sequence for the other four outer splits produces one reliable held-out loss for every observation.

Nested cross-validation separates the two questions. Let \Lambda denote the set of candidate hyperparameter configurations. In ordinary nested K-fold CV, the full sample is first partitioned into K_{\mathrm{out}} outer folds. Then:

Hold out outer fold D_k^{\mathrm{out}} and keep it completely untouched while making modeling decisions.
Using only the remaining outer training observations, run an inner cross-validation search over \lambda\in\Lambda. Every data-dependent step—including imputation, standardization, feature selection, and early stopping—must be repeated using each inner split’s training observations.
Select the configuration \hat\lambda_k with the best inner-CV criterion.
Refit the complete pipeline with \hat\lambda_k on all observations in the outer training sample.
Evaluate that fitted pipeline once on D_k^{\mathrm{out}}.
Repeat Steps 1–5 for every outer fold and pool the outer-fold losses.

For outer folds that partition the N observations, the resulting estimate is

\widehat R_{\mathrm{nested}} = \frac{1}{N} \sum_{k=1}^{K_{\mathrm{out}}} \sum_{i\in D_k^{\mathrm{out}}} L\!\left(y_i,\hat f_k^{\mathrm{tuned}}(x_i)\right),

where \hat f_k^{\mathrm{tuned}} denotes the entire rule produced from outer training sample k: inner-loop preprocessing and tuning, selection of \hat\lambda_k, and refitting with that selected configuration. The outer-fold outcome y_i has no role in constructing this rule. This separation is what makes the outer loss an assessment rather than another tuning criterion.

What does nested CV estimate? It estimates the expected performance of the complete tuning-and-fitting procedure when that procedure is trained on an outer training sample of size approximately N(K_{\mathrm{out}}-1)/K_{\mathrm{out}}. It does not estimate the risk of one prespecified configuration, nor exactly the risk of the final model fitted on all N observations. The selected configurations \hat\lambda_k may legitimately differ across outer folds: the object being evaluated is the rule choose a configuration from the available training data and then refit, not any one value of \lambda.

What happens after assessment? Once the nested estimate has been recorded, rerun the inner tuning procedure on the full development sample, choose a final configuration, and refit the pipeline on all available development observations for deployment. The nested outer loss estimate remains the performance estimate for that workflow. Do not replace it with the smaller inner-CV loss estimate obtained during the final full-sample search. If a genuinely untouched test set is available, it can instead provide the final assessment, but it must not become another dataset consulted repeatedly during tuning.

What changes with dependent data? Both layers must respect the dependence and information structure described in Section 2.6. For time series, each outer validation block must lie after its outer training window, and every inner split must preserve chronology within that outer training window. For clustered data, clusters rather than individual observations should remain intact in both layers. A time-aware outer loop cannot repair leakage created by random or otherwise invalid inner folds.

Keep the Three Roles Separate

Inner CV selects preprocessing choices, model classes, and hyperparameters.
Outer CV assesses the full procedure that performs that selection and then refits.
The final full-sample refit produces the model used for deployment; it does not produce a new reliable performance estimate.

2.9 Summary

Key Takeaways

Cross-validation compares fitting procedures using predictions for observations excluded from fitting.
With equal folds, K-fold CV targets average prediction risk for a procedure trained on N(K-1)/K observations.
Nested cross-validation uses inner folds for selection and outer folds to assess the complete tuning procedure.
Forecast validation should reproduce the timing and information available when predictions will be issued.

Common Pitfalls

Reporting the smallest tuning loss, or repeatedly inspecting the test set, makes the final assessment optimistic.
Fit preprocessing separately inside each training fold to prevent validation information from entering model fitting.
Across-fold MSE dispersion is not a standard error because the training samples overlap.
A training run that reaches the epoch limit has not necessarily converged.

2.10 Exercises

Exercise 2.1: Why Model Search Creates Optimistically Low Validation-Loss Estimates

Throughout this exercise, treat the fitted forecasting rules \hat f_A and \hat f_B as fixed. Their generalization risks are therefore fixed numbers, whereas their validation-loss estimates are random because the validation sample is random. Expectations and probabilities are over the validation sample, and in Part 5 also over the independent test sample introduced there.

For Parts 1–3, suppose the two fitted rules have the same finite generalization risk,

R(\hat{f}_A)=R(\hat{f}_B)=R_0.

Their validation-set risk estimates satisfy

\hat{R}_A = R_0+\eta_A, \qquad \hat{R}_B = R_0+\eta_B,

where \eta_A and \eta_B are the loss-estimation errors. Assume for these three parts that \eta_A and \eta_B have mean zero and that \mathbb{P}(\eta_A\neq\eta_B)>0, so the two estimates do not coincide with certainty. The errors need not be independent: candidate models evaluated on the same validation observations will generally have dependent loss-estimation errors.

The researcher selects the model with the smaller validation estimate, choosing A in a tie.

Show that \min(\hat{R}_A,\hat{R}_B) = R_0+\frac{\eta_A+\eta_B-|\eta_A-\eta_B|}{2}.
Deduce that \mathbb{E}\big[\min(\hat{R}_A,\hat{R}_B)\big] = R_0-\frac{1}{2}\mathbb{E}\big[|\eta_A-\eta_B|\big] < R_0.
Suppose in addition that (\eta_A,\eta_B) is jointly Gaussian with \operatorname{Var}(\eta_A)=\operatorname{Var}(\eta_B)=\sigma^2, \qquad \operatorname{Corr}(\eta_A,\eta_B)=c<1. Derive the exact expectation of the selected validation-loss estimate. Show that it equals \mathbb{E}\big[\min(\hat{R}_A,\hat{R}_B)\big] = R_0-\frac{\sigma\sqrt{1-c}}{\sqrt{\pi}}, and explain why stronger positive dependence between candidate loss estimates reduces the selection-induced optimism. Specialize the result to independent errors.
For Parts 4–5, drop the equal-risk restriction and write R_A=R(\hat f_A) and R_B=R(\hat f_B) for the fixed, finite generalization risks, labeling the models so that R_A\leq R_B. The random validation-loss estimates satisfy \mathbb E[\hat R_A]=R_A and \mathbb E[\hat R_B]=R_B; no Gaussian assumption is imposed. Compare \mathbb E[\min(\hat R_A,\hat R_B)] with the smaller generalization risk R_A. Prove your comparison and give a necessary and sufficient condition for the inequality to be strict. Express your condition as a probability statement about \hat R_A and \hat R_B.

Exam level. The exercise derives selection optimism for equal and unequal generalization risks and connects this to validation-based model selection.

Hint for Part 3

First compute \operatorname{Var}(\eta_A-\eta_B). Under the joint Gaussian assumption, the difference is Gaussian. Use the fact that if Z\sim\mathcal{N}(0,\tau^2), then

\mathbb{E}|Z|=\tau\sqrt{\frac{2}{\pi}}.

Solution

Part 1: Writing the Selected Validation-Loss Estimate

Use the elementary identity \min(a,b)=(a+b-|a-b|)/2, valid for all real a,b, with a=\hat{R}_A and b=\hat{R}_B:

\min(\hat{R}_A,\hat{R}_B) =\frac{\hat{R}_A+\hat{R}_B-|\hat{R}_A-\hat{R}_B|}{2}.

Substituting \hat{R}_A=R_0+\eta_A and \hat{R}_B=R_0+\eta_B gives

\min(\hat{R}_A,\hat{R}_B) =\frac{2R_0+\eta_A+\eta_B-|\eta_A-\eta_B|}{2} =R_0+\frac{\eta_A+\eta_B-|\eta_A-\eta_B|}{2}.

Part 2: Showing the Optimism Bias

Taking expectations,

\mathbb{E}\big[\min(\hat{R}_A,\hat{R}_B)\big] = R_0+\frac{\mathbb{E}[\eta_A]+\mathbb{E}[\eta_B]-\mathbb{E}[|\eta_A-\eta_B|]}{2}.

Since both errors have mean zero,

\mathbb{E}\big[\min(\hat{R}_A,\hat{R}_B)\big] = R_0-\frac{1}{2}\mathbb{E}[|\eta_A-\eta_B|].

Because |\eta_A-\eta_B|\ge 0 and, by the assumption \mathbb{P}(\eta_A\neq\eta_B)>0, it is strictly positive with positive probability,

\mathbb{E}\big[\min(\hat{R}_A,\hat{R}_B)\big] < R_0.

Part 3: Exact Gaussian Optimism

Under the additional assumption,

\eta_A-\eta_B \sim \mathcal{N}\big(0,2\sigma^2(1-c)\big).

If Z\sim\mathcal{N}(0,\tau^2), then

\mathbb{E}|Z|=\tau\sqrt{\frac{2}{\pi}}.

Here \tau=\sigma\sqrt{2(1-c)}, so

\mathbb{E}\big[|\eta_A-\eta_B|\big] =\sigma\sqrt{2(1-c)}\sqrt{\frac{2}{\pi}} =\frac{2\sigma\sqrt{1-c}}{\sqrt{\pi}}.

Substituting this into Part 2 gives

\mathbb{E}\big[\min(\hat{R}_A,\hat{R}_B)\big] = R_0-\frac{1}{2}\cdot \frac{2\sigma\sqrt{1-c}}{\sqrt{\pi}} = R_0-\frac{\sigma\sqrt{1-c}}{\sqrt{\pi}}.

As c increases, the two loss-estimation errors move together and their difference becomes less variable, so selecting the smaller estimate extracts less favorable noise. Setting c=0 gives the independent-error result R_0-\sigma/\sqrt{\pi}.

Part 4: Possibly Unequal Generalization Risks

The fitted rules remain fixed, so R_A and R_B are constants. For every validation sample, \min(\hat R_A,\hat R_B)\leq\hat R_A. Taking expectations gives

\mathbb E[\min(\hat R_A,\hat R_B)] \leq \mathbb E[\hat R_A]=R_A.

Since the models are labeled so that R_A\leq R_B, the right-hand side is the smaller generalization risk. To determine when the inequality is strict, write x^+:=\max\{x,0\} for the positive part and use the identity

\min(\hat R_A,\hat R_B) =\hat R_A-(\hat R_A-\hat R_B)^+,

which follows by considering \hat R_A\leq\hat R_B and \hat R_A>\hat R_B separately. Taking expectations gives

\mathbb E[\min(\hat R_A,\hat R_B)] =R_A-\mathbb E[(\hat R_A-\hat R_B)^+].

The nonnegative random variable (\hat R_A-\hat R_B)^+ has strictly positive expectation if and only if it is strictly positive with positive probability. Hence

\mathbb E[\min(\hat R_A,\hat R_B)]<R_A \quad\Longleftrightarrow\quad \mathbb P(\hat R_B<\hat R_A)>0.

When R_A<R_B, this is the event that validation strictly prefers the model with the larger generalization risk. If \hat R_A\leq\hat R_B almost surely, the expected minimum equals R_A. The same probability condition also applies when R_A=R_B.

Exercise 2.2: Feasible and Infeasible Time-Series Forecasts

Suppose the data-generating process is

Y_{t+1}=\beta Z_t+U_{t+1}, \qquad Z_t=\rho Z_{t-1}+\xi_t,

where |\rho|<1, \beta\neq 0, and \{\xi_t\} and \{U_t\} are mutually independent i.i.d. innovation sequences with mean zero and variances \sigma_\xi^2>0 and \sigma_U^2, respectively. Treat \beta and \rho as known, so the mean squared prediction error (MSPE) of a forecast \hat y_{t+1\mid t}, defined as \operatorname{MSPE}=\mathbb{E}\big[(Y_{t+1}-\hat y_{t+1\mid t})^2\big], is a population quantity net of parameter-estimation error.

Assume the predictor Z_t is published with a one-period delay, so at forecast origin t we observe Z_{t-1} but not Z_t. Consistently with the chapter’s definition, let the forecast-origin information set be generated by the available predictor and outcome histories,

\mathcal{F}_t =\sigma\big(\{Z_s:s\leq t-1\},\{Y_s:s\leq t\}\big).

Under the maintained independence assumptions, the past outcomes add no information about the new innovation \xi_t.

Consider the two forecasting rules

\hat{y}^{\,\text{oracle}}_{t+1\mid t}=\beta Z_t, \qquad \hat{y}^{\,\text{rt}}_{t+1\mid t}=\mathbb{E}[Y_{t+1}\mid \mathcal{F}_t].

Derive the feasible real-time forecast \hat{y}^{\,\text{rt}}_{t+1\mid t} as a function of Z_{t-1}. Then derive its forecast error and MSPE from the two equations of the data-generating process; do not leave the answer in conditional-expectation form.
Derive the oracle forecast error and its MSPE. Compute the exact MSPE gap between the feasible and oracle rules, and identify the innovation that generates this gap.
A researcher constructs the predictor–outcome pair for forecast origin t using Z_t, even though Z_t is published only at t+1, and evaluates the resulting oracle rule. State precisely which \mathcal{F}_t-measurability condition fails, explain why changing from random to chronological splits does not repair the mis-dated feature, and specify both the feature-timing correction and the split restriction needed for a feasible real-time validation exercise.

Exam level. The exercise formalizes the central econometric issue in time-series cross-validation: the validation design must match the real-time information set.

Hint for Part 1

Use linearity of conditional expectation and the independence of \xi_t and U_{t+1} from \mathcal{F}_t:

\mathbb{E}[Y_{t+1}\mid \mathcal{F}_t] = \beta\,\mathbb{E}[Z_t\mid \mathcal{F}_t].

Solution

Part 1: Feasible Forecast, Error, and MSPE

We have

\mathbb{E}[Y_{t+1}\mid \mathcal{F}_t] = \beta\,\mathbb{E}[Z_t\mid \mathcal{F}_t] = \beta\,\mathbb{E}[\rho Z_{t-1}+\xi_t\mid \mathcal{F}_t].

Since Z_{t-1} is \mathcal{F}_t-measurable and \mathbb{E}[\xi_t\mid \mathcal{F}_t]=0,

\hat{y}^{\,\text{rt}}_{t+1\mid t} = \beta\rho Z_{t-1}.

The real-time forecast error is

\varepsilon^{\text{rt}}_{t+1} = Y_{t+1}-\hat{y}^{\,\text{rt}}_{t+1\mid t} = \beta Z_t+U_{t+1}-\beta\rho Z_{t-1}.

Using Z_t=\rho Z_{t-1}+\xi_t, this becomes

\varepsilon^{\text{rt}}_{t+1} = \beta\xi_t+U_{t+1}.

Hence

\operatorname{MSPE}_{\text{rt}} = \operatorname{Var}(\beta\xi_t+U_{t+1}) = \beta^2\sigma_\xi^2+\sigma_U^2,

using independence of \xi_t and U_{t+1}.

Part 2: Oracle Forecast and the MSPE Gap

The oracle forecast error is

\varepsilon^{\text{oracle}}_{t+1} = Y_{t+1}-\hat{y}^{\,\text{oracle}}_{t+1\mid t} = \beta Z_t+U_{t+1}-\beta Z_t = U_{t+1}.

Therefore

\operatorname{MSPE}_{\text{oracle}}=\operatorname{Var}(U_{t+1})=\sigma_U^2.

The exact gap is therefore

\operatorname{MSPE}_{\text{rt}}-\operatorname{MSPE}_{\text{oracle}} =\beta^2\sigma_\xi^2>0,

where strict positivity follows from \beta\neq0 and \sigma_\xi^2>0. The gap is generated by \xi_t, the innovation in the unavailable current value Z_t. The oracle observes this innovation through Z_t, whereas the feasible forecaster cannot condition on it at origin t.

Part 3: Measurability Failure and Its Correction

The feature Z_t is not \mathcal F_t-measurable. Indeed, Z_t=\rho Z_{t-1}+\xi_t, where Z_{t-1} is \mathcal F_t-measurable but the nondegenerate innovation \xi_t is independent of \mathcal F_t. A protocol that places Z_t in the predictor–outcome pair for origin t therefore evaluates an infeasible forecasting setup and understates achievable real-time MSPE by \beta^2\sigma_\xi^2 in this data-generating process.

Switching from random to chronological splits changes which dated predictor–outcome pairs are used for estimation, but it cannot repair a pair whose feature timing is wrong. The feature-timing correction is to replace Z_t by the value actually observable at origin t, namely Z_{t-1} here. The split correction is to train only on observations available by each forecast origin and validate on a subsequent block, as in an expanding or sliding pseudo-out-of-sample design. Both corrections are necessary: feasible predictors do not excuse future training data, and chronological ordering does not excuse infeasible predictors.

2.11 References

Arlot, Sylvain, and Alain Celisse. 2010. “A Survey of Cross-Validation Procedures for Model Selection.” Statistics Surveys 4: 40–79. https://doi.org/10.1214/09-SS054.

Bates, Stephen, Trevor Hastie, and Robert Tibshirani. 2024. “Cross-Validation: What Does It Estimate and How Well Does It Do It?” Journal of the American Statistical Association 119 (546): 1434–45. https://doi.org/10.1080/01621459.2023.2197686.

Bergmeir, Christoph, Rob J. Hyndman, and Bonsoo Koo. 2018. “A Note on the Validity of Cross-Validation for Evaluating Autoregressive Time Series Prediction.” Computational Statistics & Data Analysis 120: 70–83. https://doi.org/10.1016/j.csda.2017.11.003.

Cawley, Gavin C., and Nicola L. C. Talbot. 2010. “On over-Fitting in Model Selection and Subsequent Selection Bias in Performance Evaluation.” Journal of Machine Learning Research 11: 2079–2107.

Clark, Todd E., and Michael W. McCracken. 2013. “Advances in Forecast Evaluation.” In Handbook of Economic Forecasting, edited by Graham Elliott and Allan Timmermann, 2:1107–201. Elsevier. https://doi.org/10.1016/B978-0-444-62731-5.00020-8.

Geisser, Seymour. 1975. “The Predictive Sample Reuse Method with Applications.” Journal of the American Statistical Association 70 (350): 320–28. https://doi.org/10.1080/01621459.1975.10479865.

Lopez de Prado, Marcos. 2018. Advances in Financial Machine Learning. Hoboken, NJ: Wiley.

Roberts, David R., Volker Bahn, Simone Ciuti, Mark S. Boyce, Jane Elith, Gurutzeta Guillera-Arroita, Severin Hauenstein, et al. 2017. “Cross-Validation Strategies for Data with Temporal, Spatial, Hierarchical, or Phylogenetic Structure.” Ecography 40 (8): 913–29. https://doi.org/10.1111/ecog.02881.

Stone, M. 1974. “Cross-Validatory Choice and Assessment of Statistical Predictions.” Journal of the Royal Statistical Society, Series B 36 (2): 111–47. https://doi.org/10.1111/j.2517-6161.1974.tb00994.x.

Varma, Sudhir, and Richard Simon. 2006. “Bias in Error Estimation When Using Cross-Validation for Model Selection.” BMC Bioinformatics 7: 91. https://doi.org/10.1186/1471-2105-7-91.

West, Kenneth D. 1996. “Asymptotic Inference about Predictive Ability.” Econometrica 64 (5): 1067–84. https://doi.org/10.2307/2171956.