←Back to the bookChapter 199 · Hierarchical Forecasting: Store, Region, National

Chapter 199 · Capstone 36 · Hierarchical Forecasting

Store, Region, National

Three teams forecast the same year and produce three different totals. None of them is wrong. All of them have to be replaced before anyone can be given a target.

import numpy as np, pandas as pd, warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
from statsmodels.tsa.holtwinters import ExponentialSmoothing
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False,
    "axes.grid":True,"grid.alpha":0.22,"axes.titleweight":"bold","axes.titlesize":12.5,"axes.titlelocation":"left"})
SK, DK, LT, MUT, GD, RD = "#0369a1", "#0c4a6e", "#7dd3fc", "#94a3b8", "#047857", "#dc2626"

BASE_URL = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
fn = "capstone-hierarchical-forecasting.xlsx"
def load(sheet):
    try:    return pd.read_excel("../../data/" + fn, sheet_name=sheet)
    except FileNotFoundError: return pd.read_excel(BASE_URL + fn, sheet_name=sheet)
raw   = load("Data")
notes = load("Notes")
REG = ["Northeast", "Midwest", "South", "West"]
print("rows in the planning extract:", f"{len(raw):,}")
raw.head()

rows in the planning extract: 1,360

	month	store_id	region	units
0	2019-01-01	MI1	Midwest	6124
1	2019-01-01	MI2	Midwest	5180
2	2019-01-01	MI3	Midwest	5811
3	2019-01-01	MI4	Midwest	3346
4	2019-01-01	NO1	Northeast	3706

Step 1 · The brief

for line in notes.Notes.fillna(""):
    print(line)

CHAIN PLANNING EXTRACT. Monthly units by store, January 2019 to December 2025.

THE HIERARCHY. Sixteen stores in four regions, and one chain. Region totals are the sum of
their stores and the chain total is the sum of the regions, exactly. Nothing is allocated or
estimated: units are recorded at the till and every level above the store is arithmetic.

THE JOB. Produce a plan for the next twelve months at all three levels. The regional managers
are held to the regional number, the store managers to the store number, and the board is shown
the chain number. Those three numbers have to be consistent with each other, and when each is
forecast on its own they are not.

WHAT IS IN THE EXPORT, as it arrived:
 - the August 2022 extract ran twice, so those sixteen rows appear twice
 - store MI3 has its region written three ways across the years, one of them lower case and one
   with spaces around it
 - store NO4 shows -412 units in March 2021. That is a stocktake adjustment, not a sale
 - stores SO2 (June 2020) and WE3 (September 2023) show zero units. The till export failed;
   neither store was closed

The zeros matter more than they look. Left alone they drag a store's trend estimate down and
the error propagates into every level above it.

Two things in there decide the shape of the work.

The hierarchy is exact. Units are recorded once, at the till. Everything above a store is arithmetic, so there is no measurement question about whether the levels agree. They agree by construction, in the data. It is the forecasts that will not.

Three audiences are held to three numbers. That is what makes coherence a requirement rather than a nicety. A plan whose parts do not add up to its total is not a plan, however accurate any individual part of it is.

Step 2 · Cleaning

d = raw.drop_duplicates(subset=["month", "store_id"]).copy()
print(f"raw rows                    {len(raw):,}")
print(f"after the duplicate extract {len(d):,}   (August 2022 ran twice)")

d["region"] = d.region.str.strip().str.title()
print(f"region spellings            {raw.region.nunique()} in the file, {d.region.nunique()} after normalizing")

d["month"] = pd.to_datetime(d.month)
bad = d.units <= 0
print()
print("NOT SALES, AND NOT ZERO EITHER")
print(d.loc[bad, ["month", "store_id", "units"]].to_string(index=False))
d.loc[bad, "units"] = np.nan
d = d.sort_values(["store_id", "month"])
d["units"] = d.groupby("store_id").units.transform(lambda s: s.interpolate().bfill().ffill())
d["units"] = d.units.round().astype(int)
print()
print(f"repaired by interpolating each store's own series: {int(bad.sum())} months")
print(f"analysis frame: {d.store_id.nunique()} stores, {d.region.nunique()} regions, "
      f"{d.month.nunique()} months, {d.month.min():%b %Y} to {d.month.max():%b %Y}")

