Forecast Accuracy · Challenge Solutions¶
Worked solutions to the five challenges from Chapter 135, on the same 24-month forecast comparison.
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 scipy import stats
import warnings; warnings.filterwarnings('ignore')
try: df = pd.read_excel('../../data/forecast-accuracy--demand_forecasts.xlsx', sheet_name='Data')
except FileNotFoundError: df = pd.read_excel(BASE + 'forecast-accuracy--demand_forecasts.xlsx', sheet_name='Data')
df['month']=pd.to_datetime(df['month']); df=df.set_index('month'); a=df['actual'].values
resid={c:a-df[c].values for c in ['model_a','model_b','naive']}
print('loaded', len(df), 'months')
loaded 24 months
def mae(r): return np.abs(r).mean()
def rmse(r): return np.sqrt((r**2).mean())
print(pd.DataFrame({c:{'ME':resid[c].mean(),'MAE':mae(resid[c]),'RMSE':rmse(resid[c])} for c in resid}).T.round(0))
print('\nmodel_b and naive have large positive ME -> they systematically UNDER-forecast (biased).')
ME MAE RMSE model_a 84.0 218.0 564.0 model_b 496.0 496.0 532.0 naive 827.0 827.0 1013.0 model_b and naive have large positive ME -> they systematically UNDER-forecast (biased).
ME isolates bias: model_b (+496) and naive (+827) lean low, while model_a's ME is near zero. But model_a is still the most accurate, because ME says nothing about the size of the swings, that is what MAE and RMSE are for.
for c in resid: print('%-8s RMSE - MAE = %6.1f' % (c, rmse(resid[c])-mae(resid[c])))
big = np.argmax(np.abs(resid['model_a'])); print('\nmodel_a largest single error:', round(resid['model_a'][big]), 'at', df.index[big].date(), '(the missed spike)')
print('RMSE squares errors, so that one miss dominates -> RMSE (564) sits far above MAE (218).')
model_a RMSE - MAE = 346.5 model_b RMSE - MAE = 35.8 naive RMSE - MAE = 185.3 model_a largest single error: 2685 at 2023-03-01 (the missed spike) RMSE squares errors, so that one miss dominates -> RMSE (564) sits far above MAE (218).
Because RMSE squares before averaging, RMSE is always at least MAE, and the gap widens with the spread of the errors. model_a's gap is large because of a single missed spike, one big error can dominate RMSE while barely moving MAE.
def mape(f): return np.mean(np.abs((a-f)/a))*100
def smape(f): return np.mean(2*np.abs(a-f)/(np.abs(a)+np.abs(f)))*100
for c in ['model_a','model_b','naive']: print('%-8s MAPE %5.2f%% sMAPE %5.2f%%' % (c, mape(df[c].values), smape(df[c].values)))
print('\nMAPE pitfalls: undefined/explosive when actual ~ 0, and it penalizes OVER-forecasts more than UNDER.')
print('sMAPE bounds the over-forecast penalty but is still asymmetric; MASE avoids the percentage trap entirely.')
model_a MAPE 2.85% sMAPE 3.04% model_b MAPE 7.16% sMAPE 7.46% naive MAPE 11.73% sMAPE 12.71% MAPE pitfalls: undefined/explosive when actual ~ 0, and it penalizes OVER-forecasts more than UNDER. sMAPE bounds the over-forecast penalty but is still asymmetric; MASE avoids the percentage trap entirely.
MAPE is intuitive and scale-free but breaks down when actuals approach zero and treats over- and under-forecasts asymmetrically. sMAPE averages actual and forecast in the denominator, which bounds the over-forecast penalty, though it is still asymmetric; for series with zeros or tiny values, prefer MASE instead.
nm, nr = mae(resid['naive']), rmse(resid['naive'])
for c in ['model_a','model_b','naive']:
print('%-8s MASE %.2f Theil_U %.2f -> %s' % (c, mae(resid[c])/nm, rmse(resid[c])/nr, 'beats naive' if mae(resid[c])/nm<1 else 'benchmark'))
print('\nboth models are well below 1: they add real value over doing nothing.')
model_a MASE 0.26 Theil_U 0.56 -> beats naive model_b MASE 0.60 Theil_U 0.53 -> beats naive naive MASE 1.00 Theil_U 1.00 -> benchmark both models are well below 1: they add real value over doing nothing.
MASE and Theil's U rescale error by the naive benchmark (here the seasonal-naive forecast's MAE and RMSE on the same holdout, a relative MAE and RMSE; textbook MASE scales by the in-sample naive MAE), so below 1 means you beat naive. Both models pass easily (MASE 0.26 and 0.60). This relative check is the first hurdle: a model that cannot beat seasonal-naive is not worth deploying.
def dm(a,f1,f2):
d=(a-f1)**2-(a-f2)**2; n=len(d); dbar=d.mean(); var=np.mean((d-dbar)**2)
stat=dbar/np.sqrt(var/n); return stat, 2*(1-stats.norm.cdf(abs(stat)))
for pair in [('model_a','model_b'),('model_a','naive')]:
s,p=dm(a, df[pair[0]].values, df[pair[1]].values)
print('%-8s vs %-8s DM stat %6.2f p %.4f -> %s' % (pair[0], pair[1], s, p, 'A significantly better' if p<0.05 and s<0 else 'not significant'))
model_a vs model_b DM stat 0.13 p 0.8967 -> not significant model_a vs naive DM stat -3.85 p 0.0001 -> A significantly better
The Diebold-Mariano test asks whether two forecasts' squared errors differ beyond chance. A vs B is not significant (p about 0.90, the gap rides on one spike), but A vs naive is (A is reliably better). A ranking without a significance check can crown a winner that is really a coin flip.