Advanced Forecasting

Run the code. Part 1 builds independent base forecasts for the network total and three routes (synthetic; top-down is given the true shares, a best case for it). Check that the base forecasts disagree, that bottom-up, top-down and WLS-reconciled forecasts add up, and which level each method improves. Then change sd_base so the total is the most accurate level and see how the ranking changes. Part 2 simulates 300 intermittent series for an ancillary sold on about one flight in four. Compare mean forecasts with the true rate of 0.625 per flight, and compare MAE with RMSE for the zero forecast. Finally write one sentence on which error measure you would use to choose an intermittent-demand method for stocking decisions, and why.

import numpy as np

rng = np.random.default_rng(260)

# ---- 1. A three-route hierarchy: total = ESB + ADB + AYT (synthetic weekly bookings) -------------
S = np.array([[1, 1, 1], [1, 0, 0], [0, 1, 0], [0, 0, 1]], float)  # rows: total, ESB, ADB, AYT
truth = np.array([900.0, 700.0, 450.0])                              # route means
reps = 2000
y = truth + rng.normal(0, [30, 25, 20], (reps, 3))                   # realised route demand
actual = y @ S.T                                                     # columns: total, ESB, ADB, AYT
sd_base = np.array([55.0, 40.0, 30.0, 25.0])                         # base-forecast error sd per level
base = actual + rng.normal(0, sd_base, actual.shape)                 # independent base forecasts per level

w_inv = np.diag(1 / sd_base**2)                                      # diagonal W: weight by error variance
P = np.linalg.solve(S.T @ w_inv @ S, S.T @ w_inv)                    # (S'W^-1 S)^-1 S'W^-1
reconciled = base @ (S @ P).T                                        # y_tilde = S P y_hat
bottom_up = np.column_stack([base[:, 1:].sum(1), base[:, 1:]])
top_down = np.column_stack([base[:, 0], base[:, [0]] * (truth / truth.sum())])
print("method          RMSE total  RMSE routes  mean |total - sum(routes)|")
for name, f in [("base", base), ("bottom-up", bottom_up), ("top-down", top_down), ("WLS reconciled", reconciled)]:
    err = f - actual
    gap = np.abs(f[:, 0] - f[:, 1:].sum(1)).mean()
    print(f"{name:15s} {np.sqrt((err[:, 0] ** 2).mean()):9.1f} {np.sqrt((err[:, 1:] ** 2).mean()):11.1f} {gap:12.1f}")

# ---- 2. Intermittent demand: an ancillary bought on about 1 flight in 4 ---------------------------
def croston(d, alpha, sba=False):
    nz = np.flatnonzero(d)
    size, gap, last = d[nz[0]], nz[0] + 1, nz[0]
    out = np.zeros(len(d))
    for t in range(len(d)):
        out[t] = size / gap * (1 - alpha / 2 if sba else 1)          # forecast made before seeing d[t]
        if d[t] > 0:
            size += alpha * (d[t] - size)
            gap += alpha * ((t - last) - gap)
            last = t
    return out

p, T, burn = 0.25, 300, 100
rate = p * 2.5                                                       # E[size] = 1 + 1.5
rows = {"zero": [], "mean so far": [], "croston": [], "sba": []}
for _ in range(300):
    d = rng.binomial(1, p, T) * (1 + rng.poisson(1.5, T))
    if d.sum() == 0:
        continue
    fc = {"zero": np.zeros(T), "mean so far": np.array([d[:t].mean() if t else 0.0 for t in range(T)]),
          "croston": croston(d, 0.1), "sba": croston(d, 0.1, True)}
    for k, f in fc.items():
        e = d[burn:] - f[burn:]
        rows[k].append((f[burn:].mean(), np.abs(e).mean(), np.sqrt((e ** 2).mean())))
print(f"\ntrue demand rate {rate:.3f} per flight")
for k, v in rows.items():
    m, mae, rmse = np.mean(v, axis=0)
    print(f"{k:12s} mean forecast {m:.3f}  MAE {mae:.2f}  RMSE {rmse:.2f}")

A logistics planner asks for "a 95% range" on next month's parcel volume at a sorting hub and wants to contract trucks and temporary staff at the top of it. How do you respond?