raw rows                    1,360
after the duplicate extract 1,344   (August 2022 ran twice)
region spellings            6 in the file, 4 after normalizing

NOT SALES, AND NOT ZERO EITHER
     month store_id  units
2020-06-01      SO2      0
2021-03-01      NO4   -412
2023-09-01      WE3      0

repaired by interpolating each store's own series: 3 months
analysis frame: 16 stores, 4 regions, 84 months, Jan 2019 to Dec 2025

A till failure is not a month of zero sales, and a stocktake is not a sale. Both would be read by any forecasting method as genuine collapses in demand, and both would drag the store's trend estimate down for the rest of the history. Interpolating within the store's own series is the mildest repair that removes the fiction, and it is disclosed rather than silent.

Step 3 · First look

Sixteen series that will be forecast individually, four sums of them, and one sum of those. Before any of that, what the chain actually looks like.

w = d.pivot(index="month", columns="store_id", values="units").sort_index()
STORES = list(w.columns)
region_of = d.drop_duplicates("store_id").set_index("store_id").region

ann = w.iloc[-12:].sum()
print(f"{len(STORES)} stores, {len(REG)} regions, {len(w)} months")
print(f"  store size, last twelve months   {ann.min():,} to {ann.max():,} units "
      f"(median {ann.median():,.0f})")
print(f"  chain, last twelve months        {ann.sum():,} units")
seas = w.sum(1).groupby(w.index.month).mean()
print(f"  seasonality                      December runs {seas.max()/seas.mean()-1:+.0%} against the "
      f"average month, February {seas.min()/seas.mean()-1:+.0%}")

def trend(x):
    return np.polyfit(np.arange(len(x)), np.log(x), 1)[0]*12
tr = pd.Series({s: trend(w[s]) for s in STORES})
# a store has to be compared with its NEIGHBORS, not with a regional average it is
# itself a large part of
peers = pd.Series({s: trend(w[[o for o in STORES
                               if region_of[o] == region_of[s] and o != s]].sum(1))
                   for s in STORES})
div = (tr - peers).abs().sort_values(ascending=False)
print(f"  annual trend by store            {tr.min():+.1%} to {tr.max():+.1%}")
print("  furthest from their own neighbors:")
for s in div.index[:3]:
    print(f"      {s}  {tr[s]:+.1%} a year, against {peers[s]:+.1%} for the rest of {region_of[s]}")

16 stores, 4 regions, 84 months
  store size, last twelve months   31,944 to 204,664 units (median 68,514)
  chain, last twelve months        1,296,160 units
  seasonality                      December runs +35% against the average month, February -16%
  annual trend by store            -13.3% to +18.3%
  furthest from their own neighbors:
      SO3  -13.3% a year, against +6.1% for the rest of South
      WE2  +18.3% a year, against -0.8% for the rest of West
      SO4  +10.4% a year, against -2.9% for the rest of South

What this shows. The chain is strongly seasonal and the stores are not interchangeable. Annual trends run from thirteen percent down to eighteen percent up, and the two extremes are pulling away from their own neighbors: WE2 is growing at eighteen percent a year while the rest of the West edges down, and SO3 is falling at thirteen while the rest of the South grows at six. Hold on to those two. A method that splits one national forecast down by historical share cannot represent either of them, and Step 7 shows how badly.

fig, axes = plt.subplots(2, 2, figsize=(12.6, 8.2))

ax = axes[0, 0]
tot = w.sum(1)
ax.plot(tot.index, tot.values, color=SK, lw=2.2)
ax.set_ylabel("units per month"); ax.set_title("The chain, seven years")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x/1000:,.0f}k"))

ax = axes[0, 1]
idx = w / w.iloc[:12].mean()
for s in STORES:
    if s in ("WE2", "SO3"): continue
    ax.plot(idx.index, idx[s].values, color=MUT, lw=0.9, alpha=0.55)
ax.plot(idx.index, idx["WE2"].values, color=GD, lw=2.4, label="WE2, pulling away upward")
ax.plot(idx.index, idx["SO3"].values, color=RD, lw=2.4, label="SO3, pulling away downward")
ax.axhline(1, color=DK, lw=1.2, ls="--")
ax.set_ylabel("units, indexed to each store's first year")
ax.set_title("Sixteen stores, and two going the other way")
ax.legend(fontsize=9, loc="upper left")

