Time Series & Leakage

Assumes you know from ML-101

This lesson builds on ML-101 Lesson 9 (Model Evaluation), ML-101 Lesson 3 (Linear Regression), and ML-101 Lesson 10 (Overfitting, Bias & Variance). The reader is assumed to know that a model is evaluated on a held-out set, that $k$-fold cross-validation partitions a dataset into $k$ disjoint folds and reports average fold performance, and that ordinary least squares fits a linear model by minimising squared error. The reader is also presumed to know that overfitting occurs when a model memorises noise in the training set and generalises poorly.

What ML-101 did not do is warn that cross-validation on time-ordered data silently produces optimistic scores because the future leaks into the training fold. It also did not derive why MAPE is asymmetric and biased toward low-volume series, or what autocorrelation does to the standard error of a fitted regression coefficient. This lesson derives walk-forward validation from the causal ordering, derives the autocorrelation structure that breaks naive inference, and audits every forecasting metric by the assumption it makes about the scale of the data.

Learning Objectives

1. Derive walk-forward validation as the only unbiased evaluation strategy when the data has temporal dependence, and show why random $k$-fold produces optimistic scores by future-leakage.
2. Compute the autocorrelation function and the partial autocorrelation function for a stationary series, and explain why they require stationarity and what to do when it fails.
3. Compare MAPE, RMSE and MASE as forecasting metrics, deriving the bias and scale-sensitivity of each.
4. Implement a leakage audit on a forecasting pipeline by checking every feature's date lag, every transform's fit-on-train discipline, and every model's lookback window.
5. Fit a baseline forecaster (seasonal naïve, drift, exponential smoothing) and decide when a learned model is justified over it, by out-of-sample MAPE or MASE.

Why time series breaks cross-validation

The standard cross-validation shuffle assigns rows to folds uniformly at random. For independent observations this is unbiased: each fold is a representative sample of the marginal distribution, and the average fold score estimates the held-out score.

For a series $\{y_t\}_{t=1}^{T}$ observed in time order, the random shuffle breaks the temporal structure. Each fold contains rows from many different time points, and a row in the test fold may have a near neighbour in the training fold that is earlier or later in time. When the series is positively autocorrelated — most economic and physical series are — the near neighbour in the training fold is a good predictor of the test row, and the model's score on the test fold is optimistic.

The mechanism is concrete. Take a series with lag-1 autocorrelation $\rho_1 > 0$, and let the test fold contain the row at time $t$. The training fold contains, with probability $1/k$, the row at time $t-1$ (and $t-2$, and so on for short lags). The model sees $y_{t-1}$ as a feature — either directly or through any feature that aggregates recent observations — and uses it to predict $y_t$. The cross-validated score reflects the model's ability to look one step into the past, not its ability to forecast from genuinely older data.

For a series with $\rho_1 = 0.95$ — typical of monthly economic indicators — the random-fold score overstates the true forecasting performance by an order of magnitude. Practitioners who trust those scores build models that fail on the day they are deployed, because the production setting is always a forecasting setting: the model is asked to predict $y_{t+h}$ from information available at time $t$, not from a bag of nearby time points that include the future.

The fix is to enforce temporal ordering in the fold structure. Validation rows must come from times later than every training row. The simplest version is a single train/test split with the test set in the future of the training set; the general version is walk-forward validation.

import numpy as np

def naive_kfold_cv(y: np.ndarray, k: int) -> np.ndarray:
    """
    The standard shuffled k-fold: WRONG for time series.

    Each fold contains rows from many time points, so the validation rows
    have near-time neighbours in the training fold. For an autocorrelated
    series the score is optimistic.
    """
    n = len(y)
    idx = np.arange(n)
    rng = np.random.default_rng(0)
    rng.shuffle(idx)
    return np.array_split(idx, k)

The function looks identical to the standard sklearn KFold(shuffle=True).split(X) — there is no warning that it is silently wrong on temporal data. The lesson's whole point is that the leakage is invisible in code: the function runs, returns a number, and the number is too good.

