How far ahead is this good for? A pricing lead asked me exactly that in a review, about a forecast I had already signed off on. I had a cross validation score to show her. I did not have an answer to her question, because a single averaged number over shuffled folds cannot say anything at all about horizon, and horizon was the only thing she cared about. Rebuilding that evaluation as a walk forward backtest took most of a day. Once it ran, the answer came back as roughly four months of useful accuracy, against the twelve months the model was already being used for. Nothing about the model changed. Only the way it was measured changed, and that was enough to stop a planning cycle built on the wrong numbers.

Why a shuffled split lies about time ordered data

Random cross validation rests on one assumption: that every row is an independent draw from the same distribution, so which rows land in the test fold does not matter. Time series violate that assumption twice over. Neighbouring observations are correlated, so a randomly held out point almost always sits between two training points a few days either side, and predicting it is closer to interpolation than to forecasting. Beyond that, a time series usually has a trend, which means the future is not drawn from the same distribution as the past at all.

Rather than argue this, measure it. Below I use the Mauna Loa atmospheric carbon dioxide record that ships inside statsmodels, resampled to monthly averages. Churn does not have a natural time axis in the Telco snapshot we have been modelling, so this part borrows a clearly seasonal, clearly trending public series and ties the lesson back at the end. Lag features and a month indicator go in, a gradient boosted regressor comes out, and the identical model is scored two ways: once on a shuffled five fold split, once on the final ten years held out chronologically.

# tested against: Python 3.11, scikit-learn 1.7.2, statsmodels 0.14.6,
# pandas 2.3.3, numpy 2.2.6
import pandas as pd, statsmodels.api as sm
from sklearn.model_selection import KFold
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_absolute_error

raw = sm.datasets.co2.load_pandas().data
s = raw['co2'].resample('MS').mean().interpolate()

d = pd.DataFrame({'y': s})
d['t'] = range(len(d))
for L in (1, 2, 12):
    d[f'lag{L}'] = d['y'].shift(L)
d['month'] = d.index.month
d = d.dropna()

X = d.drop(columns='y')
y = d['y']
mk = lambda: HistGradientBoostingRegressor(random_state=0, max_iter=300)

tr, te = next(iter(KFold(5, shuffle=True, random_state=0).split(X)))
print('random shuffled split MAE  %.3f ppm' %
      mean_absolute_error(y.iloc[te], mk().fit(X.iloc[tr], y.iloc[tr]).predict(X.iloc[te])))

cut = len(X) - 120
print('last 10 years held out MAE %.3f ppm' %
      mean_absolute_error(y.iloc[cut:], mk().fit(X.iloc[:cut], y.iloc[:cut]).predict(X.iloc[cut:])))

random shuffled split MAE  0.364 ppm
last 10 years held out MAE 7.637 ppm

Twenty one times worse. Anybody who reported the first number and shipped against it would be roughly two decades of carbon dioxide growth wrong within a year. Two separate mechanisms produced that gap. Shuffling let the model see 1998 while predicting 1997, which is interpolation dressed up as forecasting. Worse, gradient boosted trees cannot extrapolate: a tree splits on thresholds it saw in training, so once the series climbs past the highest training value, every prediction is capped at the highest leaf average the model ever learned. Random folds hide that failure completely, because there is no such thing as a value beyond the training range when the training range covers the entire series.

Anatomy of an honest forecast split

Where the churn project stands: through Parts 9 to 15 we trained and tuned a churn classifier on the Telco snapshot, scored on a random stratified split, and through Part 17 we rebuilt it as a small network. Nothing about that evaluation was wrong, because the Telco file is a single snapshot with no observation dates in it. This part adds the axis that snapshot lacks. Once you start scoring churn monthly rather than once, every rule below starts applying to it too, and the stratified split you have been using quietly becomes the wrong tool.

scikit-learn ships the right primitive already. TimeSeriesSplit produces folds where the training set is always strictly earlier than the test set, and it grows the training window as it walks forward. Three arguments matter. n_splits sets how many folds you get and defaults to 5. test_size pins each test block to a fixed length rather than letting scikit-learn derive it, which is what you want when a fold should mean a quarter or a year. gap drops a number of observations between the end of training and the start of testing, and it is the single most underused argument in the whole module.

from sklearn.model_selection import TimeSeriesSplit

tscv = TimeSeriesSplit(n_splits=5, test_size=24, gap=6)
for i, (tr, te) in enumerate(tscv.split(X), 1):
    print('fold %d train %s..%s (n=%d)  test %s..%s' % (
        i, X.index[tr[0]].date(), X.index[tr[-1]].date(), len(tr),
        X.index[te[0]].date(), X.index[te[-1]].date()))

