Time-Series Forecasting
A forecast uses information available at a particular time to predict a later observation. For hourly demand, that might mean using readings through 10:00 to predict 11:00, or issuing the next twenty-four hourly predictions together. Those tasks allow different inputs. A useful backtest reproduces when observations arrive, when the model is refitted, and how far ahead its predictions must reach.
Turn a history into prediction rows
Write \(y_t\) for the observation at time \(t\). A forecast origin is the time at which the prediction is issued; horizon \(h\) is the number of steps ahead. After observing \(y_t\), a one-step forecast targets \(y_{t+1}\). In code below, a row indexed by its target time \(i\) uses lag 1, \(y_{i-1}\), and lag 24, \(y_{i-24}\), among its features. If the history is [10, 12, 11], the next row starts with lag 1 = 11 and lag 2 = 12. The value being predicted is excluded.
These examples use equally spaced observations available immediately after measurement. With hourly data, lag 24 is one day; an omitted hour must not silently turn it into twenty-four irregular records. In real data, distinguish event time from publication time. Calendar variables and committed schedules can be known in advance; future measured weather or revised sales totals usually are not. Features must use the version available when the forecast would have been issued.
Start with simple forecasts and meaningful error units
The synthetic series combines a linear trend, a repeating twenty-four-step sine wave, accumulated random increments, and observation noise. The accumulated increments form a random walk. Repeating the most recent value is the naive forecast; repeating the value one period earlier is seasonal naive. Both are useful comparisons, but their strength depends on the series and horizon.
The code makes a chronological split and fits each model once. During evaluation, each new one-step prediction receives the actual observations that have arrived since the split. It does not issue the entire test sequence at the initial cutoff. Run the blocks in order with NumPy and scikit-learn installed.
Mean absolute error (MAE) retains the target’s units, such as units sold per hour. The test/naive column divides model MAE by naive MAE on the same test rows. MASE, the mean absolute scaled error, instead uses a training-history scale. With training length \(T\) and seasonal period \(s\), that scale is \(a_s=(T-s)^{-1}\sum_{t=s}^{T-1}|y_t-y_{t-s}|\), using zero-based indices; MASE is test MAE divided by \(a_s\). The example uses \(s=24\). For nonseasonal scaling use \(s=1\), and report which choice you used. A zero scale makes the ratio undefined.
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_absolute_error
def series(n=2000, seed=0, period=24):
g = np.random.default_rng(seed); t = np.arange(n)
return (0.01 * t # trend
+ 3.0 * np.sin(2 * np.pi * t / period) # seasonality
+ np.cumsum(g.normal(0, 0.25, n)) # random-walk drift
+ g.normal(0, 0.6, n)) # observation noise
y = series()
lags = (1, 2, 3, 24, 48); m = max(lags)
X = np.column_stack([y[m - l: len(y) - l] for l in lags])
target = y[m:]
split = int(0.8 * len(target))
Xtr, Xte, ytr, yte = X[:split], X[split:], target[:split], target[split:]
naive = Xte[:, 0] # lag 1: repeat the last value
seasonal = Xte[:, 3] # lag 24: repeat one period ago
naive_mae = mean_absolute_error(yte, naive)
cut = m + split
training_history = y[:cut]
seasonal_scale = np.mean(np.abs(training_history[24:] - training_history[:-24]))
print(f"training seasonal scale {seasonal_scale:.4f}")
print(f"naive (last value) MAE {naive_mae:.4f}")
print(f"seasonal naive (lag 24) MAE {mean_absolute_error(yte, seasonal):.4f}")
for name, model in (("Ridge on lags", Ridge()),
("HistGradientBoosting", HistGradientBoostingRegressor(random_state=0))):
mae = mean_absolute_error(yte, model.fit(Xtr, ytr).predict(Xte))
print(f"{name:26s} MAE {mae:.4f} test/naive {mae / naive_mae:.4f} seasonal MASE {mae / seasonal_scale:.4f}")
# training seasonal scale 1.1709
# naive (last value) MAE 0.8186
# seasonal naive (lag 24) MAE 1.2146
# Ridge on lags MAE 0.7397 test/naive 0.9037 seasonal MASE 0.6317
# HistGradientBoosting MAE 0.8593 test/naive 1.0498 seasonal MASE 0.7339Ridge’s test MAE is about 0.740, compared with 0.819 for the last-value forecast, a reduction of about 9.6% on this split. The gradient booster scores about 0.859 and loses this comparison. That measures accuracy under these settings; it does not establish a universal model ranking or the model’s business value. The random-walk increments are unpredictable in this generator, while its trend and periodic component offer structure to estimate. The experiment does not establish how close Ridge is to the best attainable forecast.
Seasonal naive scores about 1.215 here. It cancels the deterministic seasonal component at a twenty-four-step lag but accumulates more random-walk uncertainty and trend change than the one-step baseline. Comparing the two errors alone does not isolate which component caused their difference.
Thus MASE below one means the test MAE is below the training seasonal-difference scale, not necessarily below a competing forecast on the test interval. Keep raw MAE alongside scaled errors and actual test baselines. For Ridge, 0.7397 divided by the training scale 1.1709 gives seasonal MASE 0.6317, while dividing by test naive MAE 0.8186 gives 0.9037. Scaling helps comparisons across units; units and operational costs still matter.
Backtesting must match the forecast schedule
Evaluate at several historical origins, fitting only on observations and labels available by each origin. An expanding window retains the available history; a sliding window retains a chosen recent portion. Each validation block can represent one forecast for the whole block or a sequence of one-step forecasts updated with newly observed inputs. State which one you use. Exercise 2 uses the latter, refitting the model every 120 observations.
Choose lags, model settings, and window lengths on earlier rolling validation periods, then evaluate the selected procedure on a later untouched interval. The examples compare fixed settings to illustrate behavior; their test scores should not become a repeatedly consulted tuning target. Nearby errors and overlapping multi-step forecasts are dependent, so their spread is not an independent-sample confidence interval.
TimeSeriesSplit does not protect the features
A time split controls the training and evaluation rows. A feature can still contain later observations. The next example compares past-only rolling averages with a centered average that includes the target and future values. All methods use the same retained rows and chronological folds; each fold scores successive one-step predictions.
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.model_selection import TimeSeriesSplit, cross_val_score
# a slow signal buried in heavy independent noise, so smoothing is informative
def series2(n=3000, seed=0):
g = np.random.default_rng(seed); t = np.arange(n)
signal = 2 * np.sin(2 * np.pi * t / 200) + 1.2 * np.sin(2 * np.pi * t / 77)
return signal + g.normal(0, 3.0, n)
def roll(v, w, centred):
out = np.full(len(v), np.nan); h = w // 2
for i in range(len(v)):
if centred:
if i - h >= 0 and i + h < len(v): out[i] = v[i - h:i + h + 1].mean()
elif i - w >= 0: out[i] = v[i - w:i].mean()
return out
y2 = series2(); lags = (1, 2, 3, 5, 10); m = max(lags)
base = np.column_stack([y2[m - l: len(y2) - l] for l in lags]); tgt = y2[m:]
cen, cau = roll(y2, 41, True)[m:], roll(y2, 41, False)[m:]
ok = ~(np.isnan(cen) | np.isnan(cau))
base, tgt, cen, cau = base[ok], tgt[ok], cen[ok], cau[ok]
for name, Xf in (("lags only", base),
("+ trailing rolling mean (past only)", np.c_[base, cau]),
("+ centred rolling mean (uses future)", np.c_[base, cen])):
s = cross_val_score(Ridge(), Xf, tgt, cv=TimeSeriesSplit(5),
scoring="neg_mean_absolute_error")
print(f"{name:38s} MAE {-s.mean():.4f}")
print(f"oracle expected MAE under Gaussian noise {3.0 * np.sqrt(2 / np.pi):.4f}")
# lags only MAE 2.5268
# + trailing rolling mean (past only) MAE 2.5239
# + centred rolling mean (uses future) MAE 2.3507
# oracle expected MAE under Gaussian noise 2.3937The forty-one-value centered mean uses twenty observations before the target, the target itself, and twenty after it. Its reported gain over the past-only feature is about 0.17 MAE. Training-row order is chronological, but the constructed training features can cross a fold boundary, and validation features contain their own answers. This is leakage regardless of whether the resulting score looks plausible.
For this generator, an oracle knowing the signal has expected absolute error \(3\sqrt{2/\pi}\approx2.394\), because the remaining noise is independent and Gaussian. A finite test average can fall below that expectation by chance, including the oracle’s own score. The number is a population benchmark, not a hard bound on every observed MAE. The evidence of leakage here is the feature’s access to unavailable observations; a below-benchmark score is only a reason to investigate.
Check feature availability directly. Past-only rolling calculations can be prepared over a full array if each row uses only its available history. Scalers, imputation parameters, and other learned transformations must be fitted on the training portion. Centered filters and backward filling can import unavailable future values. Ordinary forward filling from an already available observation is a different operation, although stale values may still be unsuitable.
A future target is normal in forecasting. The training label must have become observable before fitting, however: a target summarizing the next seven days is not available on its row’s start date. Exclude such unresolved labels at each cutoff. TimeSeriesSplit’s gap can leave rows between training and evaluation, but the required separation depends on label windows and reporting delays. It cannot repair a feature that reads its own target. For regular samples it also makes equal-sized folds represent comparable durations; irregular timestamps need explicit time boundaries.
Choose what a multi-step forecast produces
A recursive forecast repeatedly uses a one-step model and feeds its predictions back into the lag history. At step two, lag 1 is the first prediction, while older lags may still be observed values. A direct forecast fits a separate model for each horizon using the history available at issuance. A multi-output model predicts a whole vector of horizons jointly; hybrids also exist.
Recursion can propagate errors and encounters predicted inputs at deployment after training on observed lags. Direct models avoid that particular feedback but require enough training examples for each horizon. Intermediate predictions are functions of the existing history, not newly observed information that direct forecasting necessarily loses. Their usefulness depends on the models and series.
For a fair comparison, use the same forecast origins and target timestamps and ensure every training target precedes the fitting cutoff. Exercise 1 keeps model fits fixed, then advances through common origins as actual history becomes available. Within each forty-eight-step forecast, recursion receives no later actual values. Baseline errors are measured separately at each horizon. Standard MASE can retain its declared training scale; it need not be redefined as a horizon-specific test-baseline ratio.
Model choices and uncertainty
Lag regression is one route. Exponential smoothing updates estimates of level and, in suitable variants, trend and seasonality. ARIMA models combine dependence on past values, differencing where needed, and moving-average terms in past innovations (one-step prediction errors). These moving-average terms differ from smoothing the observed values with a rolling mean. Differencing replaces levels with changes, for example y[t] − y[t−1]; it can remove some trends but does not guarantee a stable process. Boosted trees can combine lags with calendar and external predictors, while sequence models can learn across many related histories. Compare the methods under the same availability and horizon rules rather than choosing by complexity alone.
A point forecast gives one value. For decisions such as inventory planning, prediction intervals or quantiles can express uncertainty in future observations. For example, coverage of a nominal 90% interval is the fraction of evaluated observations inside it, with 90% as the target. Evaluate coverage and width at each relevant horizon on rolling held-out periods; a wider interval is easier to cover but less precise. Quantile loss can train asymmetric point targets, such as a demand level intended to exceed actual demand 90% of the time under a calibrated 90th-percentile forecast. Independent resampling of time-series rows can destroy dependence, and a fixed residual interval may lose coverage after a regime change. Forecast uncertainty generally includes future variation as well as estimation error; it is not just uncertainty about a fitted mean.
Exercises
1. Compare recursive and direct forecasts at common origins.
Solution
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_absolute_error
def forecast_series(n=2000, seed=0, period=24):
g = np.random.default_rng(seed)
t = np.arange(n)
return (0.01 * t + 3 * np.sin(2 * np.pi * t / period)
+ np.cumsum(g.normal(0, 0.25, n)) + g.normal(0, 0.6, n))
yh = forecast_series()
horizon_lags = (1, 2, 3, 24, 48)
max_lag, H = max(horizon_lags), 48
cutoff = max_lag + int(0.8 * (len(yh) - max_lag))
def lag_rows(values, first_targets):
return np.array([[values[o-lag] for lag in horizon_lags] for o in first_targets])
train_one = np.arange(max_lag, cutoff)
one_step = Ridge().fit(lag_rows(yh, train_one), yh[train_one])
# o is the first unknown timestamp: history is available through o-1.
origins = np.arange(cutoff, len(yh) - H + 1, 7)
def recursive(model, history, steps):
history = list(history)
predictions = []
for _ in range(steps):
row = np.array([[history[-lag] for lag in horizon_lags]])
prediction = model.predict(row)[0]
history.append(prediction)
predictions.append(prediction)
return np.array(predictions)
recursive_predictions = np.array([recursive(one_step, yh[:o], H) for o in origins])
print("horizon recursive direct naive seasonal")
for h in (1, 6, 12, 24, 48):
train_origins = np.arange(max_lag, cutoff-h+1)
target_times = train_origins + h - 1
assert target_times.max() < cutoff
direct_model = Ridge().fit(lag_rows(yh, train_origins), yh[target_times])
actual = yh[origins+h-1]
direct = direct_model.predict(lag_rows(yh, origins))
naive = yh[origins-1]
seasonal = yh[origins+h-1 - 24*((h+23)//24)]
errors = [mean_absolute_error(actual, pred) for pred in
(recursive_predictions[:, h-1], direct, naive, seasonal)]
print(h, *(round(value, 4) for value in errors))
# horizon recursive direct naive seasonal
# 1 0.6964 0.6964 0.8167 1.13
# 6 1.3386 2.5472 2.6884 1.1167
# 12 1.487 4.7181 3.8809 1.2633
# 24 1.5644 1.5201 1.523 1.523
# 48 2.25 1.8007 1.7674 1.7674The training cutoff is shared. A direct h-step training row is retained only if its target at o+h−1 is before that cutoff. All four methods use the same origins for every displayed horizon, with enough future observations to evaluate forty-eight steps. At horizon one, the direct and recursive fits and predictions coincide, which is a useful alignment check.
On these origins, recursive MAE is 1.4870 against direct MAE 4.7181 at twelve steps. At forty-eight steps direct improves on recursive, 1.8007 against 2.2500, yet the naive baseline scores 1.7674. Comparing the learned models alone would miss that baseline result.
At a twenty-four-step horizon, the deterministic seasonal component returns to the same phase. At twelve steps it changes sign. This explains the seasonal contribution to the last-value baseline’s non-monotone horizon errors; the amount also depends on phase, drift, and noise. It does not make twenty-four-step forecasting universally easier. Seasonal naive repeats the last available cycle, including for horizons longer than one period.
Evaluate the horizon you intend to serve. Good one-step performance alone does not establish good twelve-step performance. Comparing more model settings on this table would turn it into validation data and require a later final test.
2. Compare expanding and sliding windows around a change.
Solution
import numpy as np
from sklearn.linear_model import Ridge
def changing_series(n=3000, seed=0, period=24, regime_at=2000):
g = np.random.default_rng(seed)
t = np.arange(n)
amplitude = np.where(t < regime_at, 3.0, 1.0)
return (0.01 * t + amplitude * np.sin(2 * np.pi * t / period)
+ np.cumsum(g.normal(0, 0.25, n)) + g.normal(0, 0.6, n))
yw = changing_series()
window_lags = (1, 2, 3, 24, 48)
mw = max(window_lags)
features = np.column_stack([yw[mw-lag:len(yw)-lag] for lag in window_lags])
response = yw[mw:]
step, window = 120, 600
point_errors, block_errors = [], []
for o in range(int(0.4 * len(response)), len(response), step):
end = min(o+step, len(response))
expanding = Ridge().fit(features[:o], response[:o])
sliding = Ridge().fit(features[max(0, o-window):o], response[max(0, o-window):o])
exp_error = np.abs(response[o:end] - expanding.predict(features[o:end]))
sli_error = np.abs(response[o:end] - sliding.predict(features[o:end]))
times = np.arange(o+mw, end+mw)
point_errors.extend(zip(times, exp_error, sli_error))
block_errors.append(exp_error.mean())
errors = np.array(point_errors)
for name, mask in (("before", errors[:, 0] < 2000),
("after", errors[:, 0] >= 2000),
("overall", np.ones(len(errors), dtype=bool))):
print(name, "n", int(mask.sum()), "expanding", round(errors[mask, 1].mean(), 4),
"sliding", round(errors[mask, 2].mean(), 4))
print("block MAE range", round(min(block_errors), 4), round(max(block_errors), 4))
# before n 772 expanding 0.734 sliding 0.7412
# after n 1000 expanding 0.6861 sliding 0.6569
# overall n 1772 expanding 0.707 sliding 0.6936
# block MAE range 0.6314 0.8016The seasonal amplitude changes at timestamp 2000. Predictions use newly observed lags within each block, with model coefficients held fixed until the next refit. They are successive one-step forecasts, not 120-step forecasts from the block’s start.
Before/after groups use each target’s timestamp. A block crossing the change contributes some predictions to each group; classifying the whole block by its start would mislabel those cases. The overall result weights each predicted observation equally, including the shorter final block.
Before the change, expanding-window MAE is 0.7340 against sliding-window MAE 0.7412. After it, the comparison is 0.6861 against 0.6569. These averages show a reversal in the observed ranking, not a test of statistical significance.
A sliding window reduces the influence of older data; an expanding window can estimate stable relationships from more observations. The error comparison does not separately measure bias or variance, and a single amplitude change does not make either model’s advantage inevitable. A window can help while it still contains some old-regime observations. Validate its length on earlier cutoffs rather than assuming it must be shorter than the time since a break.
The range summarizes errors from these evaluation blocks. Changes in fitted history and in the evaluated observations both contribute, so it is not the uncertainty of one fixed test score. Inspect errors over time and relevant regimes alongside the aggregate.
3. Catch a feature that depends on unavailable values.
Solution
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.model_selection import TimeSeriesSplit, cross_val_score
g = np.random.default_rng(0)
n = 4000
signal = np.sin(2 * np.pi * np.arange(n) / 150)
noise_sd = 2.0
ycheck = signal + g.normal(0, noise_sd, n)
check_lags = (1, 2, 3, 5, 10)
times = np.arange(max(check_lags), n-5)
past = np.array([[ycheck[t-lag] for lag in check_lags] for t in times])
def window_feature(values, t, centered):
return values[t-5:t+6].mean() if centered else values[t-11:t].mean()
straddle = np.array([window_feature(ycheck, t, True) for t in times])
target = ycheck[times]
# A feature for predicting y[t] must not depend on observations at t or later.
folds = list(TimeSeriesSplit(5).split(past))
for name, features in (("past lags only", past),
("+ straddling mean", np.c_[past, straddle])):
score = cross_val_score(Ridge(), features, target, cv=folds,
scoring="neg_mean_absolute_error")
print(name, "MAE", round(-score.mean(), 4))
print("oracle expected MAE", round(noise_sd*np.sqrt(2/np.pi), 4))
held_out_rows = np.concatenate([test for _, test in folds])
print("oracle sample MAE", round(np.abs(target[held_out_rows]
- signal[times[held_out_rows]]).mean(), 4))
t = 100
altered = ycheck.copy()
altered[t:] += 100
print("past feature unchanged", np.isclose(window_feature(ycheck, t, False),
window_feature(altered, t, False)))
print("centered feature unchanged", np.isclose(window_feature(ycheck, t, True),
window_feature(altered, t, True)))
# past lags only MAE 1.6608
# + straddling mean MAE 1.4937
# oracle expected MAE 1.5958
# oracle sample MAE 1.6075
# past feature unchanged True
# centered feature unchanged FalseThe centered eleven-point window includes y[t] and five later values. It gives the model some of the realized target noise; it does not let a linear model recover that noise exactly from the average. Both the expected oracle MAE and its finite-sample MAE are printed so they cannot be mistaken for the same quantity.
The final check changes all values at t and later while preserving the history through t−1. A valid past-only feature for predicting y[t] stays unchanged; the centered feature changes. This demonstrates the leak independently of its score. Passing one perturbation test is not a proof that a full pipeline is safe: also inspect timestamps, fitted preprocessing, data revisions, and label availability.
The formula σ√(2/π) requires centered Gaussian noise. Standard deviation alone does not determine expected absolute error. Under independent centered finite-variance noise, the known-signal squared-error benchmark is σ². Neither expectation is a lower bound on every realized sample average.
Real noise estimates require their own measurement model. For example, the variance of the difference of two independent measurements of the same unchanged quantity, with equal noise variance σ², is 2σ². Correlated sensors, changing signals, and unequal precision invalidate that shortcut. A fitted model’s residuals can also include model error, so they do not automatically reveal irreducible noise.
References
Forecasting: Principles and Practice — forecast accuracy defines scaled errors. Its time-series cross-validation chapter explains evaluation at rolling origins. The chapters on ARIMA and prediction intervals develop the modeling alternatives. The TimeSeriesSplit documentation describes its row-based split and gap parameters.
Discover more from Insightful Data Lab
Subscribe to get the latest posts sent to your email.