You build a model to predict tomorrow’s air quality. It scores R² = 0.66 on your test set. You deploy it. It fails to beat a model that predicts “same as yesterday.” A hospital trains a sepsis-alert model. It validates at AUC 0.76. Clinicians start acting on the alerts; within a year, the model’s performance quietly erodes. A predictive-policing algorithm flags a neighborhood as high-crime. Police patrol more heavily. Arrests go up. The model “confirms” its own prediction.
The validation machinery of Chapter 6 — random train/test splits, cross-validation — was built for a specific world: rows are exchangeable, the deployment distribution matches the training distribution, and the model’s predictions don’t affect the outcome. That world is real, and when it holds, random-split validation is trustworthy. This chapter is about the cases where it doesn’t. There are four failure modes arranged from mildest to deepest:
Backtesting isn’t another thing that fails; for time-series forecasting, it is established and basically trustworthy methodology, with caveats for sudden regime shifts. Subgroup validation and distribution-shift testing likewise work when you do them. The thesis isn’t “validation is broken.” The thesis is: random-split validation is one tool, and you need the right tool for each failure mode.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error
import warnings
warnings.filterwarnings('ignore')
sns.set_style('whitegrid')
plt.rcParams['figure.figsize'] = (7, 4.5)
plt.rcParams['font.size'] = 12
DATA_DIR = 'data'When your validation metric looks great but deployment fails, the first place to check is the split. If your data has temporal structure — anything that evolves over time, from stock prices to air quality to NBA performance — the textbook random train/test split is silently broken. Observations sit close together in time, and those neighbors leak information across the split.
The EPA publishes daily Air Quality Index (AQI) measurements for every U.S. county. AQI measures how polluted the air is — 0 is pristine, 50 is “good,” and anything above 150 is “unhealthy.” In California, wildfire smoke can push AQI into the thousands.
Let’s load four years of data and focus on a few counties with good coverage.
# Load all four years of EPA data, filtered to California
dfs = []
for year in range(2021, 2025):
df_year = pd.read_csv(f'{DATA_DIR}/epa-air-quality/daily_aqi_by_county_{year}.csv')
df_year = df_year[df_year['State Name'] == 'California']
dfs.append(df_year)
aqi = pd.concat(dfs, ignore_index=True)
aqi['Date'] = pd.to_datetime(aqi['Date'])
aqi = aqi.sort_values(['county Name', 'Date']).reset_index(drop=True)
print(f"California records: {len(aqi):,}")
print(f"Counties: {aqi['county Name'].nunique()}")
print(f"Date range: {aqi['Date'].min().date()} to {aqi['Date'].max().date()}")California records: 75,449
Counties: 53
Date range: 2021-01-01 to 2024-12-31We’ll focus on three counties with good daily coverage and distinct air quality patterns: Los Angeles (urban smog), Sacramento (valley pollution), and Mono (extreme dust events).
# Pick 3 counties with good daily coverage and interesting patterns
counties = ['Los Angeles', 'Sacramento', 'Mono']
aqi_sub = aqi[aqi['county Name'].isin(counties)].copy()
# Check coverage
for county in counties:
n = aqi_sub[aqi_sub['county Name'] == county].shape[0]
print(f"{county}: {n} daily observations")Los Angeles: 1461 daily observations
Sacramento: 1461 daily observations
Mono: 1460 daily observationsBefore modeling anything, let’s look at the data. Notice the seasonal pattern — summer and fall have worse air quality (wildfire season). And notice the spikes: these are smoke events that dwarf normal variation.
fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
for ax, county in zip(axes, counties):
county_data = aqi_sub[aqi_sub['county Name'] == county]
ax.plot(county_data['Date'], county_data['AQI'], linewidth=0.5, alpha=0.8)
ax.set_ylabel('AQI')
ax.set_title(f'{county} County')
ax.axhline(y=150, color='red', linestyle='--', alpha=0.5, label='Unhealthy threshold')
ax.legend(loc='upper right')
axes[-1].set_xlabel('Date')
fig.suptitle('Daily AQI in California Counties (2021–2024)', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()
Mono County is wild — AQI values above 1,000 (and once above 8,000!) from extreme dust storms near Mono Lake. Los Angeles and Sacramento show the typical California pattern: seasonal wildfire spikes layered on top of baseline urban pollution.
A time series is non-stationary when its statistical properties (mean, variance) change over time. The AQI data above is non-stationary: the mean and variance shift with the seasons, and extreme events (wildfire spikes) appear unpredictably. The rest of this chapter explores why non-stationarity undermines standard evaluation methods.
Let’s try to predict tomorrow’s AQI. A natural approach: use lag features — yesterday’s AQI, the 7-day average, etc. This approach is just the feature engineering from Chapter 5, applied to time.
We’ll focus on Los Angeles for a clean demonstration.
# Work with Los Angeles
la = aqi_sub[aqi_sub['county Name'] == 'Los Angeles'][['Date', 'AQI']].copy()
la = la.sort_values('Date').reset_index(drop=True)
# Create lag features
la['aqi_lag1'] = la['AQI'].shift(1) # yesterday
la['aqi_lag7_avg'] = la['AQI'].rolling(7).mean().shift(1) # 7-day average (as of yesterday)
la['aqi_lag3_avg'] = la['AQI'].rolling(3).mean().shift(1) # 3-day average
la['day_of_year'] = la['Date'].dt.dayofyear # seasonality feature
# Drop rows with NaN from lagging
la = la.dropna().reset_index(drop=True)
print(f"LA dataset: {len(la)} days with features")
la.head(10)LA dataset: 1454 days with features| Date | AQI | aqi_lag1 | aqi_lag7_avg | aqi_lag3_avg | day_of_year | |
|---|---|---|---|---|---|---|
| 0 | 2021-01-08 | 88 | 84.0 | 96.857143 | 86.000000 | 8 |
| 1 | 2021-01-09 | 96 | 88.0 | 83.142857 | 87.333333 | 9 |
| 2 | 2021-01-10 | 71 | 96.0 | 85.857143 | 89.333333 | 10 |
| 3 | 2021-01-11 | 69 | 71.0 | 84.285714 | 85.000000 | 11 |
| 4 | 2021-01-12 | 69 | 69.0 | 83.142857 | 78.666667 | 12 |
| 5 | 2021-01-13 | 77 | 69.0 | 81.000000 | 69.666667 | 13 |
| 6 | 2021-01-14 | 86 | 77.0 | 79.142857 | 71.666667 | 14 |
| 7 | 2021-01-15 | 82 | 86.0 | 79.428571 | 77.333333 | 15 |
| 8 | 2021-01-16 | 93 | 82.0 | 78.571429 | 81.666667 | 16 |
| 9 | 2021-01-17 | 84 | 93.0 | 78.142857 | 87.000000 | 17 |
Recall from Chapter 6 that we split data randomly into training and test sets. That approach assumed each observation was independent of the others. Time-ordered observations violate that assumption.
Here’s what most tutorials (and most AI tools) would do — randomly shuffle the data and split 80/20. Let’s see how well it works.
# Random 80/20 split — the wrong way for time series
np.random.seed(42)
features = ['aqi_lag1', 'aqi_lag7_avg', 'aqi_lag3_avg', 'day_of_year']
mask = np.random.rand(len(la)) < 0.8
train_std = la[mask]
test_std = la[~mask]
model_std = LinearRegression()
model_std.fit(train_std[features], train_std['AQI'])
pred_std = model_std.predict(test_std[features])
r2_std = r2_score(test_std['AQI'], pred_std)
mae_std = mean_absolute_error(test_std['AQI'], pred_std)
print(f"Random split results:")
print(f" R-squared = {r2_std:.3f}")
print(f" MAE = {mae_std:.1f}")
print(f" Train size: {len(train_std)}, Test size: {len(test_std)}")
print()
print("Looks great! Ship it?")Random split results:
R-squared = 0.664
MAE = 15.4
Train size: 1151, Test size: 303
Looks great! Ship it?R² about 0.66 — not bad for a first attempt. But a hidden problem lurks in the evaluation itself.
Your team built an AQI forecasting model with R²=0.95 using a random train/test split. Before deploying it for public health alerts, would you trust that number? With a random split, the test set includes days from the middle of the time series — a test-set day from June 2022 has its neighbors (June 1, June 3, …) in the training set, and the lag features are computed from data right next to training points. What split strategy would you insist on, and why?
Data leakage occurs when information from the future (or from the test set) contaminates training, producing optimistically biased results. In time series, a random split causes leakage: test observations are highly correlated with nearby training observations. The lag feature “yesterday’s AQI” for a test-set day was likely in the training set.
The fix is simple: train on the past, test on the future. This protocol mimics how the model would actually be used — you’d train it on historical data and deploy it to predict tomorrow. The diagram below shows the difference: random splitting scatters test points throughout the timeline, while a temporal split draws a clean boundary.
# Visualize random vs temporal split
fig, axes = plt.subplots(2, 1, figsize=(10, 3), sharex=True)
dates = la['Date'].values
n = len(dates)
# Random split
colors_random = np.where(mask[:n], 'steelblue', 'orange')
axes[0].scatter(dates, np.ones(n), c=colors_random, s=2, alpha=0.6)
axes[0].set_yticks([])
axes[0].set_ylabel('Random\nsplit', rotation=0, ha='right', va='center')
axes[0].set_title('Train (blue) vs Test (orange): Random Split vs Temporal Split')
# Temporal split
cutoff = pd.Timestamp('2023-01-01')
colors_temporal = np.where(la['Date'] < cutoff, 'steelblue', 'orange')
axes[1].scatter(dates, np.ones(n), c=colors_temporal, s=2, alpha=0.6)
axes[1].axvline(x=cutoff, color='black', linewidth=1.5, linestyle='--', label='Train/test boundary')
axes[1].set_yticks([])
axes[1].set_ylabel('Temporal\nsplit', rotation=0, ha='right', va='center')
axes[1].legend(loc='upper right')
plt.tight_layout()
plt.show()
print("Top: random split scatters test points everywhere (leakage!).")
print("Bottom: temporal split puts all test points in the future (correct).")
Top: random split scatters test points everywhere (leakage!).
Bottom: temporal split puts all test points in the future (correct).We’re about to train on 2021–2022 and test on 2023. Will R-squared be higher, lower, or about the same as the random split? Why?
Let’s train on 2021–2022 and test on 2023.
# Temporal split — the right approach
train_temporal = la[la['Date'] < '2023-01-01']
test_temporal = la[(la['Date'] >= '2023-01-01') & (la['Date'] < '2024-01-01')]
model_temporal = LinearRegression()
model_temporal.fit(train_temporal[features], train_temporal['AQI'])
pred_temporal = model_temporal.predict(test_temporal[features])
r2_temporal = r2_score(test_temporal['AQI'], pred_temporal)
mae_temporal = mean_absolute_error(test_temporal['AQI'], pred_temporal)
print(f"Temporal split results (train: 2021-2022, test: 2023):")
print(f" R-squared = {r2_temporal:.3f}")
print(f" MAE = {mae_temporal:.1f}")
print(f" Train size: {len(train_temporal)}, Test size: {len(test_temporal)}")Temporal split results (train: 2021-2022, test: 2023):
R-squared = 0.686
MAE = 16.4
Train size: 723, Test size: 365Now let’s compare the two approaches side by side.
# Side-by-side comparison
print("=" * 50)
print(f"{'Method':<25} {'R-sq':>8} {'MAE':>8}")
print("=" * 50)
print(f"{'Random split (WRONG)':<25} {r2_std:>8.3f} {mae_std:>8.1f}")
print(f"{'Temporal split (RIGHT)':<25} {r2_temporal:>8.3f} {mae_temporal:>8.1f}")
print("=" * 50)
print(f"{'Leakage gap':<25} {r2_std - r2_temporal:>8.3f} {mae_temporal - mae_std:>8.1f}")
print()
gap = r2_std - r2_temporal
if gap > 0.01:
print(f"The random split inflated R-squared by {gap:.3f} -- that's the leakage.")
elif gap < -0.01:
print(f"Here the temporal split slightly outperforms (gap = {gap:.3f}).")
print("The leakage didn't help much -- but the correct evaluation is still the temporal one.")
else:
print(f"The gap is small ({gap:.3f}), but the temporal split is the only honest evaluation.")
print("Even when the numbers are similar, the random split's validity is compromised.")==================================================
Method R-sq MAE
==================================================
Random split (WRONG) 0.664 15.4
Temporal split (RIGHT) 0.686 16.4
==================================================
Leakage gap -0.021 1.0
Here the temporal split slightly outperforms (gap = -0.021).
The leakage didn't help much -- but the correct evaluation is still the temporal one.In this dataset, the two R² values are close — the random split didn’t dramatically inflate performance. But the random split is still the wrong evaluation. Its validity is compromised: test-set days have their neighbors in the training set, so the lag features get to “cheat.” The fact that the numbers happen to be similar here doesn’t justify the method. The temporal split is the only honest evaluation because it mimics how the model would actually be deployed.
The naive forecast predicts that tomorrow’s value equals today’s value: \(\hat{y}_{t+1} = y_t\). This baseline is the simplest possible forecast and the bar any model must clear. In the forecasting literature, other standard benchmarks include the seasonal naive method (predict the value from the same season last year) and the drift method (extrapolate the line between the first and last observations). These methods are collectively called benchmark forecasting methods.
# Naive forecast baseline: predict "same as yesterday"
naive_pred = test_temporal['aqi_lag1'] # yesterday's AQI IS the naive forecast
mae_naive_baseline = mean_absolute_error(test_temporal['AQI'], naive_pred)
r2_naive_baseline = r2_score(test_temporal['AQI'], naive_pred)
print(f"Naive baseline (same as yesterday):")
print(f" R-squared = {r2_naive_baseline:.3f}")
print(f" MAE = {mae_naive_baseline:.1f}")
print()
print(f"Linear model (temporal split):")
print(f" R-squared = {r2_temporal:.3f}")
print(f" MAE = {mae_temporal:.1f}")
print()
print("The naive baseline is surprisingly hard to beat!")
print("This is common in time series: yesterday is a strong predictor of today.")Naive baseline (same as yesterday):
R-squared = 0.662
MAE = 16.1
Linear model (temporal split):
R-squared = 0.686
MAE = 16.4
The naive baseline is surprisingly hard to beat!
This is common in time series: yesterday is a strong predictor of today.This result illustrates why backtesting is essential. In production, your model only ever sees the past. Your evaluation should reflect that constraint.
MAE tells you the average error in the same units as your target — useful for one series, but hard to compare across series with different scales. Is an MAE of 5 good for Los Angeles (typical AQI ~50) and Mono County (typical AQI ~30, but spikes to 8,000)?
MAPE (Mean Absolute Percentage Error) expresses errors as percentages: \[\text{MAPE} = \frac{1}{n} \sum_{i=1}^{n} \left| \frac{y_i - \hat{y}_i}{y_i} \right| \times 100\%.\] MAPE is intuitive — “the forecast is off by 15% on average” — but it breaks when actual values are near zero (division by zero) and penalizes over-forecasts more heavily than under-forecasts. The asymmetry comes from the bound: when you under-forecast, the percentage error is capped at 100% (you can’t predict less than zero); when you over-forecast, it’s unbounded. As a result, MAPE-optimal forecasts tend to be biased low. When the actual value is small, even a modest error produces a large percentage — MAPE disproportionately penalizes forecasts for low-value days.
MASE (Mean Absolute Scaled Error) compares your model’s error to the naive baseline’s error: \[\text{MASE} = \frac{\text{MAE of your model}}{\text{MAE of the naive forecast}}.\] MASE < 1 means you beat the naive forecast; MASE > 1 means you did worse than just predicting “same as yesterday.” This metric connects directly to the baseline comparison above.1
# Compute MAPE and MASE on the temporal test set
actual = test_temporal['AQI'].values
predicted = pred_temporal
naive = test_temporal['aqi_lag1'].values
# MAPE — watch out for zeros
nonzero = actual != 0
mape = np.mean(np.abs((actual[nonzero] - predicted[nonzero]) / actual[nonzero])) * 100
# MASE — scale by naive forecast error
mae_model = np.mean(np.abs(actual - predicted))
mae_naive_ref = np.mean(np.abs(actual - naive))
mase = mae_model / mae_naive_ref
print(f"LA County AQI forecasts (temporal split, 2023):")
print(f" MAE = {mae_model:.2f}")
print(f" MAPE = {mape:.1f}%")
print(f" MASE = {mase:.3f}")
print()
if mase < 1:
print(f"MASE = {mase:.3f} < 1: the model beats the naive baseline.")
else:
print(f"MASE = {mase:.3f} >= 1: the model does NOT beat the naive baseline.")
print("MAPE gives an intuitive percentage, but would blow up if AQI hit 0.")LA County AQI forecasts (temporal split, 2023):
MAE = 16.40
MAPE = 20.4%
MASE = 1.019
MASE = 1.019 >= 1: the model does NOT beat the naive baseline.
MAPE gives an intuitive percentage, but would blow up if AQI hit 0.A MASE near 1 means our linear model barely matches the naive “same as yesterday” baseline. That’s not a failure of the method — it’s a genuine finding. AQI is highly autocorrelated, so yesterday’s value is a strong predictor that’s hard to improve on with a simple linear model. More sophisticated models (nonlinear, with weather covariates) might do better, but the naive baseline sets a high bar.
A single temporal split gives you one number. But model performance can change over time. Walk-forward validation trains on an expanding window and tests on the next period, giving you a performance curve. In the forecasting literature, this is called time series cross-validation — the same idea as k-fold cross-validation, but respecting temporal order.
Here’s the idea. Each row is one “fold.” The training set (blue) expands forward in time, and the model is tested on the next month (orange):
Fold 1: [===== TRAIN =====][TEST]
Fold 2: [======= TRAIN =======][TEST]
Fold 3: [========= TRAIN =========][TEST]
Fold 4: [=========== TRAIN ===========][TEST]
2021 2022 2023We used an expanding window here (all past data). An alternative is a sliding window (fixed-width training set that moves forward), which can be better when very old data is from a different regime — exactly the distribution shift problem we’ll see later.
This distinction is especially important when the data has regime changes — like wildfire season.
# Walk-forward: train on expanding window, test on next month
# (freq='MS' = month start: generates the 1st of each month)
results = []
test_months = pd.date_range('2023-01-01', '2024-01-01', freq='MS')
for i, test_start in enumerate(test_months[:-1]):
test_end = test_months[i + 1]
train_wf = la[la['Date'] < test_start]
test_wf = la[(la['Date'] >= test_start) & (la['Date'] < test_end)]
if len(test_wf) < 10 or len(train_wf) < 100:
continue
model_wf = LinearRegression()
model_wf.fit(train_wf[features], train_wf['AQI'])
pred_wf = model_wf.predict(test_wf[features])
month_mae = mean_absolute_error(test_wf['AQI'], pred_wf)
month_r2 = r2_score(test_wf['AQI'], pred_wf)
results.append({
'month': test_start,
'mae': month_mae,
'r2': month_r2,
'n_test': len(test_wf),
'max_aqi': test_wf['AQI'].max()
})
print(f" {test_start.strftime('%b %Y')}: MAE = {month_mae:.1f}, R² = {month_r2:.3f}")
wf_results = pd.DataFrame(results)
wf_results['month_label'] = wf_results['month'].dt.strftime('%b %Y')
print(wf_results[['month_label', 'mae', 'r2', 'max_aqi']].to_string(index=False)) Jan 2023: MAE = 9.3, R² = 0.130
Feb 2023: MAE = 9.6, R² = -0.260
Mar 2023: MAE = 6.8, R² = -0.607
Apr 2023: MAE = 14.9, R² = 0.479
May 2023: MAE = 15.0, R² = 0.576
Jun 2023: MAE = 17.4, R² = 0.309
Jul 2023: MAE = 27.6, R² = -0.145
Aug 2023: MAE = 23.8, R² = 0.609
Sep 2023: MAE = 22.3, R² = 0.187
Oct 2023: MAE = 20.6, R² = 0.507
Nov 2023: MAE = 14.3, R² = 0.222
Dec 2023: MAE = 13.6, R² = 0.199
month_label mae r2 max_aqi
Jan 2023 9.266960 0.130164 82
Feb 2023 9.586316 -0.260122 97
Mar 2023 6.789728 -0.607446 64
Apr 2023 14.934882 0.479284 156
May 2023 14.992056 0.575711 136
Jun 2023 17.394106 0.309244 192
Jul 2023 27.599075 -0.144907 210
Aug 2023 23.828699 0.608855 208
Sep 2023 22.339099 0.186713 164
Oct 2023 20.577634 0.507455 172
Nov 2023 14.298563 0.222461 164
Dec 2023 13.556200 0.199320 153The monthly results tell a richer story than a single annual number. Let’s visualize the performance over time.
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
ax1.bar(wf_results['month_label'], wf_results['mae'], color='steelblue', alpha=0.8)
ax1.set_ylabel('MAE')
ax1.set_title('Walk-Forward Validation: Monthly MAE (Los Angeles, 2023)')
ax1.axhline(y=mae_temporal, color='red', linestyle='--', alpha=0.5,
label=f'Full-year MAE = {mae_temporal:.1f}')
ax1.legend()
ax2.bar(wf_results['month_label'], wf_results['max_aqi'], color='coral', alpha=0.8)
ax2.set_ylabel('Max AQI in Month')
ax2.set_title('Worst Air Quality Day Each Month')
ax2.axhline(y=150, color='red', linestyle='--', alpha=0.5, label='Unhealthy threshold')
ax2.legend()
plt.xticks(rotation=45, ha='right')
plt.tight_layout()
plt.show()
Notice how the model’s error spikes during summer months when wildfire smoke events push AQI to extremes. The model was trained on mostly “normal” air quality patterns and hasn’t seen enough extreme events to handle them.
This pattern is the value of walk-forward validation: it reveals when your model breaks, not just how well it does on average.
Recall from Chapter 5 that feature engineering transforms raw data into more useful predictors. Lag features are the same idea applied to time:
| Feature | What it captures |
|---|---|
| Yesterday’s AQI | Short-term momentum |
| 3-day average | Recent trend |
| 7-day average | Weekly baseline |
| Day of year | Seasonal pattern |
| Same day last year | Year-over-year pattern |
Let’s see which features matter most.
# Which features contribute most?
print("Feature coefficients (temporal model):")
for feat, coef in zip(features, model_temporal.coef_):
print(f" {feat:<20s} {coef:+.4f}")
print(f" {'intercept':<20s} {model_temporal.intercept_:+.4f}")
print()
print("Yesterday's AQI dominates -- AQI is highly autocorrelated.")
print("The seasonal feature (day_of_year) adds a small boost.")Feature coefficients (temporal model):
aqi_lag1 +0.7579
aqi_lag7_avg +0.2048
aqi_lag3_avg -0.1952
day_of_year +0.0103
intercept +18.2827
Yesterday's AQI dominates -- AQI is highly autocorrelated.
The seasonal feature (day_of_year) adds a small boost.You might be thinking: “This example is an air quality problem. My field is tech/finance/sports.” The temporal-leakage trap appears anywhere data has a time dimension. Try predicting how many points an NBA player scores tonight from their last 5 games.
# Load NBA data and keep players with full three-season coverage
nba = pd.read_csv(f'{DATA_DIR}/nba/nba_load_management.csv')
nba['GAME_DATE'] = pd.to_datetime(nba['GAME_DATE'])
nba = nba.sort_values(['PLAYER_NAME', 'GAME_DATE']).reset_index(drop=True)
active = nba.groupby('PLAYER_NAME').size()
nba_active = nba[nba['PLAYER_NAME'].isin(active[active >= 150].index)].copy()
# Rolling lag features per player
nba_active['pts_lag1'] = nba_active.groupby('PLAYER_NAME')['PTS'].shift(1)
nba_active['pts_lag5_avg'] = nba_active.groupby('PLAYER_NAME')['PTS'].transform(
lambda x: x.rolling(5).mean().shift(1))
nba_active['min_lag5_avg'] = nba_active.groupby('PLAYER_NAME')['MIN'].transform(
lambda x: x.rolling(5).mean().shift(1))
nba_model = nba_active.dropna(subset=['pts_lag1', 'pts_lag5_avg', 'min_lag5_avg']).copy()
nba_features = ['pts_lag1', 'pts_lag5_avg', 'min_lag5_avg', 'HOME']
print(f"Modeling dataset: {len(nba_model):,} player-games; seasons {list(nba_model['SEASON'].unique())}")Modeling dataset: 47,991 player-games; seasons ['2021-22', '2022-23', '2023-24']Repeat the random-vs-temporal experiment: random 80/20 split versus train on seasons 1–2, test on season 3.
np.random.seed(42)
mask_nba = np.random.rand(len(nba_model)) < 0.8
model_r = LinearRegression().fit(nba_model[mask_nba][nba_features], nba_model[mask_nba]['PTS'])
r2_r = r2_score(nba_model[~mask_nba]['PTS'], model_r.predict(nba_model[~mask_nba][nba_features]))
train_nba = nba_model[nba_model['SEASON'].isin(['2021-22', '2022-23'])]
test_nba = nba_model[nba_model['SEASON'] == '2023-24']
model_t = LinearRegression().fit(train_nba[nba_features], train_nba['PTS'])
r2_t = r2_score(test_nba['PTS'], model_t.predict(test_nba[nba_features]))
print(f"Random split R²: {r2_r:.3f}")
print(f"Temporal split R²: {r2_t:.3f} (gap {r2_r - r2_t:+.3f})")Random split R²: 0.468
Temporal split R²: 0.484 (gap -0.017)Unlike AQI, the NBA model generalizes well across seasons — player scoring patterns are more stable year to year than air quality is. When temporal structure doesn’t create strong distribution shift, the random split may not inflate performance. But the temporal split is still the correct evaluation: it mirrors how you’d actually deploy.
One more twist before we leave temporal leakage behind. Even a model with a decent R² doesn’t translate to profit in a betting market — the market already incorporates the information your model uses. If your “edge” comes from the last 5 games’ scoring average, so does everyone else’s line, and the market has already priced it in. This is the first hint of a deeper problem: the market isn’t a static environment. It adapts to public signals. Once the whole world knows a feature predicts scoring, the feature stops predicting market outcomes. The temporal split says your R² is real; the betting profit says the market isn’t sitting still. That gap between statistical performance and deployment performance is the bridge to Part 2.
Let’s go back to air quality and look at something the model has never seen.
Mono County, near Mono Lake in eastern California, occasionally experiences extreme dust/wildfire events that push AQI into the thousands — sometimes above 8,000. That’s roughly 16x the “Hazardous” threshold of 500.
If we train a model on “normal” days (AQI below 300), what would it predict for an AQI of 8,000? How far off would it be?
# Mono County extremes
mono = aqi_sub[aqi_sub['county Name'] == 'Mono'][['Date', 'AQI']].copy()
mono = mono.sort_values('Date').reset_index(drop=True)
# Create lag features
mono['aqi_lag1'] = mono['AQI'].shift(1)
mono['aqi_lag7_avg'] = mono['AQI'].rolling(7).mean().shift(1)
mono['aqi_lag3_avg'] = mono['AQI'].rolling(3).mean().shift(1)
mono['day_of_year'] = mono['Date'].dt.dayofyear
mono = mono.dropna().reset_index(drop=True)
print("Mono County AQI distribution:")
print(mono['AQI'].describe().to_string())
print(f"\nDays above 500: {(mono['AQI'] > 500).sum()}")
print(f"Max AQI: {mono['AQI'].max()}")Mono County AQI distribution:
count 1453.000000
mean 57.946318
std 351.039623
min 0.000000
25% 14.000000
50% 23.000000
75% 39.000000
max 8368.000000
Days above 500: 16
Max AQI: 8368Let’s train a linear model on “normal” days (AQI below 300) and see what it predicts during extreme events.
# Train a model on "normal" data (AQI < 300) and see what it predicts for extremes
normal_mono = mono[mono['AQI'] < 300]
extreme_mono = mono[mono['AQI'] >= 300]
model_mono = LinearRegression()
model_mono.fit(normal_mono[features], normal_mono['AQI'])
if len(extreme_mono) > 0:
pred_extreme = model_mono.predict(extreme_mono[features])
print("Model predictions vs reality during extreme events:")
print(f"{'Actual AQI':>12} {'Predicted':>12} {'Error':>12}")
print("-" * 40)
for actual, pred in list(zip(extreme_mono['AQI'].values, pred_extreme))[:10]:
print(f"{actual:>12.0f} {pred:>12.0f} {actual - pred:>12.0f}")
coef_lag1 = model_mono.coef_[features.index('aqi_lag1')]
print()
print(f"Coefficient learned on aqi_lag1 (trained on AQI < 300): {coef_lag1:.4f}")
else:
print("No extreme events found in Mono County subset.")Model predictions vs reality during extreme events:
Actual AQI Predicted Error
----------------------------------------
1196 29 1167
662 29 633
2971 30 2941
414 29 385
761 29 732
304 30 274
3404 30 3374
768 30 738
7835 30 7805
2050 38 2012
Coefficient learned on aqi_lag1 (trained on AQI < 300): -0.0004The numbers are stark, but the visualization makes the failure visceral.
# Visualize the model's failure during extreme events
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(mono['Date'], mono['AQI'], linewidth=0.5, alpha=0.7, label='Actual AQI')
# Predict on all data
pred_all = model_mono.predict(mono[features])
ax.plot(mono['Date'], pred_all, linewidth=0.5, alpha=0.7, color='orange',
label='Model prediction')
ax.axhline(y=500, color='red', linestyle='--', alpha=0.5, label='Hazardous (500)')
ax.set_ylabel('AQI')
ax.set_xlabel('Date')
ax.set_title('Mono County: Model vs Reality During Extreme Events')
ax.legend()
plt.tight_layout()
plt.show()
print("The orange line (predictions) can't reach the spikes.")
print("A linear model trained on normal conditions extrapolates nonsensically.")
The orange line (predictions) can't reach the spikes.
A linear model trained on normal conditions extrapolates nonsensically.Why does the model fail even on days when its input — yesterday’s AQI — is itself extreme? Look at the learned coefficient on aqi_lag1: essentially zero. Inside the training regime (AQI below 300), day-to-day changes in Mono County air quality are dominated by noise, so ordinary least squares learns that yesterday’s value carries almost no predictive signal and settles for “predict the training mean.” When we plug in aqi_lag1 = 8000, the prediction barely budges from the training mean — the model never had a chance. The failure isn’t that the model refuses to extrapolate; it’s that in the training distribution, the relevant signal didn’t exist to be learned.
Distribution shift occurs when future data comes from a fundamentally different regime than the training data. It is a severe form of non-stationarity — not just gradual changes in the mean, but a qualitative change the training data never captured.
Distribution shift is the deepest problem in time series forecasting:
“Prediction is very difficult, especially about the future.” — attr. Niels Bohr
A clinical trial on adults can’t tell you about children — the inference only covers the population you sampled. Same idea here: a model trained on AQI 0–300 can’t tell you what happens at AQI 8,000. You can’t extrapolate beyond the range of your data. Chapter 6 previewed this: random splits hold out observations that look like the ones you trained on, so cross-validation protects you against in-distribution noise, not against a future that doesn’t look like the past.
In Chapter 17, we’ll see that AutoML tools often make this same mistake — they use standard cross-validation by default and don’t know your data is temporal.
Distribution shift can also be gradual — styles of play, market conditions, or climate patterns that evolve slowly can erode model performance without triggering obvious alarms.
If models can fail silently when conditions change, point forecasts hide exactly what matters. The next section quantifies the uncertainty around a forecast — at least under the regime the data captures.
A point forecast — “tomorrow’s AQI will be 62” — hides everything that matters for decision-making. Should the county issue a health advisory? That depends on whether 62 means “almost certainly between 55 and 70” or “could easily be 200.” Prediction intervals give you that range.
The idea: use the bootstrap. Recall from Chapter 8 that bootstrap resampling estimates uncertainty by simulating many alternative datasets from the data we have. Here we apply the same idea to forecast uncertainty: we resample the model’s training residuals, add them to each forecast, and repeat many times. The spread of the resulting forecasts gives you a prediction interval.
# Bootstrap prediction intervals for the temporal test set
np.random.seed(42)
train_resid = train_temporal['AQI'].values - model_temporal.predict(train_temporal[features])
n_boot = 1000
# Generate bootstrap forecasts by resampling residuals
boot_preds = np.zeros((n_boot, len(test_temporal)))
point_pred = model_temporal.predict(test_temporal[features])
for b in range(n_boot):
resampled_resid = np.random.choice(train_resid, size=len(test_temporal), replace=True)
boot_preds[b, :] = point_pred + resampled_resid
# 95% prediction interval: 2.5th and 97.5th percentiles
lower = np.percentile(boot_preds, 2.5, axis=0)
upper = np.percentile(boot_preds, 97.5, axis=0)
# How often does the actual value fall inside the interval?
coverage = np.mean((test_temporal['AQI'].values >= lower) & (test_temporal['AQI'].values <= upper))
avg_width = np.mean(upper - lower)
print(f"95% bootstrap prediction interval:")
print(f" Coverage: {coverage:.1%} of test days (target: 95%)")
print(f" Average width: {avg_width:.1f} AQI points")95% bootstrap prediction interval:
Coverage: 93.7% of test days (target: 95%)
Average width: 89.5 AQI pointsLet’s visualize the prediction intervals for the first two months of the test period. The shaded band shows the 95% interval around each point forecast.
# Visualize prediction intervals for a 60-day window
fig, ax = plt.subplots(figsize=(10, 5))
window = slice(0, 60)
dates_w = test_temporal['Date'].values[window]
ax.fill_between(dates_w, lower[window], upper[window],
alpha=0.25, color='steelblue', label='95% prediction interval')
ax.plot(dates_w, point_pred[window], color='steelblue', linewidth=1.2, label='Point forecast')
ax.plot(dates_w, test_temporal['AQI'].values[window], color='black',
linewidth=0.8, alpha=0.8, label='Actual')
ax.set_ylabel('AQI')
ax.set_xlabel('Date')
ax.set_title('LA County AQI: Forecast with 95% Prediction Interval (Jan–Feb 2023)')
ax.legend()
plt.tight_layout()
plt.show()
print("Point forecasts alone are dangerous for decision-making.")
print("Always report intervals so stakeholders can plan for the plausible range.")
Point forecasts alone are dangerous for decision-making.
Always report intervals so stakeholders can plan for the plausible range.This resampling approach is completely general: you can construct prediction intervals for any model, not just linear regression. The key principle is that a forecast without a measure of uncertainty isn’t a complete forecast. (The intervals above assume stationary residuals — in reality, bands should widen during high-variance regimes like wildfire season.)
Temporal leakage is a validation technique failure — random splits saw the future. Feedback loops are a deeper failure: even with a perfectly clean temporal split, your model can fail in deployment because your model is part of the world it predicts. What happens when a prediction model doesn’t just reflect the world but changes it?
Predictive policing. A model predicts high crime in neighborhood X. Police patrol X more heavily. More patrols produce more arrests in X. The model retrains on the new data and “confirms” that X is high-crime. The prediction is self-fulfilling — the model creates the data that justifies itself. This is not confounding in the usual sense. In a confounding story, the association already exists and we misinterpret it. In a feedback loop, the model causes the association.
Credit scoring. In Weapons of Math Destruction, Cathy O’Neil describes how credit scores can trap people in cycles they can’t escape.2 A low credit score leads to denied loans. Without loans, a borrower can’t build credit history. The score stays low — not because the person is irresponsible, but because the scoring system denied them the opportunity to demonstrate otherwise. The score creates the reality it claims to measure.
Recommendation algorithms. YouTube recommends increasingly extreme content because users engage with it. Users watch what’s recommended. The algorithm learns they “prefer” extreme content and recommends more. The model shapes the preferences it claims to predict.
Clinical alerts. A hospital deploys a sepsis-prediction model trained on historical data. Once the alerts go live, clinicians treat flagged patients earlier. Fewer patients progress to sepsis. The model, evaluated after deployment, looks worse than it did on the validation set — not because it degraded, but because the thing it was predicting stopped happening. If the hospital responds by retraining on post-deployment data, the model’s ground truth is now shaped by the model itself.
Betting markets. Back to the NBA example. Your model beats a naive baseline on historical data. You deploy it against the sportsbook. The line already encodes the same information your model uses — recent scoring averages, rest days, matchup difficulty — because thousands of sharp bettors have moved the line. Your “edge” disappears the moment the market has processed the same signals. The model doesn’t fail because of distribution shift; it fails because the market is the feedback loop, adapting to public information in real time.
In each case, there is a causal arrow running from the model’s prediction back to the outcome:
fig, ax = plt.subplots(figsize=(10, 3.2))
ax.set_xlim(0, 10); ax.set_ylim(0, 3); ax.axis('off')
nodes = [('Training\ndata', 1.0, 2.3), ('Model', 2.8, 2.3), ('Prediction', 4.6, 2.3),
('Action', 6.4, 2.3), ('Outcome', 8.2, 2.3), ('New training\ndata', 9.5, 2.3)]
for name, x, y in nodes:
ax.add_patch(plt.Rectangle((x - 0.65, y - 0.3), 1.3, 0.6,
facecolor='lightblue', edgecolor='navy', lw=1.8, zorder=3))
ax.text(x, y, name, ha='center', va='center', fontsize=10, fontweight='bold', zorder=4)
for (_, x1, y1), (_, x2, y2) in zip(nodes[:-1], nodes[1:]):
ax.annotate('', xy=(x2 - 0.65, y2), xytext=(x1 + 0.65, y1),
arrowprops=dict(arrowstyle='->', lw=2, color='navy'))
# Cycle return arrow: from "New training data" back to "Training data"
ax.annotate('', xy=(1.0, 2.0), xytext=(9.5, 1.6),
arrowprops=dict(arrowstyle='->', lw=2, color='#d62728',
connectionstyle='arc3,rad=-0.3'))
ax.text(5.0, 0.6, 'cycle repeats', ha='center', va='center',
fontsize=10, color='#d62728', fontstyle='italic')
plt.tight_layout(); plt.show()
This is why “correlation is not causation” isn’t the full story. Sometimes the model is the cause. When a model’s predictions change the environment that generates future data, standard validation breaks down. The test set no longer represents the world the model will operate in, because the model itself reshapes that world.
Whenever you deploy a model, ask: how will the predictions change the data distribution in the future? If the answer is “they won’t” (a weather forecaster doesn’t affect the weather, a particle-physics model doesn’t affect the particles), you’re safe from feedback loops. If the answer is “they will” (credit scores affect who gets loans, sepsis alerts affect who gets treated), you have a feedback loop — and your validation set can’t tell you how it will behave.
The phenomenon we’re calling “feedback loop” has a formal research framework called performative prediction (Perdomo, Zrnic, Mendler-Dünner, Hardt 2020). It treats the deployed model \(f_\theta\) as inducing a data distribution \(\mathcal{D}(\theta)\) that depends on the deployed parameters, and studies when iterative retraining converges to a fixed point (a “performatively stable” \(\theta\) where the model is accurate against the distribution it induces). The framework also subsumes strategic classification, where individual agents deliberately change their features to score better. The formal machinery is mostly relevant for algorithm designers; the phenomenon is what matters for analysts, which is why we stick with “feedback loop” in this chapter.
Cathy O’Neil calls the worst cases Weapons of Math Destruction (WMDs).3 Three properties, together, turn an ordinary predictive model into a WMD:
College rankings are a WMD: high rankings attract students, faculty, and donors, which raises the school’s quality — the ranking creates the reality it claims to measure. Parole-risk models are a WMD: keeping someone in prison longer may cause them to reoffend (fewer job prospects, deeper criminal networks), confirming the model’s prediction. By contrast, CERN’s particle-physics models are not a WMD: the outcome is measurable, the experiments are rigorous, and the predictions don’t change the particles.
You build a model that predicts which students will fail a course, and the university uses it to assign tutoring resources. Is this a Weapon of Math Destruction? Check the three properties: is the outcome measurable? Could the prediction harm students (e.g., through stigma or self-fulfilling expectations)? Does it create a feedback loop? What design choices would make this a virtuous intervention rather than a vicious cycle?
“When a measure becomes a target, it ceases to be a good measure” (attributed to Charles Goodhart, 1975; popularized by Marilyn Strathern). The law describes something deeper than overfitting. Its core is game-theoretic: once a measurement is tied to consequences, the units being measured respond to the measurement. Those who can game the metric most effectively benefit disproportionately. The data the metric generates after it becomes a target is no longer the data it was calibrated on — it’s the data that results from agents trying to score well on it.
Three ingredients define a Goodhart setup:
Before the consequences are attached, \(M\) and \(G\) may be tightly linked. After, \(M\) reflects the mix of “underlying \(G\)” and “skill at gaming \(M\)” — and the more optimization pressure you apply, the more the second term dominates.
The pattern appears everywhere once you see it.
Notice the shared structure: the measurement changes the distribution of the thing being measured, and how much it changes depends on which units are best at responding to the measurement. You can’t read Goodhart as “the data was fine, we just optimized too hard.” The data itself is different once the measurement goes live.
The famous ML phenomenon — train loss goes down while test loss rises — is a shadow of Goodhart, not its core. Go back to Chapter 6. When you train a high-capacity model on a fixed training set, you’re the sole agent optimizing a proxy (training loss) for a goal (generalization loss). Your “gaming” happens entirely inside the learning algorithm: the model memorizes idiosyncratic noise in the training rows, which looks great on the training metric but doesn’t transfer. Cross-validation — reserving data the optimizer never sees — is how ML practice resists this version.
But this is the static case: the data is fixed, only the model is optimizing. The full Goodhart picture is dynamic: the data itself is generated by strategic agents who respond to the metric after it is deployed. A recidivism-risk score that detains suspected flight risks changes which crimes go unobserved. A sepsis-alert model that triggers early intervention changes which patients develop sepsis. A ranking signal for ad quality changes which ads advertisers write. In every case, what gets to the training set next round depends on what the model did last round — and which agents are best at adapting. No amount of held-out test data from before deployment tells you what the deployed system does.
If you repeatedly tune hyperparameters against the validation set — if you let the validation set become the target — even the ML-internal version of Goodhart kicks in: the validation loss stops measuring generalization. This is why serious ML practice distinguishes training / validation / test sets, and why the test set is used once, after all decisions are locked.
The pattern is the same whether the optimizer is a gradient-descent run, a hospital administrator, or a student preparing for a standardized test: the measurement is tied to consequences, the measured party responds, and some respond more effectively than others. Whenever you use a number to make decisions, ask two questions:
If the answer to the first is yes, you have a Goodhart problem. If the answer to the second is “not the intended beneficiaries,” you have an equity problem on top of a Goodhart problem. Defensive designs include auditing with a second metric the optimizer can’t see, capping how often any one party can update against the metric, or — the strongest defense, and the subject of Chapter 19 — randomizing a holdout so you can still see what the unperturbed world would look like.
Temporal leakage, feedback loops, and Goodhart’s law are distinct failure modes, but they share a recipe for defense.
Before deployment:
After deployment:
Chapter 17 closes this course with a 15-item checklist for evaluating any analysis. Items #13 (Goodhart) and #14 (feedback loops) refer back here. When you see a production model that’s been running for a year, ask: which of the three failure modes would catch it, if any?
pd.Series.shift(n) — create lag features by shifting values forward or backwardpd.Series.rolling(k).mean() — compute rolling window averagesdf[df['Date'] < cutoff])You should be able to: (1) explain why random splitting fails for time series data and implement a temporal split, (2) describe walk-forward validation (expanding and sliding window variants), (3) explain what distribution shift is and why no cross-validation method can fully address it, (4) identify the three properties of a WMD given a deployment scenario, (5) give two examples of feedback loops and diagram the prediction-action-outcome cycle, (6) state Goodhart’s law and explain how overfitting is a special case, (7) list three pre-deployment and three post-deployment practices that defend against these failure modes, and (8) explain why prediction intervals matter and how bootstrap resampling produces them.
Strictly, MASE uses the in-sample naive MAE as denominator (Hyndman & Koehler, 2006); we use the test-set version for simplicity.↩︎
O’Neil, Weapons of Math Destruction: How Big Data Increases Inequality and Threatens Democracy, Crown, 2016.↩︎
O’Neil, Weapons of Math Destruction, Crown, 2016. O’Neil’s own three properties are opacity, scale, and damage; we focus here on a related triad — inspired by her framework — that highlights feedback dynamics.↩︎