←Back to the bookChapter 131 · Take it further

Chapter 131 · Take It Further

Forecasting Retail Sales · Take It Further

Five ways to stress-test the forecast: a rolling-origin backtest, a prediction-interval coverage check, a log-transform for SARIMA, an automatic AIC order search, and how forecast error grows with the horizon.

Statistics, Data Science and AI: A Visual Handbook · John Fisher · 2026

Forecasting Retail Sales: Take It Further

Five extensions of the Chapter 131 project, each stress-testing the forecast: a rolling-origin backtest instead of one split, a check of whether the 95% interval really covers 95%, a log-transform to help SARIMA, an automatic order search, and how far the forecast can be trusted as the horizon grows. We start by rebuilding the cleaned series.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import seaborn as sns   # seaborn = high-level statistical plots (heatmaps, pairplots, count/bar plots)
from matplotlib.colors import ListedColormap
EM="#0284c7"; DEEP="#075985"; LIGHT="#bae6fd"; INK="#1a2138"; GRID="#e6e9f2"; RED="#ef4444"; AMBER="#d97706"; GREEN="#059669"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"; ORG="#0284c7"; CYAN="#0891b2"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"

from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings; warnings.filterwarnings('ignore')
def mape(a,f): return float(np.mean(np.abs((np.asarray(a)-np.asarray(f))/np.asarray(a)))*100)

RECAP · REBUILD

The cleaned series

Load, fix the bad month, split.

try: raw = pd.read_excel('../../data/case-study-forecasting-retail-sales--retail_sales.xlsx', sheet_name='Data')
except FileNotFoundError: raw = pd.read_excel(BASE + 'case-study-forecasting-retail-sales--retail_sales.xlsx', sheet_name='Data')
raw['month']=pd.to_datetime(raw['month']); s = raw.set_index('month')['sales'].asfreq('MS')
med=s.rolling(13,center=True,min_periods=3).median(); s[(s-med).abs().idxmax()]=np.nan; clean=s.interpolate()
def hw(tr): return ExponentialSmoothing(tr, trend='add', seasonal='mul', seasonal_periods=12).fit()
print('rebuilt cleaned series:', len(clean), 'months')

rebuilt cleaned series: 108 months

EXTENSION 1 · ROLLING-ORIGIN BACKTEST

One split can lie

Re-forecast at many cutoffs and average.

origins = range(len(clean)-24, len(clean)-6, 3)   # several forecast origins near the end
scores=[]
for o in origins:
    tr, te = clean.iloc[:o], clean.iloc[o:o+6]
    scores.append(mape(te, hw(tr).forecast(6)))
fig,ax=plt.subplots(figsize=(8,3.8)); ax.bar(range(len(scores)), scores, color=EM)
ax.axhline(np.mean(scores), color=RED, ls='--', label=f'mean {np.mean(scores):.2f}%'); ax.set(title='Holt-Winters MAPE across rolling origins (6-month horizon)', xlabel='backtest window', ylabel='MAPE %'); ax.legend()
plt.tight_layout(); plt.show()
print('average MAPE over %d backtests: %.2f%% (more trustworthy than a single split)' % (len(scores), np.mean(scores)))

average MAPE over 6 backtests: 2.36% (more trustworthy than a single split)

A single train/test split can flatter or punish a model by luck. Rolling-origin evaluation re-forecasts from several cutoffs and averages the error, the standard way to grade a forecaster. The average across windows is the number to report, not one lucky split.

EXTENSION 2 · DOES THE 95% INTERVAL HOLD?

An interval you can bank on

Backtest the coverage, not just the point.

hits=0; total=0
for o in range(len(clean)-30, len(clean)-6, 2):
    tr = clean.iloc[:o]; f = hw(tr); sims = f.simulate(6, repetitions=500, error='mul')
    lo, hi = sims.quantile(0.025,axis=1).values, sims.quantile(0.975,axis=1).values
    act = clean.iloc[o:o+6].values
    hits += int(((act>=lo)&(act<=hi)).sum()); total += len(act)
print('actuals inside the 95%% interval: %d of %d = %.0f%%' % (hits, total, 100*hits/total))
print('close to 95%% -> the interval is well calibrated and safe to plan against')

actuals inside the 95% interval: 70 of 72 = 97%
close to 95%% -> the interval is well calibrated and safe to plan against

A prediction interval is a promise: 95% of actuals should fall inside it. Backtesting the coverage checks whether that promise holds, if only 70% land inside, the band is too narrow and the plan is over-confident. Here coverage lands near 95%, so the interval is trustworthy.

EXTENSION 3 · A LOG-TRANSFORM FOR SARIMA

Turn multiplicative into additive

SARIMA is an additive model; give it additive data.