Walk-forward validation

Walk-forward validation enforces the temporal order. The series is split into a sequence of expanding training windows and a fixed or expanding validation window. At step $i$, the model is trained on $\{y_t : t \le t_i\}$ and validated on $\{y_t : t_i < t \le t_{i+1}\}$.

The most general form has three parameters:

• an initial training size $n_0$ — enough data to fit the model,
• a step size $\Delta$ — how often the validation window advances,
• a horizon $h$ — how far ahead the model is being asked to forecast.

The split at step $i$ is

A common practical choice is $\Delta = h$: the validation window advances by exactly one horizon, so each step evaluates a forecast at a single horizon, and the windows do not overlap. The score is the average of the per-step scores; the variance of that average is the variance across steps, not the variance within a step.

For a horizon of one, the rolling-origin or expanding-window strategy is to fit on $\{y_1, \dots, y_{t-1}\}$ and forecast $y_t$, then advance $t$ by one. This is the most expensive but the least biased evaluation: every forecast is a genuine $h=1$ ahead prediction conditioned on data the model has not seen.

import numpy as np
from typing import Callable

def walk_forward(y: np.ndarray, fit_predict: Callable,
                 n_init: int, horizon: int = 1,
                 step: int | None = None) -> list[float]:
    """
    Walk-forward validation with expanding training window.

    `fit_predict` is a callable (history: np.ndarray, h: int) -> np.ndarray
    that fits on `history` and forecasts `h` steps ahead. Returns the
    per-step forecast errors.
    """
    step = step or horizon
    errs = []
    for t in range(n_init, len(y) - horizon + 1, step):
        history = y[:t]
        forecast = fit_predict(history, horizon)
        truth = y[t: t + horizon]
        errs.append(np.abs(forecast - truth).mean())
    return errs

The crucial property of this function is that no information from $y_t$ or later is ever passed to the fit step. The forecast for time $t$ is conditioned on $\{y_1, \dots, y_{t-1}\}$ only. The score is the average across $t$ of the forecast error; under the assumption that the test distribution matches the training distribution, this average estimates the deployment score honestly.

Walk-forward is not free. For a series of length $T$ with initial training $n_0$ and horizon $h$, the loop runs $\lfloor (T - n_0)/h \rfloor$ times. For a deep model that takes a minute to fit, this can be hours. Two practical mitigations: choose $\Delta > h$ so the windows skip, and choose a fast baseline for the burn-in. The lesson's deeper point is that the cost is the price of an honest evaluation; any cheaper evaluation has leakage baked in.

Autocorrelation and stationarity

The autocorrelation function (ACF) at lag $k$ is the correlation of the series with itself shifted by $k$:

The partial autocorrelation function (PACF) at lag $k$ is the correlation between $y_t$ and $y_{t-k}$ after regressing out the intermediate lags $y_{t-1}, \dots, y_{t-k+1}$. Equivalently,

The ACF and PACF are the diagnostic plots of time-series modelling. An AR($p$) process has a PACF that cuts off after lag $p$ and an ACF that tails off; an MA($q$) process has an ACF that cuts off after lag $q$ and a PACF that tails off. Both are derived from the Yule-Walker equations and are central to identifying the order of an ARIMA fit.

The catch is stationarity. The ACF is well-defined only when the mean $\mu$ and variance $\sigma^2$ are time-invariant. If the series has a trend, the ACF does not decay — even an i.i.d. random walk with a linear drift has an ACF that approaches 1 — and the diagnostic plots are unreadable. The fix is to difference the series until the ACF decays, typically with the augmented Dickey-Fuller test as the formal check.

The standard error of a regression coefficient assumes independent residuals. Under autocorrelated residuals, the OLS estimate is still unbiased but the standard error is wrong — usually too small, which is why naive inference produces overconfident conclusions on time-series data. The fix is Newey-West standard errors, which add a correction that depends on the estimated autocorrelation structure. The lesson's point: any regression on time-series data without an autocorrelation-robust standard error is reporting a confidence interval that does not mean what it claims.

