Case Study: From Forecast to Decision¶
A forecast is not the answer, it is an input to a decision: how many units to order tonight. And the best order is almost never the forecast itself, because being short and being over cost different amounts. This case study compares two forecasts honestly (accuracy panel plus a Diebold-Mariano test), then turns the winner into an order quantity with the newsvendor model, where the optimal order is a quantile of demand set by the cost ratio, not the average.
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.tsa.holtwinters import ExponentialSmoothing
import warnings; warnings.filterwarnings('ignore')
def mape(a,f): return float(np.mean(np.abs((a.values-np.asarray(f))/a.values))*100)
def mae(a,f): return float(np.mean(np.abs(a.values-np.asarray(f))))
def rmse(a,f): return float(np.sqrt(np.mean((a.values-np.asarray(f))**2)))
A store orders a perishable product once a day. If it runs short, it loses the margin on each missed sale, the underage cost, here about 8 dollars per unit. If it over-orders, leftover units are wasted, the overage cost, about 2 dollars per unit. These are asymmetric: a stockout hurts four times as much as a leftover. Success is minimum total cost, not minimum forecast error.
try: raw = pd.read_excel('../../data/case-study-from-forecast-to-decision--product_demand.xlsx', sheet_name='Data')
except FileNotFoundError: raw = pd.read_excel(BASE + 'case-study-from-forecast-to-decision--product_demand.xlsx', sheet_name='Data')
raw['date']=pd.to_datetime(raw['date']); s = raw.set_index('date')['demand'].asfreq('D')
print('days:', len(s), '| mean %.0f, std %.0f units/day' % (s.mean(), s.std()))
days: 731 | mean 220, std 39 units/day
H=60; train, test = s[:-H], s[-H:]
fig,ax=plt.subplots(figsize=(11,3.8)); ax.plot(s.index, s.values, color=EM, lw=0.8)
ax.axvspan(test.index[0], test.index[-1], color=AMBER, alpha=0.10); ax.set(title='Daily demand: weekend peaks, and plenty of day-to-day variability', ylabel='units')
plt.tight_layout(); plt.show()
print('weekend avg %d vs weekday avg %d units' % (s[s.index.dayofweek>=5].mean(), s[s.index.dayofweek<5].mean()))
weekend avg 250 vs weekday avg 207 units
Demand peaks on weekends and carries a slow trend, but the star of this chapter is the scatter: even a perfect forecast of the average leaves real day-to-day variability. That uncertainty is exactly what the order decision must hedge against.
snaive = np.array([train.iloc[-7 + (i % 7)] for i in range(H)]) # repeat last week
hw = ExponentialSmoothing(train, trend='add', seasonal='add', seasonal_periods=7).fit()
hw_fc = hw.forecast(H).values
print('two forecasts produced for the 60-day holdout: seasonal-naive and Holt-Winters')
two forecasts produced for the 60-day holdout: seasonal-naive and Holt-Winters
rows=[]
for name,f in [('seasonal-naive',snaive),('Holt-Winters',hw_fc)]:
rows.append({'model':name,'MAE':round(mae(test,f),1),'RMSE':round(rmse(test,f),1),'MAPE_%':round(mape(test,f),2)})
print(pd.DataFrame(rows).to_string(index=False))
model MAE RMSE MAPE_% seasonal-naive 44.1 58.8 18.64 Holt-Winters 20.0 25.1 9.04
On the held-out 60 days, Holt-Winters is clearly more accurate, about 9% MAPE versus 19% for seasonal-naive. But a ranking is only trustworthy if the gap is real, so we test it.
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)))
st, p = dm(test.values, snaive, hw_fc)
print('Diebold-Mariano (seasonal-naive vs Holt-Winters): stat %.2f, p = %.4f' % (st, p))
print('p < 0.05 and stat > 0 -> Holt-Winters is SIGNIFICANTLY more accurate. We use it.')
Diebold-Mariano (seasonal-naive vs Holt-Winters): stat 4.34, p = 0.0000 p < 0.05 and stat > 0 -> Holt-Winters is SIGNIFICANTLY more accurate. We use it.
The Diebold-Mariano test confirms the difference is statistically significant (p well below 0.05), not a fluke of one holdout. Holt-Winters wins, so it feeds the decision. Now the real work begins.
resid_sd = float(np.std(train.values - hw.fittedvalues.values))
print('Holt-Winters point forecast is a best guess; actual demand scatters around it.')
print('one-step forecast error std ~ %.1f units -> that spread is what the order must cover' % resid_sd)
Holt-Winters point forecast is a best guess; actual demand scatters around it. one-step forecast error std ~ 26.6 units -> that spread is what the order must cover
The point forecast is the middle of a distribution of possible demand, not a certainty. Its spread (an error standard deviation of about 27 units) is the uncertainty the order has to hedge. If both errors cost the same, you would order the middle; because they do not, you should not.
Cu, Co = 8.0, 2.0 # underage (stockout) and overage (waste), $/unit
fractile = Cu / (Cu + Co)
z = stats.norm.ppf(fractile)
print('critical fractile = Cu / (Cu + Co) = %.0f / %.0f = %.2f' % (Cu, Cu+Co, fractile))
print('order to the %.0fth percentile of demand (z = %.2f), NOT the 50th (the mean)' % (fractile*100, z))
print('because a stockout costs %.0fx a leftover, carry a buffer ABOVE the forecast' % (Cu/Co))
critical fractile = Cu / (Cu + Co) = 8 / 10 = 0.80 order to the 80th percentile of demand (z = 0.84), NOT the 50th (the mean) because a stockout costs 4x a leftover, carry a buffer ABOVE the forecast
The newsvendor model gives the answer in one line: the cost-minimizing order is the critical fractile quantile of demand, where the fractile is Cu / (Cu + Co). Here that is 0.80, order to the 80th percentile, well above the average. The intuition is simple: when running short is four times as expensive as overstocking, you deliberately carry a safety buffer.
q_point = hw_fc # order the point forecast (50th percentile)
q_star = hw_fc + z * resid_sd # order the 80th percentile
def cost(demand, order): return Cu*np.maximum(demand-order,0) + Co*np.maximum(order-demand,0)
fr = np.linspace(0.30, 0.97, 60); avg_cost=[cost(test.values, hw_fc + stats.norm.ppf(q)*resid_sd).mean() for q in fr]
best = fr[int(np.argmin(avg_cost))]
fig,ax=plt.subplots(figsize=(9,4)); ax.plot(fr*100, avg_cost, color=EM, lw=2.4)
ax.axvline(50, color=GREY, ls=':', label='order the forecast (50%)'); ax.axvline(fractile*100, color=RED, ls='--', label='cost-optimal (%.0f%%)'%(fractile*100))
ax.set(xlabel='order quantile (percentile of demand)', ylabel='avg daily cost ($)', title='Expected cost is minimized ABOVE the point forecast'); ax.legend(); plt.tight_layout(); plt.show()
print('empirical cost-minimizing fractile on the holdout: %.0f%% (theory said %.0f%%)' % (best*100, fractile*100))
empirical cost-minimizing fractile on the holdout: 78% (theory said 80%)
Sweeping the order across quantiles traces a clear cost curve, and its minimum sits near the 80th percentile, exactly where the newsvendor formula predicted, not at the 50th where the point forecast lives. Ordering the average is a common, expensive mistake; the cost curve shows why.
c_point = cost(test.values, q_point).mean(); c_star = cost(test.values, q_star).mean()
print('order = point forecast : $%.0f/day | service level %.0f%% (days demand fully met)' % (c_point, 100*np.mean(test.values<=q_point)))
print('order = 80%% fractile : $%.0f/day | service level %.0f%%' % (c_star, 100*np.mean(test.values<=q_star)))
print('the quantile order is %.0f%% cheaper AND meets demand far more often' % (100*(c_point-c_star)/c_point))
out_pt = test.values > q_point; out_both = out_pt & (test.values > q_star)
print('stockouts if ordering the forecast: %d of %d days; the buffer prevents %d of them' % (out_pt.sum(), len(test), out_pt.sum()-out_both.sum()))
fig,ax=plt.subplots(figsize=(11,3.8)); ax.plot(test.index, test.values, color='k', lw=1.6, label='actual demand')
ax.plot(test.index, q_point, '--', color=GREY, label='order = forecast'); ax.plot(test.index, q_star, '--', color=EM, lw=2, label='order = 80% fractile')
sh = test.values > q_point; ax.plot(test.index[sh], test.values[sh], 'v', color=RED, ms=6, label='stockout if ordering the forecast')
ax.set(title='The fractile order sits above the forecast and prevents most stockouts', ylabel='units'); ax.legend(fontsize=8); plt.tight_layout(); plt.show()
order = point forecast : $100/day | service level 47% (days demand fully met) order = 80% fractile : $69/day | service level 83% the quantile order is 31% cheaper AND meets demand far more often stockouts if ordering the forecast: 32 of 60 days; the buffer prevents 22 of them
Backtested over the holdout, ordering the point forecast leaves the shelf empty on more than half the days (service level about 47%) and costs about 100 dollars a day. The 80th-percentile order lifts service to about 83% and cuts cost to about 69 dollars a day, roughly 31% cheaper. The red markers are the stockouts from ordering the forecast, and the buffer prevents most of them (the counts are printed above). Better decisions, from the same forecast.
def order_quantity(point_forecast, error_sd, Cu, Co):
z = stats.norm.ppf(Cu/(Cu+Co)); return point_forecast + z*error_sd
final = ExponentialSmoothing(s, trend='add', seasonal='add', seasonal_periods=7).fit()
tomorrow = float(final.forecast(1).iloc[0]); sd = float(np.std(s.values - final.fittedvalues.values))
q = order_quantity(tomorrow, sd, Cu, Co)
print('tomorrow: point forecast %.0f units | recommended ORDER %.0f units (a %.0f-unit safety buffer)' % (tomorrow, q, q-tomorrow))
print('expected: fully stocked ~80%% of such days, at the lowest long-run cost given the 4:1 cost ratio')
tomorrow: point forecast 198 units | recommended ORDER 220 units (a 22-unit safety buffer)
expected: fully stocked ~80%% of such days, at the lowest long-run cost given the 4:1 cost ratio
Deployment is a single rule: refit the forecast on all data, read tomorrow's point forecast and the error spread, and return the critical-fractile order. It reruns every evening. Change the cost ratio and the buffer changes automatically, which is exactly what you want when margins or waste costs shift.
For the store manager: our demand forecast is good (accurate to about 9%), but you should not order the forecast number. Because a stockout costs about four times a leftover, the cheapest policy is to order to the 80th percentile, roughly 22 units above the forecast. Doing so keeps the shelf stocked about 83% of days instead of under half, and cuts total cost by about 31%. Order the buffer, not the average, and we will revisit it if margins or waste costs change.
From forecast to decision, in one view¶
- Define the decision and its asymmetric costs (underage vs overage). 2-3. Collect and inspect the demand, noting its day-to-day scatter. 4-6. Forecast with competing models, compare accuracy, and test the gap with Diebold-Mariano.
- Quantify uncertainty: the point forecast is the middle of a distribution.
- Decision model: the newsvendor optimal order is the Cu/(Cu+Co) quantile of demand. 9-10. Compare and backtest: the quantile order beats ordering the average on cost and service. 11-12. Deploy a one-line order rule and communicate the buffer in plain English.
The one idea to keep: the most accurate forecast is not the decision, the best order is a quantile of demand set by the cost of being wrong in each direction, so accuracy earns you a good forecast and the newsvendor formula turns it into a good decision.