train, test = clean[:-12], clean[-12:]
sar = SARIMAX(train, order=(1,1,1), seasonal_order=(1,1,1,12), enforce_stationarity=False, enforce_invertibility=False).fit(disp=False)
sar_log = SARIMAX(np.log(train), order=(1,1,1), seasonal_order=(1,1,1,12), enforce_stationarity=False, enforce_invertibility=False).fit(disp=False)
print('SARIMA on raw   MAPE: %.2f%%' % mape(test, sar.forecast(12)))
print('SARIMA on log   MAPE: %.2f%%' % mape(test, np.exp(sar_log.forecast(12))))
print('logging stabilizes the growing seasonal swing, so additive SARIMA fits better')

SARIMA on raw   MAPE: 4.88%
SARIMA on log   MAPE: 3.61%
logging stabilizes the growing seasonal swing, so additive SARIMA fits better

Our seasonality is multiplicative (the swing grows with the level), but SARIMA is fundamentally an additive model. Taking logs stabilizes the swing, log(T x S) = log T + log S, so SARIMA on the logged series is a fairer fight and usually improves. Back-transform with exp to read the forecast in dollars.

EXTENSION 4 · AUTOMATIC ORDER SEARCH

Let AIC pick the orders

A small grid beats guessing (p,d,q).

import itertools
best=(None, np.inf)
for p,q,P,Q in itertools.product([0,1],[0,1],[0,1],[0,1]):
    try:
        m=SARIMAX(np.log(train), order=(p,1,q), seasonal_order=(P,1,Q,12), enforce_stationarity=False, enforce_invertibility=False).fit(disp=False)
        if m.aic<best[1]: best=((p,1,q,P,1,Q), m.aic)
    except Exception: pass
print('lowest-AIC order (p,d,q,P,D,Q):', best[0], ' AIC %.0f' % best[1])
print('this is the idea behind auto_arima: search a small grid, pick the best by AIC, no manual guessing')

lowest-AIC order (p,d,q,P,D,Q): (0, 1, 1, 0, 1, 1)  AIC -266
this is the idea behind auto_arima: search a small grid, pick the best by AIC, no manual guessing

Rather than guess the orders, search a small grid and pick the combination with the lowest AIC (which balances fit against complexity, without touching the test set). This is exactly what tools like pmdarima.auto_arima automate; here it confirms a compact, sensible order.

EXTENSION 5 · HOW FAR CAN YOU TRUST IT?

Uncertainty widens with the horizon

See the actual, the forecast, and the fanning 95% band together.

train, test = clean[:-12], clean[-12:]
f = hw(train); fc = f.forecast(12)
sims = f.simulate(12, repetitions=3000, error='mul')
lo, hi = sims.quantile(0.025,axis=1), sims.quantile(0.975,axis=1)
fig,ax=plt.subplots(figsize=(10,4.4))
ax.plot(clean.index[-30:], clean.values[-30:], color=GREY, lw=2.2, label='actual')
ax.plot(test.index, fc, 'o--', color=EM, lw=2, ms=4, label='forecast')
ax.fill_between(test.index, lo, hi, color=EM, alpha=0.20, label='95% interval')
ax.axvline(test.index[0], color=RED, ls=':', lw=1.2)
ax.text(test.index[0], ax.get_ylim()[0], ' forecast starts', color=RED, fontsize=8, va='bottom')
ax.set(title='Actual vs forecast on the holdout: the 95% band fans out the further ahead you go', ylabel='sales ($)')
ax.legend(loc='upper left'); plt.tight_layout(); plt.show()
hw_pct = ((hi - lo) / 2 / fc * 100).values
inside = ((test.values >= lo.values) & (test.values <= hi.values)).mean()*100
print('interval half-width: month 1 = +/-%.1f%%  ->  month 12 = +/-%.1f%%' % (hw_pct[0], hw_pct[-1]))
print('actuals falling inside the band: %.0f%%  -> accurate near-term, honestly less CERTAIN far out' % inside)

interval half-width: month 1 = +/-7.1%  ->  month 12 = +/-9.7%
actuals falling inside the band: 100%  -> accurate near-term, honestly less CERTAIN far out

This is the picture to keep. The forecast (blue) tracks the actual holdout year (gray) closely, and the 95% band fans out from about plus or minus 7% next month to plus or minus 9% a year ahead. The actuals stay inside the band the whole way, so the interval is honest, and its widening is the true signal of horizon risk: even when the point forecast stays accurate, the uncertainty compounds. Trust the near term more, and always plan against the widening band, not the single line.

Take-it-further summary, in plain terms

• Rolling-origin backtesting grades a forecaster over many cutoffs, not one lucky split.
• Check interval coverage: 95% intervals should contain about 95% of actuals, or the plan is over-confident.
• Log-transform turns multiplicative seasonality into additive, helping SARIMA.
• Search orders by AIC instead of guessing (the auto_arima idea).
• Error grows with the horizon: trust near-term forecasts more and widen the band far out.

---

Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher

←Back to the bookChapter 131 · Take it further