ax = axes[1, 0]
mi = w.sum(1).groupby(w.index.month)
shape = (mi.mean()/mi.mean().mean() - 1)*100
ax.bar(range(1, 13), shape.values, color=[GD if v > 0 else MUT for v in shape.values], width=0.66)
for m_, v in zip(range(1, 13), shape.values):
    ax.text(m_, v + (1.4 if v > 0 else -3.2), f"{v:+.0f}", ha="center", fontsize=9,
            fontweight="bold", color=DK)
ax.axhline(0, color=DK, lw=1.4)
ax.set_xticks(range(1, 13)); ax.set_xticklabels(list("JFMAMJJASOND"))
ax.set_ylim(shape.min()*1.5, shape.max()*1.25)
ax.set_ylabel("percent against the average month")
ax.set_title("December carries the year")

ax = axes[1, 1]
lv = {"Chain (1 series)": [w.sum(1)],
      "Region (4)": [w[[s for s in STORES if region_of[s] == r]].sum(1) for r in REG],
      "Store (16)": [w[s] for s in STORES]}
names, vals = [], []
for k, series in lv.items():
    names.append(k)
    vals.append(np.mean([np.log(s).diff(12).dropna().std() for s in series])*100)
ax.bar(range(3), vals, color=[SK, LT, MUT], width=0.6)
for i, v in enumerate(vals):
    ax.text(i, v + 0.35, f"{v:.1f}%", ha="center", fontsize=10.5, fontweight="bold", color=DK)
indep = vals[2]/np.sqrt(len(STORES))          # what sixteen independent stores would give
ax.axhline(indep, color=RD, lw=2.0, ls="--")
ax.text(2.42, indep + 0.35, f"if the sixteen were independent: {indep:.1f}%",
        ha="right", fontsize=9.5, fontweight="bold", color=RD)
ax.set_xticks(range(3)); ax.set_xticklabels(names)
ax.set_ylim(0, max(vals)*1.22)
ax.set_ylabel("sd of year-on-year change (%)")
ax.set_title("Adding stores up smooths them, only so far")
ax.grid(axis="x", alpha=0)
plt.tight_layout(); plt.show()

Top left: a clean seasonal series with a mild upward drift, the kind of thing a forecasting method handles well. Top right: the same chain from underneath, and it is not one thing. Most stores drift gently; WE2 nearly doubles and SO3 falls away, and each is in the region least able to explain it.

Bottom left: December runs about a third above an average month, February well below. Bottom right: the season differenced out, so what is left is what genuinely differs between stores. Sixteen stores at 17.6 percent become a chain at 5.6, and the dashed line is where sixteen completely independent stores would have landed. The gap between the bar and the line is what the stores share: weather, national promotions, the economy. It is small here, which means most of the store-to-store variation really does cancel when you add them up. Remember that when Step 10 asks why bottom-up won.

Step 4 · The hierarchy, written down

Every method below is a choice of how to get from sixteen numbers to twenty-one. Writing the structure as a matrix makes that choice explicit, and makes three of the four methods one line of code each.

NAMES = ["Total"] + REG + STORES
S = np.zeros((len(NAMES), len(STORES)))
S[0, :] = 1
for i, r in enumerate(REG):
    S[1+i, [STORES.index(s) for s in STORES if region_of[s] == r]] = 1
S[1+len(REG):, :] = np.eye(len(STORES))

Y = S @ w.values.T                       # 21 series x 84 months, every level
print(f"summing matrix S: {S.shape[0]} series built from {S.shape[1]} stores")
print(f"  row 0 is the chain, rows 1 to {len(REG)} the regions, the rest the stores themselves")
print()
chk = np.abs(Y[0] - Y[1:1+len(REG)].sum(0)).max()
print(f"in the DATA the levels agree exactly: largest disagreement {chk:.0f} units")
print("Everything that follows is about the forecasts, which will not.")

summing matrix S: 21 series built from 16 stores
  row 0 is the chain, rows 1 to 4 the regions, the rest the stores themselves

in the DATA the levels agree exactly: largest disagreement 0 units
Everything that follows is about the forecasts, which will not.

Step 5 · Forecast every level on its own, and watch it come apart

The obvious thing, and the thing most organizations actually do: each team forecasts the series it is responsible for. Holt-Winters on each of the twenty-one, fitted on the first six years, asked for the next twelve months.