import numpy as np

def acf(y: np.ndarray, max_lag: int) -> np.ndarray:
    """
    Sample autocorrelation function up to max_lag.

    Subtract the mean, normalise by the variance, and correlate with the
    lag-k version of itself. The first entry is 1 by construction.
    """
    yc = y - y.mean()
    var = (yc ** 2).sum()
    out = np.empty(max_lag + 1)
    for k in range(max_lag + 1):
        out[k] = (yc[k:] * yc[:len(y) - k]).sum() / var
    return out

The function is a literal transcription of the ACF definition. Its assumption is that the series is (weakly) stationary: constant mean, constant variance, autocorrelation depending only on the lag. Real economic and demand series are not stationary, and the practical workflow is to difference, log-transform, or detrend before computing the ACF — but only after fitting the transform on the training window and applying it unchanged to the validation and test windows, the same discipline as for cross-sectional feature engineering.

Forecasting metrics: MAPE, RMSE, MASE

Three metrics dominate the forecasting literature and they are not interchangeable.

Mean Absolute Percentage Error is

It has three problems. First, it is undefined when $y_i = 0$ and unstable when $y_i$ is small — a 1-unit error on a series of value 1 is 100% MAPE, while the same error on a series of value 100 is 1%. Second, it is asymmetric: a forecast that overshoots by 50% has a smaller MAPE than a forecast that undershoots by 50% of the same absolute amount, because the denominator is the realised value not the forecast. Third, it weights low-volume rows more heavily than high-volume rows, which biases model selection toward series with smaller denominators.

Root Mean Squared Error is

It penalises large errors more than small ones because of the squaring, which is appropriate when the cost of a large miss is itself large. It is on the same scale as $y$, so it is interpretable. But it is not scale-free: a model that produces a 5-unit RMSE on a series with mean 1000 is excellent, while a 5-unit RMSE on a series with mean 10 is terrible.

Mean Absolute Scaled Error (Hyndman & Koehler, 2006) is

The denominator is the mean absolute error of a seasonal naïve forecast with season $m$: predict $y_t$ by repeating the value from $m$ periods ago. MASE $\lt 1$ means the model beats seasonal naïve; MASE $\ge 1$ means it does not. The metric is symmetric, well-defined at zero (the denominator is non-zero for any non-constant seasonal process), and scale-free in the sense that it is invariant to linear rescaling of the series.

The 201-level choice is MASE on a competitive forecasting task, RMSE when the cost function is explicitly quadratic, and MAPE only when the percentage framing is required by stakeholders and the series is bounded away from zero. Reporting only MAPE is a tell that the practitioner has not thought about the metric's bias.

Leakage audit on a forecasting pipeline

Every forecasting pipeline is a chain of operations that each must be checked for time-leak. The audit has six questions, in order.

1. What is the row grain? A row at time $t$ must contain only information available at or before $t$. If the row aggregates data over a window that includes the future, every row leaks. The audit opens with the row's as-of date and walks through the construction.
2. What is the feature lookback? Each feature must use data from times before the row's as-of date. A "last 7 days average" feature computed in a way that includes the current day is a 1-step-ahead leak. The fix is to lag every aggregate by one period.
3. Are transforms fit on the training window only? A scaler fit on the full series leaks future summary statistics into the training data. The fix is to fit on the training window and apply unchanged to validation and test, exactly as in cross-sectional feature engineering.
4. Does the model use future context? Any model that uses bidirectional context (a Transformer, a bidirectional RNN, a CNN with future-frames padding) is leaking when deployed in a forecasting setting. The fix is to constrain the model to causal context only.
5. Is the horizon honest? If the model claims to forecast horizon $h$ but uses any feature that requires data closer than $h$ periods ahead, the effective horizon is shorter than claimed. The fix is to walk every feature back to its minimum lookback.
6. Is the validation order-preserving? The validation must use a temporal split, not a random split. This is the lesson's central point, but it is also the easiest check to forget in a notebook.