fold 1 train 1959-03-01..1990-12-01 (n=382)  test 1991-07-01..1993-06-01
fold 2 train 1959-03-01..1992-12-01 (n=406)  test 1993-07-01..1995-06-01
fold 3 train 1959-03-01..1994-12-01 (n=430)  test 1995-07-01..1997-06-01
fold 4 train 1959-03-01..1996-12-01 (n=454)  test 1997-07-01..1999-06-01
fold 5 train 1959-03-01..1998-12-01 (n=478)  test 1999-07-01..2001-06-01

Notice what the gap bought. Each training window ends in December and each test window opens the following July. Those six missing months are deliberate. If your label takes ninety days to confirm, then at the moment you score January you genuinely do not know October yet, and a fold that trains on October is scoring a model you could never have deployed. Setting gap to match your true label latency is the cheapest way to make a backtest honest, and it costs you nothing except a slightly smaller effective training set.

One decision remains: expanding window or rolling window. An expanding window keeps every historical observation, which is what TimeSeriesSplit does by default. A rolling window, set through max_train_size, keeps only the most recent block. Expanding is the better default because more data usually wins. Switch to rolling when the process itself changed, a pricing model rewritten in 2023 or a product relaunch, and old history is actively misleading rather than merely old.

Leakage that only exists in time series

Part 11 covered leakage in general terms. Time ordered data adds a family of leaks that a stratified split cannot even express, because they are all variations on one theme: a feature or a preprocessing step that quietly consumed information from after the prediction moment.

Four show up constantly. Fitting a scaler or an imputer on the full series hands the model the future mean and the future variance, and it is invisible because the numbers look perfectly ordinary. Centred rolling features are worse: a thirteen month centred moving average contains six months you will not have at prediction time, and it will look like your single best feature. Target derived aggregates such as a customer lifetime average computed over the whole file leak the outcome directly. Lastly, and most often missed, label availability lag: refunds settle, subscriptions confirm, and a churn flag stamped today may not be knowable for weeks.

Two defences cover most of this. Put every transformation inside a scikit-learn Pipeline so it refits within each fold rather than once over everything, a habit worth carrying from Part 11 into every forecasting project. Then write down, for each feature, the exact moment it becomes knowable. If a feature cannot state its own availability time, it does not go in the model. That rule has caught more bad features for me than any importance plot ever has, and it borrows directly from the discipline of exploratory data analysis, where you interrogate where a column came from before you trust what it says.

Walk forward backtesting with a moving origin

Cross validation folds answer a narrow question: how well does this model do on held out blocks. Forecasting asks something else. Standing at a given date with only history behind you, how good is a prediction twelve steps ahead. Walk forward backtesting, sometimes called evaluation on a rolling origin, answers that directly. Pick a series of origins, refit at each one using only data up to it, forecast the full horizon, then compare against what actually happened.

Below the same monthly series is backtested with SARIMAX, a seasonal autoregressive integrated moving average model with optional external regressors, fitted by maximum likelihood through a Kalman filter. Six origins, one per year from 1995, twelve months of horizon at each.

import numpy as np, pandas as pd, statsmodels.api as sm

H = 12
origins = pd.date_range('1995-01-01', '2000-12-01', freq='12MS')
err = {h: [] for h in range(1, H + 1)}

for o in origins:
    train = s[:o]
    model = sm.tsa.statespace.SARIMAX(
        train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12),
        enforce_stationarity=False).fit(disp=False)
    fc = model.forecast(H)
    actual = s.reindex(fc.index)
    for h, (f, a) in enumerate(zip(fc.values, actual.values), 1):
        if not np.isnan(a):
            err[h].append(abs(f - a))

for h in range(1, H + 1):
    print('%2d  %.3f   %d' % (h, np.mean(err[h]), len(err[h])))

Two failures are worth meeting deliberately, because both hit almost everyone once. Asking for more folds than the data supports raises a clear error, and a datetime index that lost its frequency raises a warning you should treat as an error.

>>> list(TimeSeriesSplit(n_splits=5, test_size=120).split(X))  # X has 514 rows
ValueError: Too many splits=5 for number of samples=514 with test_size=120 and gap=0.

>>> sm.tsa.statespace.SARIMAX(s2, order=(1,1,1), seasonal_order=(1,1,1,12)).fit(disp=False)
ValueWarning: No frequency information was provided, so inferred frequency MS will be used.