H, ORIGIN = 12, 72

def base_forecast(series, h):
    m = ExponentialSmoothing(series, trend="add", seasonal="add", seasonal_periods=12,
                             initialization_method="estimated").fit()
    return np.asarray(m.forecast(h)), np.asarray(m.fittedvalues)

train = Y[:, :ORIGIN]
base = np.zeros((len(NAMES), H)); resid = []
for i in range(len(NAMES)):
    f, fitted = base_forecast(train[i].astype(float), H)
    base[i] = f; resid.append(train[i] - fitted)
resid = np.vstack(resid)

nat  = base[0].sum()
regs = base[1:1+len(REG)].sum()
stor = base[1+len(REG):].sum()
print(f"forecast for {w.index[ORIGIN]:%b %Y} to {w.index[-1]:%b %Y}, made from {w.index[ORIGIN-1]:%b %Y}")
print()
print(f"  the board is shown            {nat:12,.0f}   (the chain series, forecast directly)")
print(f"  the four regions sum to       {regs:12,.0f}   ({regs/nat-1:+.2%})")
print(f"  the sixteen stores sum to     {stor:12,.0f}   ({stor/nat-1:+.2%})")
print(f"  widest disagreement           {max(nat,regs,stor)-min(nat,regs,stor):12,.0f} units")

forecast for Jan 2025 to Dec 2025, made from Dec 2024

  the board is shown               1,278,691   (the chain series, forecast directly)
  the four regions sum to          1,273,167   (-0.43%)
  the sixteen stores sum to        1,284,091   (+0.42%)
  widest disagreement                 10,924 units

Under one percent, and it still cannot be shipped. The gap is small enough that nobody would call any of the three forecasts wrong, and that is exactly the trap: there is no error to find and no model to fix. The three numbers disagree because they were produced by three separate fits to three different series, and nothing in the procedure ever required them to agree.

A plan is a set of commitments. If the regional targets sum to something other than the number the board was given, then somebody is going to be held to a figure that the plan does not contain. The size of the gap is beside the point. Coherence is a property the output has to have, and only one of the four methods below achieves it by accident.

Step 6 · Bottom-up

Forecast the sixteen stores and add them up. Coherent by construction, and it keeps every store's own history, including the two that are going the wrong way.

G_bu = np.hstack([np.zeros((len(STORES), len(NAMES)-len(STORES))), np.eye(len(STORES))])
f_bu = S @ (G_bu @ base)
print(f"bottom-up chain total   {f_bu[0].sum():12,.0f}")
print(f"the direct chain fit    {base[0].sum():12,.0f}   (discarded: only the store forecasts are used)")
print()
print("Coherent by construction, because every level is now a sum of the same sixteen numbers.")
print("The cost is that the chain forecast is an accumulation of sixteen individual errors.")

bottom-up chain total      1,284,091
the direct chain fit       1,278,691   (discarded: only the store forecasts are used)

Coherent by construction, because every level is now a sum of the same sixteen numbers.
The cost is that the chain forecast is an accumulation of sixteen individual errors.

Step 7 · Top-down, and the two stores it cannot see

Forecast the chain, which is the smoothest series in the file, then split it down using each store's historical share. Also coherent, and it is the method that fails hardest here.

props = train[1+len(REG):, :].sum(1) / train[0, :].sum()
f_td = S @ (props[:, None] * base[0][None, :])
print("historical share of the chain, six years")
for s, p in sorted(zip(STORES, props), key=lambda kv: -kv[1])[:3]:
    print(f"  {s}  {p:.4f}")
print("  ...")
print()
print("WHAT A FIXED SHARE COSTS THE TWO STORES FROM STEP 3")
print(f"  {'store':6s}{'six-year share':>16s}{'share last year':>17s}"
      f"{'top-down plan':>15s}{'against last year':>19s}")
recent = w.iloc[-12:].sum()
for s in ["WE2", "SO3"]:
    i = NAMES.index(s)
    print(f"  {s:6s}{props[STORES.index(s)]:>16.1%}{recent[s]/recent.sum():>17.1%}"
          f"{f_td[i].sum():>15,.0f}{f_td[i].sum()/recent[s]-1:>18.1%}")
print()
print("The share is an average over six years. WE2 has grown through that window and SO3")
print("has shrunk, so the average describes neither of them now, and the error is a level")
print("shift rather than a wrong growth rate.")

