Forecast Accuracy¶
You have built many forecasts; now grade them fairly. This chapter defines the error metrics, ME, MAE, MSE, RMSE, MAPE, sMAPE, MASE, and Theil's U, shows how they disagree on purpose, and compares two models with the Diebold-Mariano test. The lesson underneath: a metric is a loss function, so the “best” model depends on which errors cost you the most.
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
from statsmodels.stats.diagnostic import acorr_ljungbox
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
fig,ax=plt.subplots(figsize=(10,4)); ax.plot(df.index, a, color='k', lw=2.4, label='actual')
for c,col in zip(['model_a','model_b','naive'],[EM,PUR,GREY]): ax.plot(df.index, df[c], '--', color=col, label=c)
ax.set(title='Actual vs three forecasts', ylabel='demand'); ax.legend(); plt.tight_layout(); plt.show()
print('model_a tracks closely (but watch month 2023-03... it missed a spike); model_b runs low; naive lags a year behind')
model_a tracks closely (but watch month 2023-03... it missed a spike); model_b runs low; naive lags a year behind
Three forecasts of the same 24 months: model_a is accurate but missed one demand spike (a promotion), model_b runs systematically low, and naive just repeats last year. Which is best? It depends entirely on how you score them.
resid = {c: a - df[c].values for c in ['model_a','model_b','naive']}
for c,r in resid.items(): print('%-8s ME (mean error) = %+7.1f' % (c, r.mean()))
fig,ax=plt.subplots(figsize=(10,3.6))
for c,col in zip(resid,[EM,PUR,GREY]): ax.plot(df.index, resid[c], 'o-', color=col, ms=3, label=c)
ax.axhline(0,color='k'); ax.set(title='Residuals (actual - forecast): above 0 = under-forecast', ylabel='error'); ax.legend(); plt.tight_layout(); plt.show()
model_a ME (mean error) = +83.8 model_b ME (mean error) = +496.0 naive ME (mean error) = +827.2
The residual is actual - forecast. The Mean Error (ME) averages them, so it measures bias, not accuracy. model_a's ME is near zero (errors cancel), but model_b (+496) and naive (+827) sit well above zero: they systematically under-forecast. A near-zero ME can still hide big errors that cancel, so ME alone is never enough, it only catches a lean.
def mae(r): return np.abs(r).mean()
def rmse(r): return np.sqrt((r**2).mean())
tab = pd.DataFrame({c:{'MAE':mae(resid[c]),'MSE':(resid[c]**2).mean(),'RMSE':rmse(resid[c])} for c in resid}).T.round(0)
print(tab)
print('\nBy MAE: model_a wins (%.0f). By RMSE: model_b EDGES model_a (%.0f vs %.0f)!' % (mae(resid['model_a']), rmse(resid['model_b']), rmse(resid['model_a'])))
print('why: model_a\'s single missed spike is squared by RMSE, so one big miss dominates.')
fig,ax=plt.subplots(figsize=(8,4)); x=np.arange(3); w=0.38
ax.bar(x-w/2,[mae(resid[c]) for c in resid],w,label='MAE',color=EM); ax.bar(x+w/2,[rmse(resid[c]) for c in resid],w,label='RMSE',color=RED)
ax.set_xticks(x); ax.set_xticklabels(list(resid)); ax.set(title='RMSE punishes big errors more than MAE', ylabel='error'); ax.legend(); plt.tight_layout(); plt.show()
MAE MSE RMSE model_a 218.0 318210.0 564.0 model_b 496.0 282799.0 532.0 naive 827.0 1025332.0 1013.0 By MAE: model_a wins (218). By RMSE: model_b EDGES model_a (532 vs 564)! why: model_a's single missed spike is squared by RMSE, so one big miss dominates.
These metrics use the error's own units. MAE (mean absolute error) treats every miss equally. MSE and its square root RMSE square the errors first, so a single large miss weighs far more. That is why the ranking flips: model_a wins on MAE (218 vs 496) but its one missed spike inflates its RMSE above model_b's (564 vs 532). RMSE is always at least MAE, and the gap grows with the spread of the errors. Choosing between them is really choosing how much a big miss should hurt.
naive_mae = mae(resid['naive']); naive_rmse = rmse(resid['naive'])
def mape(f): r=a-f; return np.mean(np.abs(r/a))*100
def smape(f): r=a-f; return np.mean(2*np.abs(r)/(np.abs(a)+np.abs(f)))*100
rows=[]
for c in ['model_a','model_b','naive']:
f=df[c].values; r=resid[c]
rows.append({'model':c,'MAPE_%':round(mape(f),2),'sMAPE_%':round(smape(f),2),'MASE':round(mae(r)/naive_mae,2),'Theil_U':round(rmse(r)/naive_rmse,2)})
print(pd.DataFrame(rows).to_string(index=False))
print('\nMASE-style relative MAE / Theil U below 1 => beats the seasonal-naive benchmark. Both models do; naive is 1.00 by construction.')
model MAPE_% sMAPE_% MASE Theil_U model_a 2.85 3.04 0.26 0.56 model_b 7.16 7.46 0.60 0.53 naive 11.73 12.71 1.00 1.00 MASE-style relative MAE / Theil U below 1 => beats the seasonal-naive benchmark. Both models do; naive is 1.00 by construction.
To compare across products or scales you need unit-free scores. MAPE averages the percentage error, intuitive, but it blows up near zero and punishes over- and under-forecasts asymmetrically; sMAPE bounds the over-forecast penalty but is still asymmetric and explodes when actual and forecast are both near zero. MASE and Theil's U are relative to the naive benchmark. Strictly, MASE scales by the naive forecast's MAE on the training data. With only the holdout here we scale by the seasonal-naive forecast's MAE on the same 24 months, a relative MAE, so naive scores exactly 1. below 1 means you beat naive, above 1 means you did worse than doing nothing. Both models clear the bar (MASE 0.26 and 0.60), which is the first thing any forecast must prove.
def dm_test(a, f1, f2, h=1):
d=(a-f1)**2 - (a-f2)**2; n=len(d); dbar=d.mean(); var=np.mean((d-dbar)**2)
for k in range(1,h): var += 2*(1-k/h)*np.mean((d[k:]-dbar)*(d[:-k]-dbar))
stat=dbar/np.sqrt(var/n); return stat, 2*(1-stats.norm.cdf(abs(stat)))
s,p = dm_test(a, df['model_a'].values, df['model_b'].values)
print('Diebold-Mariano A vs B (squared-error loss): stat %.2f, p = %.3f' % (s,p))
print('p > 0.05 -> the RMSE gap between A and B is NOT statistically significant (driven by one spike).')
lb = acorr_ljungbox(resid['model_a'], lags=[6])['lb_pvalue'].iloc[0]
print('\nmodel_a residual Ljung-Box p = %.2f -> %s' % (lb, 'no evidence of leftover autocorrelation' if lb>0.05 else 'structure remains'))
Diebold-Mariano A vs B (squared-error loss): stat 0.13, p = 0.897 p > 0.05 -> the RMSE gap between A and B is NOT statistically significant (driven by one spike). model_a residual Ljung-Box p = 0.99 -> no evidence of leftover autocorrelation
A ranking is only as good as its reliability. The Diebold-Mariano test asks whether two forecasts' errors differ significantly. Here A vs B gives p about 0.90: the RMSE gap is not statistically distinguishable, it hangs on that one spike, so you should not crown B on RMSE alone. Good practice also checks the residuals: model_a's show no evidence of leftover autocorrelation in a Ljung-Box test (p about 0.99 on 24 months). That is not proof of white noise, the test covers only autocorrelation, and the missed promotion spike is itself signal the model did not use.
panel=[]
for c in ['model_a','model_b','naive']:
f=df[c].values; r=resid[c]
panel.append({'model':c,'ME':round(r.mean()),'MAE':round(mae(r)),'RMSE':round(rmse(r)),'MAPE%':round(mape(f),1),'MASE':round(mae(r)/naive_mae,2)})
print(pd.DataFrame(panel).to_string(index=False))
print('\nTypical-accuracy view (MAE, MAPE): model_a wins clearly.')
print('Large-miss-averse view (RMSE): A and B cannot be told apart on this sample, because one spike dominates.')
print('Both beat naive. If missing rare spikes is costly, invest in catching them; otherwise ship model_a.')
model ME MAE RMSE MAPE% MASE model_a 84 218 564 2.8 0.26 model_b 496 496 532 7.2 0.60 naive 827 827 1013 11.7 1.00 Typical-accuracy view (MAE, MAPE): model_a wins clearly. Large-miss-averse view (RMSE): A and B cannot be told apart on this sample, because one spike dominates. Both beat naive. If missing rare spikes is costly, invest in catching them; otherwise ship model_a.
There is no single “accuracy” number. On typical accuracy (MAE, MAPE) model_a is the clear winner; if a single large miss is what hurts (RMSE), A and B cannot be told apart on these 24 months. The right call follows the cost of your errors: report several metrics, check they agree, test the difference, and pick the metric that mirrors the decision.
The accuracy toolkit, in one view¶
- ME measures bias (a systematic lean), not accuracy; near-zero can still hide canceling errors.
- MAE / MSE / RMSE are scale-dependent; RMSE punishes big misses hardest (RMSE always at least MAE).
- MAPE / sMAPE are percentage errors, comparable across scales but fragile near zero.
- MASE / Theil's U score you against the naive benchmark (here as relative MAE and RMSE on the holdout), below 1 is the bar every forecast must clear.
- Diebold-Mariano tests whether one model is significantly better, a ranking gap can be noise.
The one idea to keep: a forecast metric is a loss function in disguise, so choosing a metric is choosing which errors you refuse to tolerate, and the honest report shows several, confirms they agree, and tests the difference.