First one is arithmetic: five folds of 120 observations each plus a training set will not fit inside 514 rows, so either drop to three folds or shrink test_size. Second one is subtler and bites hardest in production. When a datetime index carries no frequency, statsmodels guesses, and here it guessed correctly. Guess wrongly and forecast returns timestamps at the wrong spacing, so your predictions silently align to the wrong months and every error you compute afterwards is meaningless. Call asfreq on the index yourself and never let the library infer it.

Accuracy decay across the forecast horizon

Here is what the backtest actually returned. Reading it as a single averaged number, 0.473 ppm, would throw away the only finding that matters.

Growth is not smooth, and that irregularity is real rather than noise to be sanded off. Horizon 8 is worse than horizon 9 because with only six origins each point averages six errors, so individual bad years show through. That is a genuine limitation of the setup and it belongs in whatever you report. More origins would settle the curve, at the cost of a refit each. Even at six origins the shape is unmistakable, and it converts directly into a policy anybody non technical can act on: this model is trustworthy to about four months, usable with caveats to eight, and beyond that you are quoting a trend line with decoration.

Statistical models against machine learning models

Section one showed gradient boosting collapsing at 7.637 ppm on a trending series while SARIMAX held below 0.72 ppm at a full year out. Do not read that as a general ranking. It is one specific weakness meeting one specific data shape. Trees cannot extrapolate beyond their training range, and this series does nothing but climb, so trees were the wrong instrument. Difference the target first, predicting month over month change instead of level, and boosting becomes competitive again, because the differenced series no longer trends.

One nuance deserves stating because the blanket rule against shuffling is not quite true. Bergmeir and Benitez showed that for a purely autoregressive model whose residuals are uncorrelated, standard K fold cross validation is in fact valid and often more efficient than a single holdout, since the dependence that breaks cross validation lives in the errors rather than in the ordering itself. Worth knowing, and worth being careful with. Conditions are strict, residual independence has to be tested rather than assumed, and the moment external regressors or non stationary structure enter, it stops holding. My position after several projects: use time ordered splits by default, and reach for shuffled folds only when you have checked the residuals and can defend the choice in writing.

Back to churn. Everything above lands on the project the moment it goes live, which it does in Part 25. A churn model scored once on a snapshot is a photograph. A churn model in production is a series, scored every month, and the honest question is not what its AUC was in July but whether July still tells you anything about November. Same distinction, same fix: score by origin, watch the curve, and do not average away the thing you needed to see. Correlation between a feature and churn is also not stability of that relationship over time, which is the same trap examined from a different angle in correlation versus causation.

Split protocol I recommend for forecasting work

Six steps, in order, and I have not found a forecasting project where skipping any of them paid off. Sort by time and confirm the index carries an explicit frequency. Fit a seasonal naive baseline and record its error, because a model that cannot beat last year same month is not a model. Build only trailing features and write down when each becomes knowable. Score with TimeSeriesSplit, setting gap to your real label latency and test_size to a business meaningful block. Backtest across at least six origins and report error per horizon rather than averaged. Then, and only then, compare model families.

If you take one habit from this part, take the fold boundary printout. Four lines of code, run before you look at any score, and it would have saved me those eleven days. Most bad forecasts I have reviewed were not caused by a weak model. They were caused by an evaluation that answered an easier question than the one being asked, and nobody noticed because the number looked good.

Run the backtest on whatever forecast you currently have in production, print the horizon curve, and see whether the horizon you are actually using is still inside the range the model can support. Most are not. Part 20 turns to recommender systems, where collaborative filtering brings its own version of this problem, since a user who has rated nothing at all cannot be split into any fold honestly.

References

• TimeSeriesSplit API reference, including gap and max_train_size, scikit-learn
• SARIMAX API reference, statsmodels 0.14.6
• Bergmeir and Benitez, On the use of cross validation for time series predictor evaluation, Information Sciences, 2012
• Bergmeir, Hyndman and Koo, A note on the validity of cross validation for evaluating autoregressive time series prediction
• Telco Customer Churn dataset, IBM repository

About The Author

Dr. Pranay Jha

Dr. Pranay Jha is a Cloud and AI Consultant with 18+ years of experience in hybrid cloud, virtualization, and enterprise infrastructure transformation. He specializes in VMware technologies, multi-cloud strategy, and Generative AI solutions. He holds a PhD in Computer Applications with research focused on Cloud and AI, has published multiple research papers, and has been a VMware vExpert since 2016 and a VMUG Community Leader.

See author's posts