historical share of the chain, six years
  WE2  0.1007
  SO4  0.0851
  NO4  0.0799
  ...

WHAT A FIXED SHARE COSTS THE TWO STORES FROM STEP 3
  store   six-year share  share last year  top-down plan  against last year
  WE2              10.1%            15.8%        128,722            -37.1%
  SO3               4.5%             2.5%         57,251             79.2%

The share is an average over six years. WE2 has grown through that window and SO3
has shrunk, so the average describes neither of them now, and the error is a level
shift rather than a wrong growth rate.

A share is a constant, and these two stores are not. Top-down is not merely less accurate for WE2 and SO3, it is structurally incapable of representing them. WE2 has climbed from a tenth of the chain to a sixth while SO3 has fallen away, and averaging over six years splits the difference on both, so the plan cuts a growing store by a third and hands a shrinking one a large increase. The rest of the chain is well enough behaved that the method still looks defensible in aggregate, which is how it survives in practice.

Step 8 · Reconciliation

Both methods above throw information away: bottom-up ignores the chain fit, top-down ignores the store fits. Reconciliation keeps all twenty-one forecasts and finds the coherent set closest to them, which is a projection and therefore one line of linear algebra.

OLS reconciliation treats every series as equally reliable. MinT weights them by how large their errors have been and by how those errors move together, which is the information OLS discards.

G_ols = np.linalg.solve(S.T @ S, S.T)
f_ols = S @ (G_ols @ base)

W = np.cov(resid)
LAM = 0.35                                  # shrink toward the diagonal; 21 series, 72 months
W_sh = LAM*np.diag(np.diag(W)) + (1-LAM)*W
Wi = np.linalg.inv(W_sh + 1e-6*np.eye(len(NAMES)))
G_mint = np.linalg.solve(S.T @ Wi @ S, S.T @ Wi)
f_mint = S @ (G_mint @ base)

print("EVERY METHOD IS NOW COHERENT")
print(f"  {'method':12s}{'chain total':>14s}{'regions sum to':>17s}{'stores sum to':>16s}")
for nm, f in [("independent", base), ("bottom-up", f_bu), ("top-down", f_td),
              ("OLS", f_ols), ("MinT", f_mint)]:
    print(f"  {nm:12s}{f[0].sum():>14,.0f}{f[1:1+len(REG)].sum():>17,.0f}{f[1+len(REG):].sum():>16,.0f}")
print()
print(f"how far each moved the original forecasts (mean absolute change, units per series-month)")
for nm, f in [("bottom-up", f_bu), ("top-down", f_td), ("OLS", f_ols), ("MinT", f_mint)]:
    print(f"  {nm:12s}{np.abs(f-base).mean():>10,.0f}")

EVERY METHOD IS NOW COHERENT
  method         chain total   regions sum to   stores sum to
  independent      1,278,691        1,273,167       1,284,091
  bottom-up        1,284,091        1,284,091       1,284,091
  top-down         1,278,691        1,278,691       1,278,691
  OLS              1,277,896        1,277,896       1,277,896
  MinT             1,280,324        1,280,324       1,280,324

how far each moved the original forecasts (mean absolute change, units per series-month)
  bottom-up          103
  top-down         1,266
  OLS                 95
  MinT                93

Only the first row fails to add up, and it is the one nobody reconciled. The other four are coherent by construction or by projection.

The last block is the honest summary of what reconciliation does: it is a nudge, not a rebuild. MinT moves the original forecasts less than bottom-up or top-down do, because it is looking for the nearest coherent set rather than imposing a structure. Whether that nudge is worth having is an empirical question, and the only way to answer it is to backtest.

Step 9 · Backtest at every level

Five origins, twelve months ahead from each, every method refitted at every origin. Scored with MASE, which divides by the seasonal naive error on that series' own training data, so a store doing four thousand units a month and a chain doing a hundred thousand can appear in the same table.

