Forecasting Retail Sales: Take It Further¶
Five extensions of the Chapter 136 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)
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
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.
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.
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.
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.
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 = +/-6.8% -> month 12 = +/-9.9% 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 about plus or minus 10% 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.