A forecast is a range with a date on it. Most forecasting that goes wrong went wrong by producing a single number, by fitting the noise, or by never checking the error on data the model had not seen.
The three parts of most business series
import numpy as np
import pandas as pd
rng = np.random.default_rng(7)
weeks = pd.date_range('2024-01-07', periods=156, freq='W-SUN')
t = np.arange(len(weeks))
trend = 1000 + 7.5 * t
seasonal = 180 * np.sin(2 * np.pi * t / 52) + 60 * np.sin(2 * np.pi * t / 13)
noise = rng.normal(0, 45, len(weeks))
series = pd.Series(trend + seasonal + noise, index=weeks, name='revenue')
parts = pd.DataFrame({'trend': trend, 'seasonal': seasonal, 'noise': noise},
index=weeks)
print(parts.describe().loc[['mean', 'std', 'min', 'max']].round(1).to_string())
print()
print('share of total variation:')
total_variance = series.var()
for name in parts.columns:
print(f' {name:<9} {parts[name].var() / total_variance:.1%}')
Open this example in the editor and change it. Nothing you type is saved.
Almost every business series decomposes into those three: a trend that is the actual news, a seasonal pattern that repeats, and noise that is nothing at all. The work of forecasting is separating them, and the work of reporting is not mistaking the third for the first.
A quick visual decomposition, with no library:
import numpy as np
import pandas as pd
rng = np.random.default_rng(7)
weeks = pd.date_range('2024-01-07', periods=156, freq='W-SUN')
t = np.arange(len(weeks))
series = pd.Series(
1000 + 7.5 * t + 180 * np.sin(2 * np.pi * t / 52) + rng.normal(0, 45, len(weeks)),
index=weeks)
# Trend: a centred rolling mean over one full season.
trend = series.rolling(52, center=True, min_periods=26).mean()
detrended = series - trend
# Seasonal: the average detrended value for each week of the year.
seasonal = detrended.groupby(detrended.index.isocalendar().week).transform('mean')
residual = detrended - seasonal
print(f'series std {series.std():7.1f}')
print(f'trend std {trend.std():7.1f}')
print(f'seasonal std {seasonal.std():7.1f}')
print(f'residual std {residual.std():7.1f}')
print()
print('The residual is what no model can predict. If your forecast')
print('error is near it, you are done; if it is far above, there is')
print('structure left to find.')
Open this example in the editor and change it. Nothing you type is saved.
Start with the dumb baselines
Before any model, compute what three trivial rules would have predicted. A model that cannot beat them is not worth deploying.
import numpy as np
import pandas as pd
rng = np.random.default_rng(13)
weeks = pd.date_range('2024-01-07', periods=156, freq='W-SUN')
t = np.arange(len(weeks))
series = pd.Series(
1000 + 7.5 * t + 180 * np.sin(2 * np.pi * t / 52) + rng.normal(0, 45, len(weeks)),
index=weeks)
HORIZON = 26
train, test = series[:-HORIZON], series[-HORIZON:]
baselines = {
'last value': np.repeat(train.iloc[-1], HORIZON),
'mean of train': np.repeat(train.mean(), HORIZON),
'same week last year': series.shift(52)[-HORIZON:].to_numpy(),
'drift': train.iloc[-1] + (train.iloc[-1] - train.iloc[0]) / (len(train) - 1)
* np.arange(1, HORIZON + 1),
}
def mae(actual, predicted):
return float(np.mean(np.abs(np.asarray(actual) - np.asarray(predicted))))
def mape(actual, predicted):
actual = np.asarray(actual)
return float(np.mean(np.abs(actual - np.asarray(predicted)) / np.abs(actual)))
for name, prediction in baselines.items():
print(f'{name:<22} MAE {mae(test, prediction):7.1f} MAPE {mape(test, prediction):6.2%}')
Open this example in the editor and change it. Nothing you type is saved.
The four baselines, and what each is for:
- Last value (naive). The right baseline for anything close to a random walk - prices, stock levels.
- Seasonal naive - the same period last year. The right baseline for anything with a strong yearly pattern, and it is remarkably hard to beat.
- Mean. The right baseline when there is no trend at all.
- Drift. Last value plus the average change so far. A trend baseline that costs one line.
Publish the baseline's error beside your model's, always. "Our model has 6% error" is meaningless; "our model has 6% error against the seasonal naive's 9%" is a result.
Measure the error on data the model has not seen
import numpy as np
import pandas as pd
rng = np.random.default_rng(17)
weeks = pd.date_range('2024-01-07', periods=156, freq='W-SUN')
t = np.arange(len(weeks))
series = pd.Series(
1000 + 7.5 * t + 180 * np.sin(2 * np.pi * t / 52) + rng.normal(0, 45, len(weeks)),
index=weeks)
def fit_predict(train_values, horizon, degree):
"""A polynomial through the training data, extended forwards."""
x = np.arange(len(train_values))
coefficients = np.polyfit(x, train_values, degree)
future = np.arange(len(train_values), len(train_values) + horizon)
return np.polyval(coefficients, x), np.polyval(coefficients, future)
HORIZON = 26
train, test = series[:-HORIZON], series[-HORIZON:]
print(f'{"degree":>7} {"train MAE":>11} {"test MAE":>11}')
for degree in [1, 2, 3, 6, 12, 20]:
fitted, predicted = fit_predict(train.to_numpy(), HORIZON, degree)
train_error = float(np.mean(np.abs(train.to_numpy() - fitted)))
test_error = float(np.mean(np.abs(test.to_numpy() - predicted)))
print(f'{degree:>7} {train_error:>11.1f} {test_error:>11.1f}')
print()
print('Training error falls with every added degree. Test error does')
print('not, and past a point it explodes. That gap is overfitting.')
Open this example in the editor and change it. Nothing you type is saved.
A model is evaluated only on data it did not see. For a time series that means a split in time, never a random one - a random split lets the model see the future, and the error it reports is a fiction.
Rolling-origin evaluation is the honest version, because it tests the model at several points rather than one:
import numpy as np
import pandas as pd
rng = np.random.default_rng(19)
weeks = pd.date_range('2024-01-07', periods=156, freq='W-SUN')
t = np.arange(len(weeks))
series = pd.Series(
1000 + 7.5 * t + 180 * np.sin(2 * np.pi * t / 52) + rng.normal(0, 45, len(weeks)),
index=weeks)
def seasonal_naive(train_values, horizon, period=52):
tail = train_values[-period:]
return np.array([tail[i % period] for i in range(horizon)])
def drift(train_values, horizon):
step = (train_values[-1] - train_values[0]) / (len(train_values) - 1)
return train_values[-1] + step * np.arange(1, horizon + 1)
HORIZON = 13
errors = {'seasonal naive': [], 'drift': [], 'both averaged': []}
for origin in range(104, len(series) - HORIZON, 13):
train = series.iloc[:origin].to_numpy()
actual = series.iloc[origin:origin + HORIZON].to_numpy()
a = seasonal_naive(train, HORIZON)
b = drift(train, HORIZON)
errors['seasonal naive'].append(np.mean(np.abs(actual - a)))
errors['drift'].append(np.mean(np.abs(actual - b)))
errors['both averaged'].append(np.mean(np.abs(actual - (a + b) / 2)))
print(f'{len(errors["drift"])} origins, {HORIZON}-week horizon')
for name, values in errors.items():
print(f' {name:<16} MAE {np.mean(values):7.1f} (worst {np.max(values):7.1f})')
print()
print('Averaging two weak forecasts often beats either. That is the')
print('cheapest improvement available in forecasting.')
Open this example in the editor and change it. Nothing you type is saved.
The band is the output
import numpy as np
import pandas as pd
rng = np.random.default_rng(23)
weeks = pd.date_range('2024-01-07', periods=156, freq='W-SUN')
t = np.arange(len(weeks))
series = pd.Series(
1000 + 7.5 * t + 180 * np.sin(2 * np.pi * t / 52) + rng.normal(0, 45, len(weeks)),
index=weeks)
HORIZON = 13
train = series.iloc[:-HORIZON]
# Estimate the error at each horizon from history, then use it as the
# band. No distributional assumption beyond "the past is a guide".
residuals = {h: [] for h in range(1, HORIZON + 1)}
for origin in range(104, len(series) - HORIZON):
history = series.iloc[:origin].to_numpy()
step = (history[-1] - history[0]) / (len(history) - 1)
for h in range(1, HORIZON + 1):
predicted = history[-1] + step * h
residuals[h].append(series.iloc[origin + h - 1] - predicted)
step = (train.iloc[-1] - train.iloc[0]) / (len(train) - 1)
rows = []
for h in range(1, HORIZON + 1):
point = train.iloc[-1] + step * h
spread = np.percentile(residuals[h], [10, 90])
rows.append({'weeks ahead': h, 'point': point,
'low': point + spread[0], 'high': point + spread[1]})
band = pd.DataFrame(rows).round(0)
band['width'] = band['high'] - band['low']
print(band.to_string(index=False))
print()
print(f'the band widens from {band["width"].iloc[0]:.0f} at one week to '
f'{band["width"].iloc[-1]:.0f} at thirteen')
print('A forecast whose band does not widen with the horizon is lying.')
Open this example in the editor and change it. Nothing you type is saved.
Three things follow:
Always publish a range. A point forecast is a number somebody will treat as a commitment. A range is a forecast.
The range must widen with the horizon. Predicting next week is easier than predicting next quarter, and the picture should say so.
Say what the range means. An 80% interval means you expect to be outside it one time in five - which is often, and saying so in advance prevents the conversation where being outside it once is treated as failure.
When not to forecast
import pandas as pd
print(pd.DataFrame([
{'situation': 'fewer than two full seasonal cycles of history',
'why': 'you cannot separate seasonality from trend'},
{'situation': 'the thing being forecast is a decision',
'why': 'the forecast will change the decision and invalidate itself'},
{'situation': 'a structural break just happened',
'why': 'history is no longer a guide; the model will miss it entirely'},
{'situation': 'the series is dominated by a few huge events',
'why': 'you are forecasting whether an event occurs, not a level'},
{'situation': 'nobody will act differently on any value',
'why': 'it is a number for a slide, not a forecast'},
]).to_string(index=False))
Open this example in the editor and change it. Nothing you type is saved.
The structural-break case deserves emphasis. Every model in this lesson assumes the future resembles the past. A pricing change, a competitor launch, a regulation or a pandemic breaks that assumption, and no amount of model sophistication helps. Say so in the forecast rather than discovering it afterwards.
A worked example
import numpy as np
import pandas as pd
rng = np.random.default_rng(seed=211)
# Three years of weekly orders: a trend, a yearly season, a December
# spike, and noise.
weeks = pd.date_range('2023-01-01', periods=157, freq='W-SUN')
t = np.arange(len(weeks))
week_of_year = weeks.isocalendar().week.to_numpy()
trend = 2400 + 11 * t
seasonal = 300 * np.sin(2 * np.pi * (t - 8) / 52)
december = np.where((week_of_year >= 48) & (week_of_year <= 52), 900, 0)
noise = rng.normal(0, 110, len(weeks))
orders = pd.Series((trend + seasonal + december + noise).round(), index=weeks, name='orders')
print(f'{len(orders)} weeks, {orders.index.min().date()} to {orders.index.max().date()}')
print(f'mean {orders.mean():,.0f} std {orders.std():,.0f}')
print()
# ---- 1. Decompose, to see what is there. ---------------------------
trend_estimate = orders.rolling(52, center=True, min_periods=26).mean()
detrended = orders - trend_estimate
seasonal_estimate = detrended.groupby(week_of_year).transform('mean')
residual = (detrended - seasonal_estimate).dropna()
print('variation accounted for:')
print(f' trend std {trend_estimate.std():7.1f}')
print(f' seasonal std {seasonal_estimate.std():7.1f}')
print(f' residual std {residual.std():7.1f} <- the floor on any forecast')
print()
# ---- 2. Rolling-origin evaluation of four methods. -----------------
HORIZON = 13
def naive(history, horizon):
return np.repeat(history[-1], horizon)
def drift(history, horizon):
step = (history[-1] - history[0]) / (len(history) - 1)
return history[-1] + step * np.arange(1, horizon + 1)
def seasonal_naive(history, horizon, period=52):
tail = history[-period:]
return np.array([tail[i % period] for i in range(horizon)])
def seasonal_with_drift(history, horizon, period=52):
base = seasonal_naive(history, horizon, period)
recent = history[-period:].mean()
previous = history[-2 * period:-period].mean() if len(history) >= 2 * period else recent
yearly_growth = recent - previous
return base + yearly_growth * np.arange(1, horizon + 1) / period
METHODS = {'naive': naive, 'drift': drift,
'seasonal naive': seasonal_naive, 'seasonal + drift': seasonal_with_drift}
results = {name: [] for name in METHODS}
origins = list(range(110, len(orders) - HORIZON, 4))
for origin in origins:
history = orders.iloc[:origin].to_numpy()
actual = orders.iloc[origin:origin + HORIZON].to_numpy()
for name, method in METHODS.items():
predicted = method(history, HORIZON)
results[name].append(np.mean(np.abs(actual - predicted)))
print(f'rolling-origin evaluation: {len(origins)} origins, {HORIZON}-week horizon')
print(f'{"method":<20} {"MAE":>8} {"vs baseline":>12}')
baseline_mae = np.mean(results['seasonal naive'])
for name, errors in sorted(results.items(), key=lambda kv: np.mean(kv[1])):
mae = np.mean(errors)
print(f'{name:<20} {mae:>8.1f} {mae / baseline_mae:>11.2f}x')
print()
print(f'the residual floor is {residual.std():.0f}, so even the best method')
print('is close to the irreducible noise')
print()
# ---- 3. The chosen method, with an empirical band. -----------------
best = min(results, key=lambda k: np.mean(results[k]))
print(f'chosen: {best}')
per_horizon = {h: [] for h in range(1, HORIZON + 1)}
for origin in origins:
history = orders.iloc[:origin].to_numpy()
predicted = METHODS[best](history, HORIZON)
actual = orders.iloc[origin:origin + HORIZON].to_numpy()
for h in range(1, HORIZON + 1):
per_horizon[h].append(actual[h - 1] - predicted[h - 1])
history = orders.to_numpy()
point = METHODS[best](history, HORIZON)
future = pd.date_range(orders.index[-1] + pd.Timedelta('7D'), periods=HORIZON, freq='W-SUN')
rows = []
for h in range(1, HORIZON + 1):
low, high = np.percentile(per_horizon[h], [10, 90])
rows.append({
'week': future[h - 1].date(), 'point': round(point[h - 1]),
'low': round(point[h - 1] + low), 'high': round(point[h - 1] + high),
})
forecast = pd.DataFrame(rows)
forecast['width'] = forecast['high'] - forecast['low']
print(forecast.to_string(index=False))
print()
print(f'band width: {forecast["width"].iloc[0]:,} at 1 week, '
f'{forecast["width"].iloc[-1]:,} at {HORIZON} weeks')
print()
# ---- 4. Check the band is honest, on history. ----------------------
inside = 0
tested = 0
for origin in origins:
history_slice = orders.iloc[:origin].to_numpy()
predicted = METHODS[best](history_slice, HORIZON)
actual = orders.iloc[origin:origin + HORIZON].to_numpy()
for h in range(1, HORIZON + 1):
low, high = np.percentile(per_horizon[h], [10, 90])
tested += 1
if predicted[h - 1] + low <= actual[h - 1] <= predicted[h - 1] + high:
inside += 1
print(f'coverage of the 80% band on history: {inside / tested:.1%} of {tested} points')
print('(it should be near 80%; well above means the band is too wide,')
print(' well below means it is too narrow and will embarrass you)')
print()
# ---- 5. The sentence. ----------------------------------------------
quarter_total = forecast['point'].sum()
quarter_low, quarter_high = forecast['low'].sum(), forecast['high'].sum()
print(f'Orders for the next {HORIZON} weeks are forecast at '
f'{quarter_total:,.0f} in total (80% interval {quarter_low:,.0f} to '
f'{quarter_high:,.0f}), using {best} evaluated over {len(origins)} '
f'rolling origins. This assumes no structural change: a pricing '
f'change, a competitor launch or a supply interruption would put '
f'the outcome outside this range.')
Open this example in the editor and change it. Nothing you type is saved.
Five things make that a forecast rather than a guess. The decomposition establishes a floor on achievable error. Four methods are compared on data none of them saw, at many origins rather than one. The band comes from the model's own historical errors rather than a distributional assumption. The band's coverage is checked against history. And the final sentence names the assumption that would invalidate it.
When it goes wrong
| Symptom | Cause |
|---|---|
| Excellent fit, terrible forecast | Overfitting; evaluate out of sample |
| The error was computed on training data | Split in time, never randomly |
| A complex model lost to "last year's number" | The seasonal naive is hard to beat |
| The band never widens | The uncertainty is not being modelled |
| The forecast was treated as a commitment | A point estimate; publish a range |
| The model missed a step change | No model extrapolates a break; say so |
| Seasonality was mistaken for trend | Fewer than two full cycles of history |
| Being outside the band caused a crisis | Nobody said an 80% band is outside 20% of the time |
A check you can run
Take your current forecast and compute what "the same period last year" would have predicted.
If the seasonal naive is within a few percent of your model, the model is adding nothing, and the honest move is to ship the one-line baseline and spend the effort elsewhere. In my experience it beats the model more often than anyone expects, and finding that out costs five lines.