ORIGINS = [48, 54, 60, 66, 72]
rows = []
for o in ORIGINS:
    tr = Y[:, :o]
    b = np.zeros((len(NAMES), H)); rs = []
    for i in range(len(NAMES)):
        f, fitted = base_forecast(tr[i].astype(float), H)
        b[i] = f; rs.append(tr[i] - fitted)
    Wk = np.cov(np.vstack(rs))
    Wik = np.linalg.inv(LAM*np.diag(np.diag(Wk)) + (1-LAM)*Wk + 1e-6*np.eye(len(NAMES)))
    Gm = np.linalg.solve(S.T @ Wik @ S, S.T @ Wik)
    Go = np.linalg.solve(S.T @ S, S.T)
    pr = tr[1+len(REG):, :].sum(1) / tr[0, :].sum()
    cand = {"independent": b, "bottom-up": S @ b[1+len(REG):],
            "top-down": S @ (pr[:, None] * b[0][None, :]),
            "OLS": S @ (Go @ b), "MinT": S @ (Gm @ b)}
    act = Y[:, o:o+H]
    scale = np.array([np.mean(np.abs(tr[i][12:] - tr[i][:-12])) for i in range(len(NAMES))])
    for nm, f in cand.items():
        e = np.abs(f - act)
        for i in range(len(NAMES)):
            lvl = "Chain" if i == 0 else ("Region" if i <= len(REG) else "Store")
            rows.append(dict(origin=o, method=nm, level=lvl, series=NAMES[i],
                             mase=e[i].mean()/scale[i]))
B = pd.DataFrame(rows)
ORDER = ["independent", "bottom-up", "top-down", "OLS", "MinT"]
tbl = B.groupby(["level", "method"]).mase.mean().unstack().reindex(["Chain", "Region", "Store"])[ORDER]
print(f"MASE by level, {len(ORIGINS)} origins x {H} months (below 1 beats seasonal naive)")
display(tbl.round(3))

MASE by level, 5 origins x 12 months (below 1 beats seasonal naive)

method	independent	bottom-up	top-down	OLS	MinT
level
Chain	0.840	0.747	0.840	0.839	0.802
Region	0.799	0.782	0.999	0.796	0.788
Store	0.792	0.792	1.198	0.798	0.796

Read the top-down column first. It is the only one that is worse than a seasonal naive rule anywhere, and at store level it is roughly half again as bad as everything else. That is the cost of expressing sixteen different businesses as one number times sixteen constants.

Then read across the other four. They are close. Bottom-up is the best single choice on this chain at every level, MinT is a step behind it, and OLS trails MinT everywhere, which is the argument for using the error covariance rather than assuming it away.

worst = B[B.method == "top-down"].groupby("series").mase.mean().sort_values(ascending=False)
bu = B[B.method == "bottom-up"].groupby("series").mase.mean()
print("WHERE TOP-DOWN GOES WRONG, worst five series")
print(f"  {'series':8s}{'top-down':>10s}{'bottom-up':>11s}{'penalty':>10s}")
for s in worst.index[:5]:
    print(f"  {s:8s}{worst[s]:>10.2f}{bu[s]:>11.2f}{worst[s]/bu[s]:>9.1f}x")

WHERE TOP-DOWN GOES WRONG, worst five series
  series    top-down  bottom-up   penalty
  WE2           3.18       0.74      4.3x
  SO3           2.61       0.69      3.8x
  SO4           1.64       0.71      2.3x
  West          1.26       0.85      1.5x
  WE1           1.23       0.88      1.4x

The two stores identified in Step 3 are at the top of that list, and they are there for the reason the picture predicted rather than for anything the backtest discovered on its own. This is what it looks like when a diagnostic done before any modeling tells you where the model will fail.

Step 10 · Which number to plan on

# only the coherent methods are candidates: "independent" cannot be shipped as a plan
COH = ["bottom-up", "top-down", "OLS", "MinT"]
ct = tbl[COH]
best, gap = ct.idxmin(axis=1), ct["MinT"] - ct.min(axis=1)
print(f"  {'level':8s}{'best coherent':>15s}{'its MASE':>11s}{'MinT':>9s}{'MinT gives up':>15s}")
for lvl in ["Chain", "Region", "Store"]:
    print(f"  {lvl:8s}{best[lvl]:>15s}{ct.loc[lvl].min():>11.3f}"
          f"{ct.loc[lvl,'MinT']:>9.3f}{gap[lvl]:>15.3f}")
print()
print(f"choosing top-down instead would cost "
      f"{tbl.loc['Store','top-down']/ct.loc['Store'].min()-1:.0%} at store level.")