The order matters. Question 1 is the most likely place for a silent leak to enter, because the row construction is often a one-line SQL query that does not announce its as-of date. Question 4 is the most expensive to retrofit, because converting a bidirectional model to a causal one changes its capacity and often its accuracy.

Baselines and when they win

A learned forecasting model is only justified when it beats a baseline. The two baselines that cover most practical settings are seasonal naïve (predict $y_t$ by repeating $y_{t-m}$) and drift (predict $y_t$ by extrapolating the average slope from the start of the series to the most recent observation).

The drift baseline is

It is the simplest model of a non-stationary series with a linear trend. The seasonal-naïve baseline is

It is the simplest model of a series with a known seasonal period $m$. Either baseline can be competitive with a sophisticated deep model on series where the seasonal or trend structure dominates the dynamics.

The 201-level judgement is that a learned model must beat the appropriate baseline on the appropriate metric, on a walk-forward validation, before it is shipped. The "appropriate baseline" is not always seasonal naïve — for a series with a strong trend, drift is the correct baseline. The "appropriate metric" is MASE or RMSE, not MAPE. The "walk-forward validation" is not optional. A model that beats seasonal naïve on random $k$-fold but loses on walk-forward is a model whose apparent skill is leakage, not forecasting ability.

import numpy as np

def seasonal_naive(y: np.ndarray, horizon: int, season: int) -> np.ndarray:
    """
    Predict y_{t+h} = y_{t+h-m*ceil(h/m)} for h = 1..horizon.

    The baseline is the most recent value from the same season; it wins on
    strongly seasonal series with no trend. Anything a learned model beats
    this on is justified skill.
    """
    out = np.empty(horizon)
    for h in range(1, horizon + 1):
        out[h - 1] = y[-(season * ((h - 1) // season + 1) - (season - (h % season)) % season) - 1] \
            if (h % season) != 0 else y[-season * (h // season)]
    return out

The function is intentionally verbose: the indexing for the seasonal-naïve lookup is the kind of detail that is silently wrong in clean-looking vectorised code. The point is that the baseline must be correctly implemented before it is used as a benchmark; a buggy baseline that under-forecasts is not a real baseline.

Stationarity testing and the Dickey-Fuller derivation

Walk-forward validation is necessary but not sufficient. If the series is non-stationary — its mean or variance drifts over time — the validation score is comparing a model fit on one distribution to a model evaluated on a different distribution. The mean of $y_t$ might be drifting up while the variance is constant; a model that learned the level of the training window predicts at the wrong level on the validation window, not because the model is bad but because the world moved.

The Augmented Dickey-Fuller (ADF) test is the formal check. The null hypothesis is that the series has a unit root — i.e. $y_t = y_{t-1} + \varepsilon_t$ for some innovation $\varepsilon_t$. The test statistic comes from the regression

and the test is on $\gamma$. If $\gamma < 0$ and significantly so (the test statistic is compared to a non-standard distribution tabulated by Dickey and Fuller), the null of a unit root is rejected and the series is stationary around a linear trend.

The intuition worth holding onto: $\gamma < 0$ is what makes $y_{t-1}$ a useful predictor of $\Delta y_t$, which is what makes $y_t$ mean-reverting. If $\gamma = 0$, then $y_{t-1}$ carries no information about the next change and the series is a random walk — the variance grows without bound and there is no "level" to forecast. The ADF's non-standard distribution arises because under the null the regression is a spurious regression: $t$-statistics on $y_{t-1}$ do not have the standard Gaussian limit, and the Dickey-Fuller critical values are larger in magnitude than the usual 5% threshold of $-1.96$.

import numpy as np
from statsmodels.tsa.stattools import adfuller

def adf_pvalue(y: np.ndarray, max_lag: int | None = None) -> float:
    """
    Return the ADF p-value for the null of a unit root.

    A low p-value (e.g. < 0.05) rejects the null and supports stationarity.
    A high p-value means the series is consistent with a random walk under
    the test's assumptions.
    """
    result = adfuller(y, maxlag=max_lag or 0, autolag="AIC")
    return float(result[1])

The failure mode of the ADF test is low power against alternatives that are "almost" unit-root processes — series with $\gamma$ close to zero but not exactly zero. The KPSS test inverts the null (stationarity is the null, unit root is the alternative) and a series that passes ADF and fails KPSS is one where the evidence is genuinely mixed; the 201-level practice is to report both and treat the conclusion as "probably stationary" rather than "stationary".

ARIMA and the role of integration

Once the series is differenced to stationarity, the natural model is the ARIMA($p, d, q$) — $p$ autoregressive lags, $d$ differences, $q$ moving-average lags. The model is

where $L$ is the lag operator ($L y_t = y_{t-1}$), $(1-L)^d y_t$ is the $d$-th difference, and $\varepsilon_t$ is white noise. The $p$ and $q$ are chosen by minimising an information criterion (AIC or BIC) over a small grid.

The integration order $d$ is what ARIMA owes to the unit-root pre-processing: $d = 0$ means the original series is stationary, $d = 1$ means the first difference is stationary, and so on. The estimator combines Yule-Walker equations for the AR part, an innovation algorithm for the MA part, and a Kalman filter for the joint estimation. The standard errors are computed from the Hessian of the log-likelihood at the optimum.

The 201-level judgement: ARIMA is the right baseline for univariate forecasting on series with tens to a few hundred observations, where a learned neural model is over-parameterised and a non-parametric kernel smoother is too data-hungry. The Box-Jenkins methodology — identify, estimate, diagnose — is a discipline, not an algorithm: identification is from the ACF and PACF after stationarising; estimation is by maximum likelihood; diagnosis is by residual plots and the Ljung-Box test for residual autocorrelation. A model with autocorrelated residuals is mis-specified; the fix is more lags, a different $d$, or a transformation.

Autocorrelation-robust inference: Newey-West standard errors

The OLS standard error assumes independent residuals. Under autocorrelated residuals the OLS coefficient is still unbiased, but the variance estimate is wrong:

with $\sigma^2 = \frac{1}{n}\sum_i \hat{\varepsilon}_i^2$ underestimates the true variance when the residuals are positively autocorrelated. The Newey-West (1987) estimator replaces $\sigma^2$ with a long-run variance estimate that sums the autocovariances up to a bandwidth $L$:

where $\hat{\Omega}_L = \hat{\Gamma}_0 + \sum_{k=1}^{L}\bigl(1 - \frac{k}{L+1}\bigr)\bigl(\hat{\Gamma}_k + \hat{\Gamma}_k^{\top}\bigr)$ and $\hat{\Gamma}_k = \sum_{t=k+1}^{n}\hat{\varepsilon}_t\hat{\varepsilon}_{t-k}\,x_t x_{t-k}^{\top}$.

The bandwidth $L$ controls the bias-variance trade-off. Small $L$ ignores long-range dependence and the estimator is biased downward; large $L$ sums noisy autocovariance estimates and the variance blows up. The default in practice is $L = \lfloor 4(n/100)^{2/9}\rfloor$, the Newey-West rule of thumb.

The 201-level practice: any time-series regression whose conclusions are reported with OLS standard errors is over-confident. The fix is to always use Newey-West (or a HAC variant) when the regression's residuals show significant autocorrelation. The hypothesis test is the same — the $t$-statistic is the coefficient divided by the Newey-West standard error — but the standard error is honest.

Spectral analysis and the periodogram

For a stationary series, the ACF and the spectral density are Fourier pairs:

The spectral density $f(\omega)$ says how much of the series's variance is concentrated at frequency $\omega$; the periodogram is the sample estimate $\hat{f}(\omega) = \frac{1}{n}\lvert \sum_t y_t e^{-i\omega t}\rvert^2$.

A spike in the periodogram at frequency $\omega_0$ is the signature of a periodic component at period $T = 2\pi / \omega_0$. A daily demand series has spikes at the daily frequency and its harmonics; an economic series has spikes at the quarterly frequency; a meteorological series has spikes at the annual frequency. The periodogram is the diagnostic plot for detecting seasonality that the ACF/PACF pair cannot easily see: the seasonal lag might be hundreds of periods and the ACF at lag 100 is noisy.

The failure mode is that the periodogram is a biased estimator of the spectral density (it is inconsistent) and the spikes leak across frequencies (spectral leakage). The fix is a windowed periodogram (multiply the series by a taper before transforming) and Welch's method (average periodograms across overlapping windows). The 201-level practice is to use Welch's periodogram as a diagnostic, then test whether the dominant frequencies are stationary across the window.

Worked example: monthly retail demand, end to end

The following walks the full pipeline on a synthetic monthly demand series with trend, seasonality and noise. It exercises stationarity testing, ARIMA fitting, walk-forward validation and the MASE benchmark.

import numpy as np
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error

rng = np.random.default_rng(0)
T = 120  # 10 years of monthly data
t = np.arange(T)
trend = 100 + 0.5 * t
season = 10 * np.sin(2 * np.pi * t / 12)
noise = rng.normal(scale=3.0, size=T)
y = trend + season + noise

# Step 1: stationarity
p_orig = adfuller(y, autolag="AIC")[1]
p_diff = adfuller(np.diff(y), autolag="AIC")[1]
print(f"ADF p-value (raw):     {p_orig:.3f}")
print(f"ADF p-value (diff=1):  {p_diff:.3f}")

# Step 2: seasonal-naive baseline
seasonal_naive_pred = y[12:T]  # shift by 12 months
seasonal_naive_mase = np.abs(y[12:] - seasonal_naive_pred).mean()

# Step 3: walk-forward ARIMA
errs = []
for i in range(60, T - 12):
    train = y[:i]
    model = ARIMA(train, order=(1, 1, 1),
                  seasonal_order=(0, 0, 0, 0)).fit()
    fc = model.forecast(steps=12)
    errs.append(np.abs(fc - y[i: i + 12]).mean())

arima_mase = np.mean(errs) / seasonal_naive_mase
print(f"ARIMA walk-forward MASE: {arima_mase:.3f}")

The printout is the lesson's argument in three numbers. The raw ADF p-value is large (the series is non-stationary — it has a trend); the differenced p-value is small (the first difference is stationary — the trend has been removed by $1-L$). The walk-forward ARIMA achieves a MASE well below 1, meaning it beats the seasonal-naïve baseline by a margin that is meaningful under the framework. The 201-level judgement is that the comparison is honest precisely because no information from the future leaked into the training window at any step; the seasonal-naïve baseline is correctly implemented; and the metric is MASE, not MAPE.

The 201-level mistakes the example guards against: computing MAPE on a series that includes near-zero values (the metric blows up); fitting ARIMA on the full series before walk-forward (the validation is contaminated); and reporting a single number without the per-step variance. Walk-forward validation is a Monte Carlo estimator of the deployment score; its variance is the variance across steps, and a one-number summary that ignores it is over-confident.

Key Takeaways

• Random $k$-fold cross-validation on time-ordered data silently produces optimistic scores by leaking nearby time points across folds; walk-forward validation is the only honest evaluation.
• The ACF and PACF are diagnostic plots that require stationarity; differencing or detrending is the fix, and Newey-West standard errors are required for any regression on autocorrelated residuals.
• MAPE is asymmetric, undefined at zero and biased toward low-volume series; MASE benchmarks against the seasonal-naïve baseline and is the right default for cross-series comparison.
• A forecasting pipeline must be audited at the row grain, the feature lookback, the transform-fit discipline, the model causality, the horizon claim, and the validation ordering — every step is a possible leak.
• A learned model must beat the seasonal-naïve or drift baseline on walk-forward MASE before it is shipped; beating the baseline on random $k$-fold MAPE is not a justification, it is leakage.
• ARIMA, ADF, Newey-West and the periodogram together cover the canonical univariate time-series toolkit; the same walk-forward discipline applies to every one of them, and the only honest summary is MASE relative to a correctly-implemented baseline.