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.7339

Ridge’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.

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.3937

The 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. 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.

References

• Hyndman and Athanasopoulos. Forecasting: Principles and Practice: Evaluating point forecast accuracy defines scaled errors.
• Hyndman and Athanasopoulos. Forecasting: Principles and Practice: Time series cross-validation explains evaluation at rolling origins.
• Hyndman and Athanasopoulos. Forecasting: Principles and Practice: ARIMA models develops the modeling alternative.
• Hyndman and Athanasopoulos. Forecasting: Principles and Practice: Prediction intervals explains interval construction.
• scikit-learn documentation. TimeSeriesSplit describes its row-based split and gap parameters.