print(f"choosing MinT instead of the winner costs at most {gap.max():.3f} MASE, at the chain level.")

  level     best coherent   its MASE     MinT  MinT gives up
  Chain         bottom-up      0.747    0.802          0.055
  Region        bottom-up      0.782    0.788          0.007
  Store         bottom-up      0.792    0.796          0.005

choosing top-down instead would cost 51% at store level.
choosing MinT instead of the winner costs at most 0.055 MASE, at the chain level.

Bottom-up won, and that could not have been known in advance. It won because these sixteen stores are individually forecastable and their errors partly cancel. On a chain of two hundred small stores, or on weekly data where each store is much noisier, the same backtest routinely goes the other way and the aggregate forecast carries the day.

That is the case for reconciliation, and it is a modest one. MinT was never the best method at any level here and never more than 0.055 MASE behind the one that was, while requiring no judgment about which level to trust. It is insurance rather than an improvement, and it is cheap. What the backtest does establish beyond doubt is which method not to use.

Step 11 · The pictures

fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.4))

ax = axes[0]
lv = ["Chain", "Region", "Store"]
xp = np.arange(len(lv)); wd = 0.16
cols = {"independent": MUT, "bottom-up": SK, "top-down": RD, "OLS": LT, "MinT": GD}
for j, m_ in enumerate(ORDER):
    v = [tbl.loc[l, m_] for l in lv]
    ax.bar(xp + (j-2)*wd, v, width=wd, color=cols[m_], label=m_)
ax.axhline(1, color=DK, lw=1.6, ls="--")
ax.text(2.42, 1.02, "seasonal naive", ha="right", fontsize=9, fontweight="bold", color=DK)
ax.set_xticks(xp); ax.set_xticklabels(lv)
ax.set_ylabel("MASE (lower is better)"); ax.set_ylim(0, 1.35)
ax.set_title("Top-down is the only one that loses to naive")
ax.legend(fontsize=8.5, ncol=3, loc="upper left"); ax.grid(axis="x", alpha=0)

ax = axes[1]
sel = ["WE2", "SO3"]
xp2 = np.arange(len(sel)); wd2 = 0.34
for j, m_ in enumerate(["bottom-up", "top-down"]):
    v = [B[(B.method == m_) & (B.series == s)].mase.mean() for s in sel]
    ax.bar(xp2 + (j-0.5)*wd2, v, width=wd2, color=cols[m_], label=m_)
    for x_, val in zip(xp2 + (j-0.5)*wd2, v):
        ax.text(x_, val + 0.06, f"{val:.2f}", ha="center", fontsize=10,
                fontweight="bold", color=DK)
ax.axhline(1, color=DK, lw=1.6, ls="--")
ax.set_xticks(xp2)
ax.set_xticklabels(["WE2, share up\na tenth to a sixth", "SO3, share falling\naway"])
ax.set_ylabel("MASE"); ax.set_ylim(0, 3.7)
ax.set_title("The two stores a share cannot describe")
ax.legend(fontsize=9, loc="upper left"); ax.grid(axis="x", alpha=0)
plt.tight_layout(); plt.show()

Left: four of the five methods sit close together and comfortably inside the naive benchmark. Top-down is outside it at store level, which is a way of saying that for individual stores you would have been better off assuming next December equals last December.

Right: the two stores from the first look, scored. Top-down is between three and four times worse on both, and the reason was visible in Step 3 rather than discovered here.

Step 12 · What this does not settle

Coherence is not accuracy. Every method after Step 5 produces numbers that add up, including the worst one in the table. Adding up is a property of the output, and it says nothing at all about whether the output is any good.

The winner is a property of this chain. Sixteen stores, monthly, with partly independent errors. Change any of those and the ranking can invert, which is why the backtest is the deliverable and not the model.

One structure, one moment. The hierarchy here is geographic and fixed. Real chains open and close stores, move them between regions and sell through channels that cut across geography, so a store can belong to two hierarchies at once. Reconciling across several groupings at the same time is a harder problem than the one worked here.

Nothing above produced an interval. Every number in this notebook is a point, and a plan is built on a range. Prediction intervals for a reconciled forecast are not the intervals of the base forecasts, because the projection mixes them, and coverage has to be measured rather than assumed. That is the subject of the last capstone in this part.

←Back to the bookChapter 199 · Hierarchical Forecasting: Store, Region, National