ai_bust_probability.py
raw file ยท 1,143 lines
"""
ai_bust_probability.py
======================
Part I of the working paper: HOW LIKELY is an AI bust? Part III: expected damage and an odds tracker.
Two crash definitions, reported separately:
* MARKET CRASH - AI-linked stocks (proxy: the PHLX Semiconductor index, SOX) fall >= 40% from a peak
(the crash definition used by Greenwood, Shleifer & You, "Bubbles for Fama", 2019).
* ECONOMIC BUST - end-customer AI spending falls >= 15% below the PLAN path at any quarter-end.
PLAN = revenue keeps pace with the AI capital stock (constant revenue per dollar of AI
capital, i.e. returns on the buildout stop improving AND fall 15% from today's level).
15% is where the firm-level model (Model F) produces neocloud failures with certainty.
Four independent estimates for each definition, as of October 2026, cumulative to end-2027 / 2028 / 2029:
1. FUNDAMENTALS (Model G) - simulated revenue paths: hypergrowth that decays at historical rates,
growth noise, and random demand stalls; compared with the plan path.
2. HISTORY - economic bust: past privately financed investment booms (onset -> bust timing);
market crash: Greenwood-Shleifer-You crash odds after >=100% sector run-ups.
3. MARKET PRICES - market crash: option-implied probability of a 40% fall (barrier formula);
economic bust: credit-spread-implied neocloud default odds, translated to a bust.
4. INDICATORS - economic bust: Greenwood-Hanson-Shleifer-Sorensen "R-zone" (credit growth + price run-up);
market crash: GSY crash characteristics (volatility, issuance, acceleration) present.
Each method is a Monte Carlo over its own uncertain inputs. Methods are combined by weighted
logarithmic pooling (weighted average of log-odds), sampled jointly to give a range.
Every input is labelled [REAL] (sourced, Sept/Oct 2026), [HIST] (historical record) or [ASSUMED].
Run: python3 ai_bust_probability.py (~3.5 min including the Model D and F reruns and sensitivities)
python3 ai_bust_probability.py --fast (~30 s; skips the Model F rerun, writes probability_results_fast.json)
Needs numpy and scipy; the expected-damage step imports ai_bust_models.py and ai_bust_abm.py.
"""
from __future__ import annotations
import argparse
import json
import math
from dataclasses import replace
import numpy as np
from scipy.stats import norm
N = 40_000 # Monte Carlo draws per method
RNG = np.random.default_rng(20261005)
# Decimal calendar time: 2026.75 = 1 Oct 2026; the END of 2027 is t = 2028.0.
HORIZONS = {"end-2027": 2028.0, "end-2028": 2029.0, "end-2029": 2030.0} # cumulative probabilities
NOW = 2026.75 # as of early October 2026
BUST_SHORTFALL = 0.15 # economic-bust threshold
def logit(p):
p = np.clip(p, 1e-6, 1 - 1e-6)
return np.log(p / (1 - p))
def inv_logit(x):
return 1 / (1 + np.exp(-x))
# =============================================================================
# 1. FUNDAMENTALS - MODEL G: revenue paths vs the plan path
# =============================================================================
# Plan path: AI capital stock = cumulative gross AI capex since 2024, $T, at each YEAR-END (Model A):
# end-2025 1.05, end-2026 2.05 [REAL-anchored]; end-2027 3.40, end-2028 5.00 (Goldman Sachs capex path),
# end-2029 6.80 (2029 capex ~$1.8T, extrapolated) [ASSUMED]. Keys are decimal times (end-2026 = 2027.0).
CAPITAL = {2026.0: 1.05, 2027.0: 2.05, 2028.0: 3.40, 2029.0: 5.00, 2030.0: 6.80}
def plan_index(t):
"""Revenue the plan requires at time t, relative to Q4-2026 (log-linear in the capital stock)."""
ks = sorted(CAPITAL)
t = np.asarray(t, dtype=float)
logk = np.interp(t, ks, [math.log(CAPITAL[k]) for k in ks])
return np.exp(logk - math.log(CAPITAL[2026.0]))
def model_g_legacy(n=N, rng=RNG, overrides=None, return_paths=False):
"""Simulate end-customer AI spending R(t) quarterly from October 2026, normalised so R/plan = 1 today:
today's revenue per dollar of AI capital is the plan's reference level.
Annual log growth g follows, per quarter:
g <- g_inf + phi^(1/4) (g - g_inf) + sigma_g * 0.5 * eps (decaying hypergrowth + growth noise)
log R <- log R + g/4 + sigma_l * eta (level noise, run-rate measurement)
Demand stalls arrive with hazard lam per year: growth drops to g_stall for D quarters, then restarts at g_inf.
Priors (per path):
G0 = e^g0 - 1 ~ Tri(0.5, 0.9, 1.5) [REAL-anchored] current growth: OpenAI booked revenue +94% annualised
(Q1->Q2 2026), hyperscaler AI run-rates ~x2.4/yr, Google Cloud +63%
G_inf ~ U(0.08, 0.25) [ASSUMED] long-run growth of a general-purpose technology's spend
phi (annual) ~ U(0.40, 0.80) [HIST] persistence of excess growth: smartphones 2010-15 ~0.75-0.8,
cloud IaaS 2015-22 ~0.75; faster decay from higher starting growth
sigma_g ~ U(0.08, 0.20) [ASSUMED] dispersion of annual growth
lam ~ U(0.05, 0.20) /yr [HIST] semiconductor demand downturns ~every 4-5 yrs; enterprise IT
stalls mostly in recessions (NY Fed yield-curve model: 13.9% for
the 12 months to Aug-2027)
g_stall ~ U(-0.20, 0.05) [ASSUMED] annualised log growth during a stall
D ~ U{2..6} quarters
"""
o = overrides or {}
def draw(name, f):
v = o.get(name)
return np.full(n, float(v)) if v is not None else f()
g0 = np.log1p(draw("G0", lambda: rng.triangular(0.5, 0.9, 1.5, n)))
ginf = np.log1p(draw("G_inf", lambda: rng.uniform(0.08, 0.25, n)))
phi = draw("phi", lambda: rng.uniform(0.40, 0.80, n))
sig_g = draw("sigma_g", lambda: rng.uniform(0.08, 0.20, n))
lam = draw("lam", lambda: rng.uniform(0.05, 0.20, n))
sig_l = 0.02
times = np.round(np.arange(NOW + 0.25, 2030.0 + 1e-9, 0.25), 2) # quarter-ends, end-2026 ... end-2029
T = len(times)
logR = np.zeros(n)
g = g0.copy()
stall_left = np.zeros(n, dtype=int)
g_stall = rng.uniform(-0.20, 0.05, n)
ratio = np.zeros((n, T)) # R / plan at each quarter-end
stalled = np.zeros(n, dtype=bool)
phi_q = phi ** 0.25
for j, t in enumerate(times):
# stall arrivals
start = (stall_left == 0) & (rng.random(n) < 1 - np.exp(-lam * 0.25))
stall_left = np.where(start, rng.integers(2, 7, n), stall_left)
g_stall = np.where(start, rng.uniform(-0.20, 0.05, n), g_stall)
stalled |= start
in_stall = stall_left > 0
g_eff = np.where(in_stall, g_stall, g)
logR = logR + g_eff / 4 + sig_l * rng.standard_normal(n)
# growth dynamics (outside stalls); after a stall ends growth restarts at g_inf
g = np.where(in_stall, ginf, ginf + phi_q * (g - ginf) + sig_g * 0.5 * rng.standard_normal(n))
stall_left = np.where(in_stall, stall_left - 1, 0)
ratio[:, j] = np.exp(logR) / plan_index(t) * plan_index(NOW)
shortfall = 1 - ratio # >0 = below plan
out = {"times": times, "ratio": ratio, "shortfall": shortfall, "stalled": stalled,
"params": {"g0": g0, "ginf": ginf, "phi": phi, "lam": lam}}
return out
# -----------------------------------------------------------------------------
# MODEL G v2: tiered demand (flighty / sticky / sovereign) + a physical power ceiling on capex
# -----------------------------------------------------------------------------
# SPEC switches every upgrade on or off, so each can be ablated. All four off = the v1 model (up to Monte Carlo noise).
SPEC = {"tiers": True, # split demand into flighty vs sticky private spending
"flighty_hazard": True, # extra flighty-only "disillusionment" events (no macro recession needed)
"sovereign": True, # government/state-backed demand as a third, low-beta tier
"power": True} # energisable-GW ceiling on capex, plus a revenue cap from powered capacity
# Tier calibration (research brief, Oct 2026). [S] = sourced proxy, [I] = judgment.
# share of Q3-2026 end-customer AI spend that is flighty (consumer subs, weak-adoption seat copilots, POCs, unbudgeted
# token overrun): 38% (30-50%) [I from Menlo/a16z/OpenAI mix]; sovereign share 4-15% [I: strict ~4%, Nvidia-broad ~15%]
F_SHARE = (0.30, 0.38, 0.50) # Tri(lo, mode, hi) of PRIVATE spend that is flighty
W_SOV = (0.04, 0.15) # U: sovereign share of total
STICKY_OVER_FLIGHTY_GROWTH = (1.5, 3.0) # ratio of sticky to flighty current annual growth (arithmetic)
G_SOV0 = (0.6, 1.2, 2.2) # Tri: current sovereign spend growth (Nvidia: sovereign >3x y/y)
GINF = {"f": (0.02, 0.15), "s": (0.10, 0.28), "v": (0.10, 0.25)}
PHI = {"f": (0.30, 0.65), "s": (0.45, 0.85), "v": (0.50, 0.85)}
# Stall sensitivities (multipliers on a macro demand stall):
X_F = (0.20, 0.35, 0.50) # flighty level falls 20-50% over <=4 quarters (2002/2009 new-project spend, AI GRR/POC data)
KAPPA_S = (0.40, 0.80) # sticky growth slows by 40-80% of its pre-shock rate (2009 maintenance, 2022-23 usage cos.)
X_S = (0.0, 0.04, 0.15) # sticky level change in a stall 0 to -15%
KAPPA_V = (0.10, 0.35) # sovereign growth slows 10-35% in a private stall (Stargate UK/Norway: private offtakers move)
LAM_F = (0.05, 0.25) # flighty-only disillusionment events per year
XEV_F = (0.10, 0.30) # their size
LAM_V = (0.15, 0.35) # sovereign programme delays per year [30%/flagship/yr in the brief, diluted by breadth]
XEV_V = (0.05, 0.20) # level lost per delay
# Power ceiling (research brief): energisable new AI capacity, GW per year, global (low, central, high); $B per GW all-in.
GW_TRI = {2027: (14, 20, 27), 2028: (16, 23, 32), 2029: (18, 26, 38)}
COST_PER_GW = {2027: 50.0, 2028: 53.0, 2029: 55.0} # $B/GW central; common scalar Tri(0.84, 1, 1.2)
REFRESH_HEADROOM = (0.05, 0.20) # $T/yr: refreshing chips in halls already powered
STRANDED_SHARE = (0.20, 0.50) # share of over-ceiling spend that is bought anyway and parked
REV_HEADROOM = (0.05, 0.30) # revenue per powered $ can rise this much (utilisation, scarcity pricing)
def _tri_ppf(u, lo, mode, hi):
c = (mode - lo) / (hi - lo)
return np.where(u < c, lo + np.sqrt(u * (hi - lo) * (mode - lo)), hi - np.sqrt((1 - u) * (hi - lo) * (hi - mode)))
def _interp_knots(logK, t, knots=(2026.0, 2027.0, 2028.0, 2029.0, 2030.0)):
"""Log-linear interpolation of a per-path capital stock (n x 5 knots) at time t."""
k = np.asarray(knots)
i = int(np.clip(np.searchsorted(k, t, side="right") - 1, 0, len(k) - 2))
w = (t - k[i]) / (k[i + 1] - k[i])
return logK[:, i] * (1 - w) + logK[:, i + 1] * w
def power_capital(n, rng, on=True, overrides=None):
"""Capital stock at the 5 year-end knots: SPENT (what firms have bought) and POWERED (what can earn revenue).
Plan capex 2027-29 = diff of CAPITAL. Ceiling_y = GW_y x $/GW_y + refresh headroom. Spend above the ceiling is
deferred (never spent in the horizon) except a stranded share that is bought and parked, energised a year later."""
o = overrides or {}
ks = sorted(CAPITAL)
plan = np.array([CAPITAL[ks[i + 1]] - CAPITAL[ks[i]] for i in range(1, 4)]) # capex in 2027, 2028, 2029 ($T)
base = np.log([CAPITAL[k] for k in ks])
if not on:
lk = np.tile(base, (n, 1))
return lk, lk.copy(), {"ceiling": None}
rho = 0.85
z0 = rng.standard_normal(n)
cost_scale = _tri_ppf(rng.random(n), 0.84, 1.0, 1.2)
scale = float(o.get("power_scale", 1.0))
stranded = rng.uniform(*STRANDED_SHARE, n)
refresh = rng.uniform(*REFRESH_HEADROOM, n)
ceil = np.zeros((n, 3))
for j, y in enumerate((2027, 2028, 2029)):
z = rho * z0 + math.sqrt(1 - rho ** 2) * rng.standard_normal(n)
gw = _tri_ppf(norm.cdf(z), *GW_TRI[y])
ceil[:, j] = scale * (gw * COST_PER_GW[y] * cost_scale / 1000.0 + refresh)
spent_y = np.zeros((n, 3)); powered_y = np.zeros((n, 3)); parked = np.zeros(n)
for j in range(3):
within = np.minimum(plan[j], ceil[:, j]); excess = np.maximum(plan[j] - ceil[:, j], 0.0)
spent_y[:, j] = within + stranded * excess
powered_y[:, j] = within + parked # parked chips from last year get energised
parked = stranded * excess
Ks = np.concatenate([np.tile(np.exp(base[:2]), (n, 1)), CAPITAL[2027.0] + np.cumsum(spent_y, axis=1)], axis=1)
Kp = np.concatenate([np.tile(np.exp(base[:2]), (n, 1)), CAPITAL[2027.0] + np.cumsum(powered_y, axis=1)], axis=1)
return np.log(Ks), np.log(Kp), {"ceiling": ceil, "plan": plan, "excess_share": np.maximum(plan - ceil, 0) / plan}
def model_g(n=N, rng=RNG, overrides=None, spec=None):
"""Model G v2. Quarterly simulation of end-customer AI spending by tier from October 2026.
Tiers: flighty (consumer subs, weak-adoption seats, POCs, token overrun), sticky (embedded in products and
pipelines), sovereign (government / state-backed buyers). Each tier has its own growth, persistence and long-run
growth; a macro demand stall (hazard lam) hits them very differently:
flighty level falls x_f (20-50%) within <=4 quarters (2002/2009 new-project spend; AI-app GRR 23-45%)
sticky growth slows by kappa_s (40-80%), level -0..-15% (2009 maintenance; 2022-23 usage-based software)
sovereign growth slows by kappa_v (10-35%) (budget-funded; private-offtaker projects move)
After any macro stall growth restarts at the tier's long-run rate, as in v1.
Optional extras: flighty-only disillusionment events, sovereign programme delays, a power ceiling on capex.
With all four spec switches off this is a single-tier model like v1 (plus nothing else)."""
sp = {**SPEC, **(spec or {})}
o = overrides or {}
def draw(name, f):
v = o.get(name)
return np.full(n, float(v)) if v is not None else f()
G_agg = draw("G0", lambda: rng.triangular(0.5, 0.9, 1.5, n))
u_ginf = (np.asarray(o["G_inf"]) - 0.08) / 0.17 if "G_inf" in o else rng.random(n)
u_phi = (np.asarray(o["phi"]) - 0.40) / 0.40 if "phi" in o else rng.random(n)
u_ginf = np.broadcast_to(u_ginf, (n,)); u_phi = np.broadcast_to(u_phi, (n,))
sig_g = draw("sigma_g", lambda: rng.uniform(0.08, 0.20, n))
lam = draw("lam", lambda: rng.uniform(0.05, 0.20, n))
sig_l = 0.02
times = np.round(np.arange(NOW + 0.25, 2030.0 + 1e-9, 0.25), 2)
T = len(times)
tiers = ["f", "s", "v"] if sp["tiers"] and sp["sovereign"] else (["f", "s"] if sp["tiers"] else ["a"])
K = len(tiers)
# ---- shares and current growth by tier
if sp["tiers"]:
w_v = draw("w_v", lambda: rng.uniform(*W_SOV, n)) if sp["sovereign"] else np.zeros(n)
f_sh = draw("f_share", lambda: rng.triangular(*F_SHARE, n))
w = {"f": (1 - w_v) * f_sh, "s": (1 - w_v) * (1 - f_sh)}
if sp["sovereign"]:
w["v"] = w_v
r = rng.uniform(*STICKY_OVER_FLIGHTY_GROWTH, n)
Gv = rng.triangular(*G_SOV0, n) if sp["sovereign"] else np.zeros(n)
Gf = np.maximum((G_agg - w_v * Gv) / (w["f"] + w["s"] * r), 0.05)
G0 = {"f": Gf, "s": r * Gf, "v": Gv}
ginf = {k: np.log1p(GINF[k][0] + u_ginf * (GINF[k][1] - GINF[k][0])) for k in tiers}
phi = {k: PHI[k][0] + u_phi * (PHI[k][1] - PHI[k][0]) for k in tiers}
else:
w = {"a": np.ones(n)}
G0 = {"a": G_agg}
ginf = {"a": np.log1p(0.08 + u_ginf * 0.17)}
phi = {"a": 0.40 + u_phi * 0.40}
g = {k: np.log1p(G0[k]) for k in tiers}
logR = {k: np.log(w[k]) for k in tiers}
phi_q = {k: phi[k] ** 0.25 for k in tiers}
# ---- stall machinery (macro, shared)
stall_left = np.zeros(n, dtype=int); D = np.zeros(n, dtype=int)
stalled = np.zeros(n, dtype=bool)
g_stall_a = rng.uniform(-0.20, 0.05, n)
kap_s = np.zeros(n); x_s = np.zeros(n); x_f = np.zeros(n); kap_v = np.zeros(n)
ev_left = np.zeros(n, dtype=int); ev_D = np.ones(n, dtype=int); ev_x = np.zeros(n)
lam_f = rng.uniform(*LAM_F, n); lam_v = rng.uniform(*LAM_V, n)
# ---- capital paths
logKs, logKp, pinfo = power_capital(n, rng, on=sp["power"], overrides=o)
h_rev = draw("headroom", lambda: rng.uniform(*REV_HEADROOM, n))
logKs0 = _interp_knots(logKs, NOW); logKp0 = _interp_knots(logKp, NOW)
rho = 0.6
ratio = np.zeros((n, T)); Rtot = np.zeros((n, T)); Rdem = np.zeros((n, T)); plan = np.zeros((n, T)); tierR = {k: np.zeros((n, T)) for k in tiers}
lvl_noise = np.zeros(n)
for j, t in enumerate(times):
start = (stall_left == 0) & (rng.random(n) < 1 - np.exp(-lam * 0.25))
newD = rng.integers(2, 7, n)
stall_left = np.where(start, newD, stall_left); D = np.where(start, newD, D)
g_stall_a = np.where(start, rng.uniform(-0.20, 0.05, n), g_stall_a)
kap_s = np.where(start, o.get("kappa_s", rng.uniform(*KAPPA_S, n)), kap_s)
x_s = np.where(start, o.get("x_s", _tri_ppf(rng.random(n), *X_S)), x_s)
x_f = np.where(start, o.get("x_f", _tri_ppf(rng.random(n), *X_F)), x_f)
kap_v = np.where(start, rng.uniform(*KAPPA_V, n), kap_v)
stalled |= start
in_stall = stall_left > 0
elapsed = D - stall_left # 0-based quarter of the stall
if sp["tiers"] and sp["flighty_hazard"]:
evs = (ev_left == 0) & ~in_stall & (rng.random(n) < 1 - np.exp(-lam_f * 0.25))
newE = rng.integers(2, 5, n)
ev_left = np.where(evs, newE, ev_left); ev_D = np.where(evs, newE, ev_D)
ev_x = np.where(evs, rng.uniform(*XEV_F, n), ev_x)
in_ev = ev_left > 0
z_c = rng.standard_normal(n)
lvl_noise = lvl_noise + sig_l * rng.standard_normal(n)
for k in tiers:
eps = math.sqrt(rho) * z_c + math.sqrt(1 - rho) * rng.standard_normal(n)
if k == "a":
g_eff = np.where(in_stall, g_stall_a, g[k])
elif k == "f":
dur = np.minimum(D, 4)
fall = np.where(elapsed < dur, np.log1p(-np.clip(x_f, 0, 0.9)) * 4 / np.maximum(dur, 1), 0.0)
evf = np.log1p(-np.clip(ev_x, 0, 0.9)) * 4 / np.maximum(ev_D, 1)
g_eff = np.where(in_stall, fall, np.where(in_ev, evf, g[k]))
elif k == "s":
g_eff = np.where(in_stall, (1 - kap_s) * g[k] + np.log1p(-np.clip(x_s, 0, 0.9)) * 4 / np.maximum(D, 1), g[k])
else: # sovereign
g_eff = np.where(in_stall, (1 - kap_v) * g[k], g[k])
logR[k] = logR[k] + g_eff / 4
if k == "v": # programme delays: one-off level loss
dly = rng.random(n) < 1 - np.exp(-lam_v * 0.25)
logR[k] = logR[k] + np.where(dly, np.log1p(-rng.uniform(*XEV_V, n)), 0.0)
restart = in_stall | (in_ev if k == "f" else False)
g[k] = np.where(restart, ginf[k], ginf[k] + phi_q[k] * (g[k] - ginf[k]) + sig_g * 0.5 * eps)
stall_left = np.where(in_stall, stall_left - 1, 0)
ev_left = np.where(in_ev & ~in_stall, ev_left - 1, np.where(in_stall, 0, ev_left))
R_dem = np.exp(lvl_noise) * sum(np.exp(logR[k]) for k in tiers)
for k in tiers:
tierR[k][:, j] = np.exp(logR[k])
# power: capital the plan requires revenue to support, and revenue the powered capital can earn
plan_t = np.exp(_interp_knots(logKs, t) - logKs0)
cap_t = (1 + h_rev) * np.exp((_interp_knots(logKp, t) - logKp0) - (_interp_knots(logKs, t) - logKs0)) * plan_t \
if sp["power"] else np.full(n, np.inf)
R_obs = np.minimum(R_dem, cap_t)
Rdem[:, j] = R_dem; Rtot[:, j] = R_obs; plan[:, j] = plan_t; ratio[:, j] = R_obs / plan_t
shortfall = 1 - ratio
return {"times": times, "ratio": ratio, "shortfall": shortfall, "R": Rtot, "R_dem": Rdem, "plan": plan, "stalled": stalled,
"tier_R": tierR, "tiers": tiers, "power": pinfo, "spec": sp,
"params": {"g0": np.log1p(G_agg), "ginf": ginf[tiers[0]], "phi": 0.4 + 0.4 * u_phi, "lam": lam}}
def g_probabilities(sim, threshold=BUST_SHORTFALL):
"""Cumulative probability that the shortfall has reached `threshold` at any quarter-end by each horizon."""
res = {}
for name, h in HORIZONS.items():
cols = sim["times"] <= h + 1e-9
res[name] = float(np.mean(sim["shortfall"][:, cols].max(axis=1) >= threshold))
return res
def max_shortfall(sim, horizon=2029.0):
cols = sim["times"] <= horizon + 1e-9
return np.clip(sim["shortfall"][:, cols].max(axis=1), 0.0, 0.6)
# Model F mapping (Part II): share of seeds with a neocloud failure at each shortfall (from the 32-seed cliff run; no-shock baseline is not zero).
NC_FAIL_CURVE = ([0.0, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 1.0], [0.0, 0.0, 0.094, 0.188, 0.219, 0.469, 0.656, 0.781, 0.938, 0.938, 0.969, 1.0, 1.0, 1.0])
# [MODEL F v2 output: 32 seeds per shock, all upgrades on; regenerate with python3 -c "import ai_bust_abm as m; m.cliff_curve()"]
# Model F v2: shortfall at which AI-linked stocks fall 40% (inter-quartile range of 12 per-seed crossings: 12.3% to 13.6%; it was
# 7.5-10% in v1 before the Fed loop, sovereign demand and the power ceiling). Chips are more volatile than F's broad AI index, so for
# the SOX we keep the v1 adjustment (-2.5 pts at the low end, -1 pt at the high end): 10% to 12.5%. [MODEL F output]
MARKET_CRASH_SHORTFALL = (0.10, 0.125)
# HOW A MODEL G SHORTFALL MAPS ONTO MODEL F (a key judgment call; see section 1 of the paper).
# Model G has no financial feedback, so the CENTRAL reading treats its shortfall as the underlying demand shock that
# Model F takes as input (Model F's `roi_shock`); Model F then adds the feedback from failures, stocks and the macro
# economy. Each model measures the shortfall against its OWN plan (Part I: revenue keeps pace with the capital stock,
# +66% then +47%; Model F: +35% a year), because each model's firms built for their own plan.
# The ALTERNATIVE reading treats Model G's shortfall as already including the feedback, i.e. as the realised bottom
# of demand vs plan in Model F. Model F's realised bottom for each input shock (median over seeds):
SHOCK_TO_REALISED = ([0.0, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.25, 0.3, 0.4, 0.55],
[0.0, 0.04, 0.08, 0.1235, 0.1615, 0.197, 0.23, 0.26, 0.29, 0.3435, 0.394, 0.488, 0.6215])
# [MODEL F v2 output, median of 12 seeds: regenerate with python3 ai_bust_probability.py --calibrate-mapping]
def _cal_job(args):
import ai_bust_abm as F
s, sd = args
base, scen = F.simulate(replace(F.Params(roi_shock=s, seed=sd)), return_base=True)
sm = F.summarize(scen)
e = F.excess(base, scen)
return {"shock": s, "seed": sd, "realised": 1 - sm["demand_index_min"], "ai_min": sm["ai_index_min"],
"nc_fail": e["excess_defaults"]["neoclouds"] > 0}
def calibrate_mapping(seeds=tuple(range(1, 13)), workers=2):
"""Rebuild the three Model F lookups from Model F v2 (all upgrades on), 12 seeds per shock:
SHOCK_TO_REALISED (realised bottom of demand vs plan), NC_FAIL_CURVE (share of seeds with a neocloud failure) and the
shortfall at which AI-linked stocks fall 40% (per-seed crossing; inter-quartile range)."""
from multiprocessing import Pool
shocks = [0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.15, 0.175, 0.20, 0.25, 0.30, 0.40, 0.55]
jobs = [(s, sd) for s in shocks for sd in seeds]
with Pool(workers) as pool:
rows = pool.map(_cal_job, jobs, chunksize=3)
out = {"shocks": shocks, "n_seeds": len(seeds), "by_shock": []}
cross = []
for sd in seeds:
rs = sorted([r for r in rows if r["seed"] == sd], key=lambda r: r["shock"])
ai = np.minimum.accumulate([r["ai_min"] for r in rs]) # AI index minimum, made monotone in the shock
below = [i for i, v in enumerate(ai) if v <= 0.60]
if below and below[0] > 0:
i = below[0]
x0, x1, y0, y1 = rs[i - 1]["shock"], rs[i]["shock"], ai[i - 1], ai[i]
cross.append(x0 + (0.60 - y0) / (y1 - y0) * (x1 - x0))
elif below:
cross.append(rs[0]["shock"])
for sh in shocks:
rs = [r for r in rows if r["shock"] == sh]
out["by_shock"].append({"shock": sh, "realised_bottom": float(np.median([r["realised"] for r in rs])),
"lo": float(min(r["realised"] for r in rs)), "hi": float(max(r["realised"] for r in rs)),
"p_nc_fail": float(np.mean([r["nc_fail"] for r in rs])),
"ai_index_min_median": float(np.median([r["ai_min"] for r in rs]))})
xs = [b["shock"] for b in out["by_shock"]]
out["shock_to_realised"] = [xs, [b["realised_bottom"] for b in out["by_shock"]]]
pn = np.maximum.accumulate([b["p_nc_fail"] for b in out["by_shock"]])
out["nc_fail_curve"] = [xs[:-1] + [1.0], [float(v) for v in pn[:-1]] + [1.0]]
out["ai_crash_shock_F"] = {"per_seed": [round(float(c), 4) for c in cross], "q25": float(np.percentile(cross, 25)) if cross else None,
"median": float(np.median(cross)) if cross else None, "q75": float(np.percentile(cross, 75)) if cross else None}
return out
def realised_to_shock(x):
"""Model F input shock whose realised demand bottom equals x (inverse of SHOCK_TO_REALISED)."""
return np.interp(x, SHOCK_TO_REALISED[1], SHOCK_TO_REALISED[0])
def shock_to_realised(s):
return np.interp(s, *SHOCK_TO_REALISED)
def fundamentals(rng=RNG):
sim = model_g(rng=rng)
econ = g_probabilities(sim)
thr = rng.uniform(*MARKET_CRASH_SHORTFALL, size=len(sim["shortfall"]))
mkt = {}
for name, h in HORIZONS.items():
cols = sim["times"] <= h + 1e-9
mkt[name] = float(np.mean(sim["shortfall"][:, cols].max(axis=1) >= thr))
# per-draw uncertainty for pooling: bootstrap groups of 400 paths -> spread of the estimate
return {"econ": econ, "market": mkt, "sim": sim}
# =============================================================================
# 2. HISTORY
# =============================================================================
# 2a. Privately financed investment booms: years from the start of the investment acceleration to the year the
# financial bust began (panic, crash or default wave). [HIST; onset years are judgment calls, see paper]
BOOMS = [
# (boom, onset, bust, note)
("US canals", 1834, 1837, "Panic of 1837; state canal-debt defaults 1841-42"),
("UK railway mania", 1844, 1847, "Railway share crash and panic of 1847"),
("US railroads I", 1868, 1873, "Panic of 1873 (Jay Cooke & Co.)"),
("US railroads II", 1879, 1884, "Panic of 1884; receiverships peak 1893"),
("US electric utilities", 1922, 1929, "Utility holding-company collapse 1929-32"),
("US telecom / fiber", 1996, 2001, "Telecom Act 1996; default wave 2001-02 (WorldCom, Global Crossing)"),
("US housing", 2002, 2007, "Subprime crisis 2007-09"),
("US shale oil", 2011, 2015, "Oil-price collapse late 2014; shale defaults 2015-16"),
]
GAPS = np.array([b[2] - b[1] for b in BOOMS], dtype=float)
def history_econ(rng=RNG, n=N, o_lo=2023.5, o_hi=2024.5):
"""P(bust by horizon) = P(eventual bust) x F(horizon - onset).
P(eventual) ~ U(0.60, 0.85): 8 of 8 listed private booms busted (Laplace: 9/10), shaded down for the selection
bias of remembering busts (adding 3-4 unremembered non-busts gives ~0.6). F = lognormal fitted to the 8 gaps.
AI onset ~ U(2023.5, 2024.5): the same definition as the historical rows (start of the investment acceleration):
data-centre capex took off in 2023 and the Big 4 passed $200B in 2024 [REAL]. Sensitivity: dating the onset to
when AI capex passed ~1% of GDP (2025-26) instead lowers this estimate a lot (see paper)."""
mu, sd = np.log(GAPS).mean(), np.log(GAPS).std(ddof=1)
p_ev = rng.uniform(0.60, 0.85, n)
onset = rng.uniform(o_lo, o_hi, n)
# parameter uncertainty in the timing distribution (small sample): bootstrap the fit
mu_b = mu + rng.standard_normal(n) * sd / math.sqrt(len(GAPS))
out = {}
for name, h in HORIZONS.items():
elapsed = np.maximum(h - onset, 1e-3)
F = norm.cdf((np.log(elapsed) - mu_b) / sd)
out[name] = p_ev * F
return out, {"gap_years": GAPS.tolist(), "median_gap": float(np.median(GAPS)), "lognormal_mu": mu, "lognormal_sd": sd}
# 2b. Greenwood, Shleifer & You (2019): P(40% crash within 24 months of a run-up), US 1928-2012 / intl 1987-2012.
GSY = {"us": {0.5: 0.20, 1.0: 0.53, 1.5: 0.80}, "intl": {0.5: 0.36, 1.0: 0.50, 1.5: 0.67}}
# SOX run-up at identification (spring 2026): +181% over 12 months to 23 Jun 2026 [REAL]; two-year raw return
# roughly +160%, S&P 500 ~+40% -> net-of-market ~+120% [REAL-anchored, approximate]
SOX_NET_RUNUP = (1.05, 1.40)
IDENTIFIED = (2026.25, 2026.45) # run-up first identifiable ~Apr-Jun 2026 [ASSUMED from the price path]
BASE_HAZARD = 0.07 # GSY unconditional ~14% per 2 years (US) [HIST]
def gsy_prob(runup, table):
xs = sorted(table)
return float(np.interp(runup, xs, [table[x] for x in xs]))
def crash_timing_share(months):
"""GSY: crashes come after prices peak (~6 months after identification) inside the 24-month window.
Share of eventual crashes that have happened by `months` after identification, uniform over months 6-24."""
return np.clip((months - 6) / 18, 0, 1)
def history_market(rng=RNG, n=N, upper=False):
runup = rng.uniform(*SOX_NET_RUNUP, n)
w = rng.uniform(0, 1, n) # weight on US vs international estimates
if upper: # INDICATORS version: crash-episode characteristics present -> use the 150% row (GSY Table 4)
p_win = w * GSY["us"][1.5] + (1 - w) * GSY["intl"][1.5]
else:
p_win = np.array([w_i * gsy_prob(r, GSY["us"]) + (1 - w_i) * gsy_prob(r, GSY["intl"]) for r, w_i in zip(runup, w)])
ident = rng.uniform(*IDENTIFIED, n)
out = {}
for name, h in HORIZONS.items():
months = (h - ident) * 12
inside = p_win * crash_timing_share(np.minimum(months, 24))
after = np.maximum(h - (ident + 2), 0) # beyond the window: base hazard
out[name] = 1 - (1 - inside) * np.exp(-BASE_HAZARD * after)
return out
# =============================================================================
# 3. MARKET PRICES
# =============================================================================
# LIVE market inputs. The data pipeline (ai_bust_live.py) overwrites these from fresh quotes; the defaults are the
# 5 Oct 2026 values used in the paper. Model G (revenue) is NOT live: it needs lab revenue data that no free feed provides.
LIVE = {"r_free": 0.044, # 3-month bill
"sigma": (0.32, 0.45), # SOX volatility range: NVDA implied .. realised + 6 pts
"sox_drawdown": (0.12, 0.22), # current drawdown from the June-2026 peak
"nc_spread": (0.045, 0.065), # CoreWeave credit spread
"oracle_cds": 0.0227}
R_FREE = 0.044 # 3-month bill ~4.4% (10-yr 5.29% minus the 0.87pt spread, NY Fed Aug-2026) [REAL]
NVDA_IV = 0.32 # NVDA 30-day implied vol 32.2%, 3 Oct 2026; realised 38.9% [REAL]
def barrier_prob(sigma, T, B, r=None):
"""Risk-neutral P(min_{0..T} S <= B * S0) for geometric Brownian motion (reflection principle)."""
r = LIVE["r_free"] if r is None else r
nu = r - 0.5 * sigma ** 2
b = np.log(B)
s = sigma * np.sqrt(T)
return norm.cdf((b - nu * T) / s) + np.exp(2 * nu * b / sigma ** 2) * norm.cdf((b + nu * T) / s)
def market_market(rng=RNG, n=N):
"""Option-implied odds that the SOX falls 40% below its June-2026 peak.
sigma ~ U(0.32, 0.45): NVDA implied 32%, realised 39%; SOX realised volatility at dot-com-era levels [REAL].
Current drawdown d0 ~ U(0.12, 0.22): fell as much as 29% from the late-June peak, rebounded ~11% in September [REAL].
Physical / risk-neutral ratio kappa ~ U(0.6, 0.9): option prices embed a crash-risk premium [ASSUMED]."""
sigma = rng.uniform(*LIVE["sigma"], n)
d0 = rng.uniform(*LIVE["sox_drawdown"], n)
kappa = rng.uniform(0.6, 0.9, n)
B = 0.60 / (1 - d0)
out = {}
for name, h in HORIZONS.items():
out[name] = kappa * barrier_prob(sigma, h - NOW, B)
return out
def credit_implied(rng=RNG, n=N, shortfall_sample=None):
"""Credit-spread-implied odds of a neocloud default, translated into an economic-bust probability.
CoreWeave spread s ~ U(4.5%, 6.5%): term loan repriced at SOFR+550bp (all-in 10.4%, Aug-2026); 5-yr CDS 4.5%
(Jun-2026) after an 8.8% peak (Dec-2025) [REAL]. Hazard = s / (1 - recovery), recovery ~ U(0.3, 0.5).
Physical / risk-neutral kappa ~ U(0.5, 0.8) for high-yield credit [ASSUMED].
Translation: P(bust) = P(neocloud default) x P(shortfall >= 15%) / P(neocloud default | shortfall dist),
using Model F's failure curve and Model G's shortfall SHAPE (this couples the methods; see paper)."""
s = rng.uniform(*LIVE["nc_spread"], n)
rec = rng.uniform(0.3, 0.5, n)
kappa = rng.uniform(0.5, 0.8, n)
lam = s / (1 - rec)
sf = shortfall_sample
p_nc_given = np.interp(sf, *NC_FAIL_CURVE)
ratio = np.mean(sf >= BUST_SHORTFALL) / max(np.mean(p_nc_given), 1e-9)
out_nc, out = {}, {}
for name, h in HORIZONS.items():
pd = kappa * (1 - np.exp(-lam * (h - NOW)))
out_nc[name] = pd
out[name] = np.clip(pd * ratio, 0, 1)
return out, out_nc, float(ratio)
def oracle_implied(cds=None, rec=0.4, kappa=0.65):
"""Cross-check only: Oracle 5-yr CDS 227bp (28 Sep 2026, BBB-) [REAL]."""
cds = LIVE["oracle_cds"] if cds is None else cds
lam = cds / (1 - rec)
return {name: kappa * (1 - math.exp(-lam * (h - NOW))) for name, h in HORIZONS.items()}
# =============================================================================
# 4. WARNING INDICATORS
# =============================================================================
def indicators_econ(rng=RNG, n=N):
"""Greenwood, Hanson, Shleifer & Sorensen (2022): when business credit growth is in its top quintile and equity
prices rose into their top tercile ("R-zone"), the chance of a financial crisis within 3 years is 45% [HIST].
AI-sector analogue: Big-5 bond issuance $121B in 2025 (>4x the 2020-24 average), private credit to AI from
~0 to $200B+, CoreWeave debt $21B -> $35B in six months; SOX +181% in a year [REAL].
The paper's result is country-level, so we use 25-45% over 3 years for the sector analogue [ASSUMED],
with R-zone entry ~ U(2025.5, 2026.0)."""
p3 = rng.uniform(0.25, 0.45, n)
lam = -np.log(1 - p3) / 3
entry = rng.uniform(2025.5, 2026.0, n)
return {name: 1 - np.exp(-lam * (h - entry)) for name, h in HORIZONS.items()}
SCORECARD = [
# indicator, today, pre-bust reference, reading
("Capex vs revenue", "Amazon, Alphabet and Microsoft spend 102% of cloud revenue on capex (2026)",
"Telecom capex ~1.4% of GDP at its 2000 peak; AI capex 1.8% of US GDP in 2026", "Above telecom peak"),
("Debt funding", "33-37% of hyperscaler capex debt-funded; ~$660B of off-balance-sheet leases",
"Telecom 1997-2001 was largely debt-financed", "Rising"),
("Credit growth", "Big-5 issuance $121B in 2025, >4x the 2020-24 average; AI private credit ~$0 -> $200B+",
"R-zone: top-quintile credit growth", "In the zone"),
("Price run-up", "SOX +181% in 12 months to June 2026", "GSY run-up: >=100% in 2 years", "In the zone"),
("Volatility", "SOX: nine 5%+ up-days in 60 sessions (Jun 2026), matched in Jan 2009 and exceeded only in the dot-com era",
"GSY: rising volatility precedes crashes", "Warning"),
("Issuance", "CoreWeave $3.7B convertible (Sep 2026), lab IPO filings", "GSY: issuance rises before crashes", "Warning"),
("Credit spreads", "Oracle 5-yr CDS record 227bp; CoreWeave loan SOFR+550bp with covenants", "Spreads widen before busts", "Widening"),
("Crowding", "53% of fund managers: chips are the most crowded trade (Sep 2026)", "Crowding precedes sharp reversals", "Warning"),
("Rates", "10-yr Treasury 5.29%, highest since 2007", "Tightening often triggers busts", "Warning"),
]
# =============================================================================
# 5. COMBINATION (logarithmic pooling)
# =============================================================================
def pool(method_draws: dict, weights: dict, rng=RNG, n=N, random_weights=True):
"""Weighted average of log-odds across methods (logarithmic pooling). Each draw takes a random value from each
method's own uncertainty AND random method weights (base weight x Gamma(2,1), renormalised), so the range
reflects both input uncertainty and uncertainty about how much to trust each method.
Returns the median and 10th/90th percentiles of the pooled probability."""
names = list(method_draws)
base = np.array([weights[k] for k in names], dtype=float)
if random_weights:
W = base[None, :] * rng.gamma(2.0, 1.0, size=(n, len(names)))
else:
W = np.tile(base, (n, 1))
W = W / W.sum(axis=1, keepdims=True)
L = np.zeros(n)
for j, k in enumerate(names):
d = np.asarray(method_draws[k], dtype=float).ravel()
if d.size == 1:
d = np.full(n, float(d[0]))
L += W[:, j] * logit(rng.choice(d, n))
p = inv_logit(L)
return {"median": float(np.median(p)), "p10": float(np.percentile(p, 10)), "p90": float(np.percentile(p, 90)),
"mean": float(np.mean(p))}
def summarize_method(d):
d = np.asarray(d, dtype=float).ravel()
if d.size == 1:
return {"median": float(d[0]), "p10": float(d[0]), "p90": float(d[0])}
return {"median": float(np.median(d)), "p10": float(np.percentile(d, 10)), "p90": float(np.percentile(d, 90))}
def g_draws(sim, threshold_fn, horizon, groups=400):
"""Turn Model G's path-level outcomes into a distribution of the probability estimate by bootstrapping
groups of paths that share parameter draws (captures prior uncertainty)."""
cols = sim["times"] <= horizon + 1e-9
hit = (sim["shortfall"][:, cols].max(axis=1) >= threshold_fn(len(sim["shortfall"]))).astype(float)
# group by quantiles of the stall hazard and persistence (the two most influential priors)
lam, phi = sim["params"]["lam"], sim["params"]["phi"]
key = np.floor((lam - 0.05) / 0.15 * 10).clip(0, 9) * 10 + np.floor((phi - 0.4) / 0.4 * 10).clip(0, 9)
vals = [hit[key == k].mean() for k in np.unique(key) if (key == k).sum() > 50]
return np.array(vals)
# =============================================================================
# 6. ODDS TRACKER, SENSITIVITY, BACKTEST
# =============================================================================
def tracker_from_g(sim):
"""Condition Model G on what will be observable by mid-2027: annualised growth of end-customer AI
spending over Q4-2026 -> Q2-2027."""
times = sim["times"]
j = int(np.where(np.isclose(times, 2027.5))[0][0]) # end of June 2027
rev_growth = sim["R"][:, j] # R(2027.5)/R(2026.75)
ann = rev_growth ** (1 / 0.75) - 1
cols = times <= HORIZONS["end-2028"] + 1e-9
hit = sim["shortfall"][:, cols].max(axis=1) >= BUST_SHORTFALL
bins = [(-1, 0.30, "below 30%"), (0.30, 0.50, "30-50%"), (0.50, 0.70, "50-70%"), (0.70, 0.90, "70-90%"), (0.90, 10, "above 90%")]
rows = []
for lo, hi, label in bins:
m = (ann >= lo) & (ann < hi)
rows.append({"growth_bin": label, "lo": lo, "hi": hi, "share_of_paths": float(m.mean()),
"p_bust_2028": float(hit[m].mean()) if m.any() else None})
return rows, float(np.median(ann))
def tracker_combined(draws_e, weights_e, draws_m, weights_m, g_rows, rng):
"""How the POOLED end-2028 odds move when one input changes (others held at today's values)."""
out = []
for r in g_rows:
if r["p_bust_2028"] is None:
continue
d = dict(draws_e); d["fundamentals"] = np.array([r["p_bust_2028"]])
out.append({"signal": f"AI spending growth Q4-2026 to mid-2027 (annualised) {r['growth_bin']}",
"definition": "economic bust", "combined": pool(d, weights_e, rng)["median"]})
# credit spreads: CoreWeave spread at its Dec-2025 peak (8.8%) or tighter (4%)
for spread, label in ((0.088, "CoreWeave credit spread back to its Dec-2025 peak (8.8%)"), (0.04, "CoreWeave credit spread tightens to 4%")):
lam = spread / (1 - 0.4)
pd = 0.65 * (1 - math.exp(-lam * (HORIZONS["end-2028"] - NOW)))
d = dict(draws_e); d["market"] = np.array([min(1, pd * CREDIT_RATIO[0])])
out.append({"signal": label, "definition": "economic bust", "combined": pool(d, weights_e, rng)["median"]})
# R-zone exit: AI-related borrowing falls back to normal growth -> indicator estimate drops to the GHSS base
d = dict(draws_e); d["indicators"] = np.array([0.15])
out.append({"signal": "AI borrowing growth falls back to normal (leaves the R-zone)", "definition": "economic bust",
"combined": pool(d, weights_e, rng)["median"]})
# market crash: drawdown deepens to 30% from peak, or SOX makes a new high
for d0, label in ((0.30, "SOX already 30% below its June peak"), (0.0, "SOX back at a new high")):
sigma = 0.38
p = 0.75 * barrier_prob(sigma, HORIZONS["end-2028"] - NOW, 0.60 / (1 - d0))
d = dict(draws_m); d["market"] = np.array([p])
out.append({"signal": label, "definition": "market crash", "combined": pool(d, weights_m, rng)["median"]})
return out
CREDIT_RATIO = [0.77]
def g_tornado(base_overrides=None, n=20_000):
"""Fix one Model G input at its low / high end; report P(economic bust by end-2028)."""
ranges = {"G0": (0.5, 1.5), "G_inf": (0.08, 0.25), "phi": (0.40, 0.80), "sigma_g": (0.08, 0.20), "lam": (0.05, 0.20),
"f_share": (0.30, 0.50), "w_v": (0.04, 0.15), "kappa_s": (0.40, 0.80), "x_f": (0.20, 0.50),
"power_scale": (0.7, 1.3), "headroom": (0.0, 0.30)}
base = g_probabilities(model_g(n=n, rng=np.random.default_rng(7)))["end-2028"]
out = {"base": base}
for k, (lo, hi) in ranges.items():
plo = g_probabilities(model_g(n=n, rng=np.random.default_rng(7), overrides={k: lo}))["end-2028"]
phi_ = g_probabilities(model_g(n=n, rng=np.random.default_rng(7), overrides={k: hi}))["end-2028"]
out[k] = {"low_value": lo, "high_value": hi, "p_at_low": plo, "p_at_high": phi_}
# stress scenarios on the new assumptions
sc = {"sticky tier is not sticky (growth slows 100%, level -25%)": {"kappa_s": 1.0, "x_s": 0.25},
"hard power cap (no revenue headroom)": {"headroom": 0.0},
"power ceiling 30% lower (grid slower than central)": {"power_scale": 0.7},
"flighty share 50%, sovereign 4%": {"f_share": 0.50, "w_v": 0.04}}
out["stress"] = {k: g_probabilities(model_g(n=n, rng=np.random.default_rng(7), overrides=v))["end-2028"] for k, v in sc.items()}
return out
_OFF = {"tiers": False, "flighty_hazard": False, "sovereign": False, "power": False}
def first_bust_quarter(sim, thr=BUST_SHORTFALL):
hit = sim["shortfall"] >= thr
ever = hit.any(axis=1)
first = np.where(ever, sim["times"][hit.argmax(axis=1)], np.nan)
return first, ever
def upgrade_ablation(n=40_000):
"""Effect of each Model G upgrade on P(economic bust) by horizon: a cumulative ladder and leave-one-out from all-on."""
def one(spec):
sim = model_g(n=n, rng=np.random.default_rng(7), spec=spec)
p = g_probabilities(sim)
first, ever = first_bust_quarter(sim)
thr = np.random.default_rng(3).uniform(*MARKET_CRASH_SHORTFALL, size=n)
mc = float(np.mean(sim["shortfall"][:, sim["times"] <= 2029.0 + 1e-9].max(axis=1) >= thr))
return {**{k: round(v, 4) for k, v in p.items()}, "market_crash_end-2028": round(mc, 4),
"median_shortfall_end-2028": float(np.median(max_shortfall(sim, 2029.0))),
"median_bust_date_given_bust_by_2030": float(np.nanmedian(first[ever])) if ever.any() else None}
ladder = [("v1: one demand tier, plan capex", _OFF),
("+ flighty / sticky tiers", {**_OFF, "tiers": True}),
("+ sovereign tier", {**_OFF, "tiers": True, "sovereign": True}),
("+ flighty-only disillusionment events", {**_OFF, "tiers": True, "sovereign": True, "flighty_hazard": True}),
("+ power ceiling (all on)", {**_OFF, "tiers": True, "sovereign": True, "flighty_hazard": True, "power": True})]
out = {"ladder": [{"step": a, **one(b)} for a, b in ladder]}
allon = {"tiers": True, "sovereign": True, "flighty_hazard": True, "power": True}
out["leave_one_out"] = {
"all on": one(allon),
"without tiers (single tier, keeps power)": one({**_OFF, "power": True}),
"without sovereign": one({**allon, "sovereign": False}),
"without flighty events": one({**allon, "flighty_hazard": False}),
"without power ceiling": one({**allon, "power": False})}
out["power_only"] = one({**_OFF, "power": True})
return out
def tier_diagnostics(sim):
times = sim["times"]; out = {}
idx = {str(int(y)): int(np.where(np.isclose(times, y))[0][0]) for y in (2027.0, 2028.0, 2029.0, 2030.0)}
tot = sum(sim["tier_R"][k] for k in sim["tiers"])
out["median_revenue_index_by_tier"] = {k: {y: float(np.median(sim["tier_R"][k][:, j])) for y, j in idx.items()} for k in sim["tiers"]}
out["median_share_by_tier"] = {k: {y: float(np.median(sim["tier_R"][k][:, j] / tot[:, j])) for y, j in idx.items()} for k in sim["tiers"]}
out["median_total_demand_index"] = {y: float(np.median(tot[:, j])) for y, j in idx.items()}
pw = sim["power"]
if pw.get("ceiling") is not None:
out["power"] = {"ceiling_T_median": [float(x) for x in np.median(pw["ceiling"], axis=0)],
"plan_capex_T": [float(x) for x in pw["plan"]],
"p_ceiling_binds": [float(x) for x in (pw["excess_share"] > 0).mean(axis=0)],
"median_excess_share": [float(x) for x in np.median(pw["excess_share"], axis=0)],
"p_revenue_capped_by_power": {y: float(np.mean(sim["R_dem"][:, j] > sim["R"][:, j] + 1e-12)) for y, j in idx.items()}}
kp = np.median(sim["plan"], axis=0)
out["power"]["median_plan_index"] = {y: float(kp[j]) for y, j in idx.items()}
out["power"]["v1_plan_index"] = {y: float(plan_index(float(y)) / plan_index(NOW)) for y in idx}
return out
def softer_capex_plan():
"""Sensitivity: if 2028-29 capex is cut 20% vs today's plans, the plan path is lower."""
global CAPITAL
saved = dict(CAPITAL)
CAPITAL = {2026.0: 1.05, 2027.0: 2.05, 2028.0: 3.40, 2029.0: 3.40 + 0.8 * 1.60, 2030.0: 3.40 + 0.8 * 1.60 + 0.8 * 1.80}
p = g_probabilities(model_g(n=20_000, rng=np.random.default_rng(7)))
CAPITAL = saved
return p
BACKTEST = [
# case, as-of, history rule, GSY run-up rule, R-zone rule, fundamentals rule, outcome
{"case": "Telecom / dot-com", "as_of": "Dec 1999",
"history": "Telecom investment boom began ~1996 -> 3 yrs in; earlier booms busted after a median 5 yrs",
"gsy": "Nasdaq ~+160% in 2 yrs, ~+105% net of S&P -> 53-80%",
"rzone": "Corporate credit growth and equity prices both high -> yes",
"fundamentals": "Capacity built for traffic 'doubling every 100 days' (~4x/yr) vs ~2x/yr actual -> large shortfall",
"signal": "High on all four", "outcome": "Market crash (Nasdaq -78%) and economic bust (2001-02)"},
{"case": "Cloud buildout", "as_of": "Oct 2018",
"history": "Hyperscaler capex ~0.4% of GDP -> below the boom threshold",
"gsy": "Big-tech 2-yr run-up well under +100% net of market -> base rate ~14%",
"rzone": "Mostly cash-funded, credit growth normal -> no",
"fundamentals": "Cloud revenue growing ~40%, capacity tracked demand",
"signal": "Low on all four", "outcome": "No bust; SOX fell ~25% in Q4 2018, short of 40%"},
{"case": "Cloud software", "as_of": "Dec 2020",
"history": "Not an investment boom",
"gsy": "Cloud-software stocks well over +100% in 2 yrs -> 53-80%",
"rzone": "Equity high, credit growth not top quintile -> no",
"fundamentals": "Revenue grew; no capacity overhang",
"signal": "High for a market crash, low for an economic bust", "outcome": "Market crash (cloud indices -60%+ in 2022), no economic bust"},
]
# =============================================================================
# 7. EXPECTED DAMAGE (Part III): rerun Models D and F on the estimated shortfall distribution
# =============================================================================
def reweighted_sampler(sf, p_bust, rng):
"""Sample shortfalls from Model G's max-shortfall distribution, reweighted so that P(shortfall >= 15%)
equals the pooled economic-bust probability."""
hi = sf[sf >= BUST_SHORTFALL]
lo = sf[sf < BUST_SHORTFALL]
def sample(r=None):
g = r if r is not None else rng
return float(g.choice(hi)) if g.random() < p_bust else float(g.choice(lo))
return sample
def expected_damage(sf, p_bust, n_d=10_000, n_f=400, run_f=True, workers=2):
import ai_bust_models as D
out = {}
rng = np.random.default_rng(11)
samp = reweighted_sampler(sf, p_bust, rng)
rows = D.monte_carlo(n_d, shock_sampler=lambda: samp())
cl = np.array([r["credit_losses"] for r in rows])
out["model_d"] = {
"n": n_d,
"p_bust_outcome": float(np.mean([r["outcome"] in ("Bust", "Systemic crisis") for r in rows])),
"p_sp500_30": float(np.mean([r["sp500_drawdown"] >= 0.30 for r in rows])),
"expected_credit_losses": float(cl.mean()),
"p_credit_gt_100": float(np.mean(cl > 100)), "p_credit_gt_250": float(np.mean(cl > 250)),
"expected_gdp_hit": float(np.mean([r["gdp_hit_pct"] for r in rows])),
"p_neocloud_default": float(np.mean([r["neocloud_default"] for r in rows])),
"p_lab_default": float(np.mean([r["lab_default"] for r in rows])),
"p_private_credit_default": float(np.mean([r["pc_default"] for r in rows])),
"p_bank_default": float(np.mean(["Banks" in r["defaults"] for r in rows])),
"expected_bank_loss": float(np.mean([r["bank_loss"] for r in rows])),
"expected_pension_insurer_loss": float(np.mean([r["pension_loss"] for r in rows])),
"exceedance": [{"loss": x, "p": float(np.mean(cl > x))} for x in (0, 10, 25, 50, 100, 150, 200, 250, 300, 350)],
}
# Model D with stochastic legal friction (own sampler copy, same shortfall distribution)
samp2 = reweighted_sampler(sf, p_bust, np.random.default_rng(11))
D.RNG = np.random.default_rng(20260930)
rows_l = D.monte_carlo(n_d, shock_sampler=lambda: samp2(), contagion_kwargs={"legal_on": True})
cl2 = np.array([r["credit_losses"] for r in rows_l]); fz = np.array([r["frozen_peak"] for r in rows_l])
out["model_d_legal"] = {"n": n_d, "expected_credit_losses_booked": float(cl2.mean()),
"expected_frozen_peak": float(fz.mean()), "p90_frozen_peak": float(np.percentile(fz, 90)),
"p_frozen_gt_100": float(np.mean(fz > 100)),
"p_bust_outcome": float(np.mean([r["outcome"] in ("Bust", "Systemic crisis") for r in rows_l]))}
if run_f:
import ai_bust_abm as F
res = F.monte_carlo(n_f, seed=5, shock_sampler=lambda r: samp(r), workers=workers)
out["model_f"] = summarize_model_f(res["rows"])
res1 = F.monte_carlo(n_f, seed=5, shock_sampler=lambda r: samp(r), params_override=F._V1_OFF, workers=workers)
out["model_f_v1"] = summarize_model_f(res1["rows"])
out["model_f_by_shock"] = res["by_shock"]
return out
def summarize_model_f(rf):
"""Part III summary of a Model F Monte Carlo (outcomes net of each draw's no-shock twin; runs Jan 2027-end 2029)."""
xs = np.array([max(0.0, r["excess_credit_losses"]) for r in rf])
return {
"n": len(rf),
"p_neocloud_default": float(np.mean([r["excess_defaults"]["neoclouds"] > 0 for r in rf])),
"p_lab_default": float(np.mean([r["excess_defaults"]["labs"] > 0 for r in rf])),
"p_fund_failure": float(np.mean([r["excess_defaults"]["pc_funds"] > 0 for r in rf])),
"p_bank_failure": float(np.mean([r["excess_defaults"]["banks"] > 0 for r in rf])),
"expected_credit_losses": float(xs.mean()),
"p_credit_gt_50": float(np.mean(xs > 50)), "p_credit_gt_100": float(np.mean(xs > 100)),
"p_sp500_30": float(np.mean([r["sp500_drawdown"] >= 0.30 for r in rf])),
"p_unemp_6": float(np.mean([r["unemp_peak"] >= 6.0 for r in rf])),
"expected_unemp_peak": float(np.mean([r["unemp_peak"] for r in rf])),
"expected_sovereign_spend": float(np.mean([max(0.0, r["excess_sovereign_spend"]) for r in rf])),
"expected_frozen_claims_peak": float(np.mean([r["frozen_claims_peak"] for r in rf])),
"p90_frozen_claims_peak": float(np.percentile([r["frozen_claims_peak"] for r in rf], 90)),
"mean_cfo_actions": float(np.mean([r["cfo"]["n_actions"] for r in rf])),
"expected_debt_exchange_haircut": float(np.mean([r["cfo"]["haircut_total"] for r in rf])),
"expected_equity_injected": float(np.mean([r["cfo"]["injected"] for r in rf])),
"expected_interest_credited_to_lenders": float(np.mean([r["sfc"]["interest_credited_to_lenders"] for r in rf if r.get("sfc")] or [0.0])),
"expected_household_credit_losses": float(np.mean([r["sfc"]["household_credit_losses"] for r in rf if r.get("sfc")] or [0.0])),
"sfc_max_claim_imbalance": float(max([r["sfc"]["max_claim_imbalance"] for r in rf if r.get("sfc")] or [0.0])),
"sfc_max_networth_error": float(max([r["sfc"]["max_networth_reconciliation_error"] for r in rf if r.get("sfc")] or [0.0])),
"network_top1_fund_share": float(np.mean([r["network"]["top1_fund_share"] for r in rf])),
"exceedance": [{"loss": x, "p": float(np.mean(xs > x))} for x in (0, 10, 25, 50, 75, 100, 125, 150)],
}
# =============================================================================
# 8. EXTRA SENSITIVITIES (own random streams, so the headline results above are unchanged)
# =============================================================================
def mapping_sensitivity(sf28, p_bust_28, draws_e, weights_e, draws_m, weights_m, sim, run_f=True, n_f=400):
"""ALTERNATIVE reading: Model G's shortfall already includes the financial feedback (= Model F's realised
demand bottom), so the matching Model F input shock is smaller. Recomputes everything that uses Model F:
the 15% line, the credit-to-bust ratio, the market-crash line for Model G, the pooled odds and Model F's damage."""
rng = np.random.default_rng(99)
out = {}
nc_real = lambda x: np.interp(realised_to_shock(x), *NC_FAIL_CURVE)
out["p_neocloud_failure_at_15pct"] = {"central": float(np.interp(0.15, *NC_FAIL_CURVE)),
"alternative": float(nc_real(0.15))}
# credit -> bust ratio
ratio_alt = float(np.mean(sf28 >= BUST_SHORTFALL) / np.mean(nc_real(sf28)))
out["credit_to_bust_ratio"] = {"central": CREDIT_RATIO[0], "alternative": ratio_alt}
# market-crash line for Model G, in realised terms
lo, hi = shock_to_realised(MARKET_CRASH_SHORTFALL[0]), shock_to_realised(MARKET_CRASH_SHORTFALL[1])
out["market_crash_line"] = {"central": list(MARKET_CRASH_SHORTFALL), "alternative": [float(lo), float(hi)]}
gm_alt = g_draws(sim, lambda n: rng.uniform(lo, hi, size=n), HORIZONS["end-2028"])
de = dict(draws_e); de["market"] = np.clip(np.asarray(draws_e["market"]) * ratio_alt / CREDIT_RATIO[0], 0, 1)
dm = dict(draws_m); dm["fundamentals"] = gm_alt
out["end_2028"] = {
"credit_implied_bust": {"central": float(np.median(draws_e["market"])), "alternative": float(np.median(de["market"]))},
"fundamentals_market_crash": {"central": float(np.median(draws_m["fundamentals"])), "alternative": float(np.median(gm_alt))},
"economic_bust_combined": pool(de, weights_e, rng),
"market_crash_combined": pool(dm, weights_m, rng),
}
if run_f:
import ai_bust_abm as F
samp = reweighted_sampler(sf28, p_bust_28, np.random.default_rng(11))
res = F.monte_carlo(n_f, seed=5, shock_sampler=lambda r: float(realised_to_shock(samp(r))), workers=2)
out["model_f_alternative"] = summarize_model_f(res["rows"])
return out
def stall_decomposition():
"""How much of Model G's bust risk comes from demand stalls (after a stall, growth restarts at the long-run rate)."""
with_s = g_probabilities(model_g(n=40_000, rng=np.random.default_rng(7)))
no_s = g_probabilities(model_g(n=40_000, rng=np.random.default_rng(7), overrides={"lam": 1e-9}))
return {"with_stalls": with_s, "no_stalls": no_s}
def mid_2027_odds(sim, weights_e):
"""Combined economic-bust odds by mid-2027, to compare with Polymarket's 'AI industry downturn by 30 Jun 2027'."""
global HORIZONS
saved = HORIZONS
HORIZONS = {"mid-2027": 2027.5}
try:
rng = np.random.default_rng(123)
sf = max_shortfall(sim, 2027.5)
ge = g_draws(sim, lambda n: np.full(n, BUST_SHORTFALL), 2027.5)
he, _ = history_econ(rng)
ce, _, _ = credit_implied(rng, shortfall_sample=max_shortfall(sim, 2029.0))
ie = indicators_econ(rng)
d = {"fundamentals": ge, "history": he["mid-2027"], "market": ce["mid-2027"], "indicators": ie["mid-2027"]}
return {"methods": {k: float(np.median(v)) for k, v in d.items()},
"fundamentals_point": float(np.mean(sf >= BUST_SHORTFALL)),
"combined": pool(d, weights_e, rng)}
finally:
HORIZONS = saved
def run_live(n=10_000, seed=20261006, n_d=4_000, n_f=40, f_overrides=None, with_f=True):
"""Lean refresh used by the live pipeline: re-pool the four methods with the current LIVE market inputs.
Model G (revenue) is held at its paper calibration; Model F is re-run on a small Monte Carlo with `f_overrides`."""
global N, RNG
saved = (N, RNG)
N, RNG = n, np.random.default_rng(seed)
try:
rng = RNG
fund = fundamentals(rng)
sim = fund["sim"]
sf28 = max_shortfall(sim, HORIZONS["end-2028"])
hist_e, _ = history_econ(rng, n=n)
hist_m = history_market(rng, n=n)
ind_m = history_market(rng, n=n, upper=True)
mkt_m = market_market(rng, n=n)
cred_e, cred_nc, ratio = credit_implied(rng, n=n, shortfall_sample=sf28)
CREDIT_RATIO[0] = ratio
ind_e = indicators_econ(rng, n=n)
thr_m = lambda m: rng.uniform(*MARKET_CRASH_SHORTFALL, size=m)
thr_e = lambda m: np.full(m, BUST_SHORTFALL)
we = {"fundamentals": 1, "history": 1, "market": 1, "indicators": 1}
wm = {"fundamentals": 1, "history": 0.5, "market": 1, "indicators": 0.5}
out = {"economic_bust": {}, "market_crash": {}, "methods_end2028": {}}
for name, h in HORIZONS.items():
ge, gm = g_draws(sim, thr_e, h), g_draws(sim, thr_m, h)
de = {"fundamentals": ge, "history": hist_e[name], "market": cred_e[name], "indicators": ind_e[name]}
dm = {"fundamentals": gm, "history": hist_m[name], "market": mkt_m[name], "indicators": ind_m[name]}
out["economic_bust"][name] = pool(de, we, rng, n=n)
out["market_crash"][name] = pool(dm, wm, rng, n=n)
if name == "end-2028":
out["methods_end2028"] = {"bust": {k: float(np.median(v)) for k, v in de.items()},
"crash": {k: float(np.median(v)) for k, v in dm.items()}}
p_bust = out["economic_bust"][name]["median"]
out["expected_damage"] = expected_damage(sf28, p_bust, n_d=n_d, n_f=n_f, run_f=False)
if with_f:
import ai_bust_abm as F
samp = reweighted_sampler(sf28, p_bust, np.random.default_rng(11))
res = F.monte_carlo(n_f, seed=5, shock_sampler=lambda r: float(samp(r)), params_override=f_overrides)
out["expected_damage"]["model_f"] = summarize_model_f(res["rows"])
out["credit_to_bust_ratio"] = float(ratio)
out["live_inputs"] = {k: (list(v) if isinstance(v, tuple) else v) for k, v in LIVE.items()}
return out
finally:
N, RNG = saved
# =============================================================================
# MAIN
# =============================================================================
def run(fast=False):
rng = RNG
res = {"as_of": "2026-10-05"}
fund = fundamentals(rng)
sim = fund["sim"]
sf28 = max_shortfall(sim, HORIZONS["end-2028"])
hist_e, hist_info = history_econ(rng)
hist_m = history_market(rng)
ind_m = history_market(rng, upper=True)
mkt_m = market_market(rng)
cred_e, cred_nc, ratio = credit_implied(rng, shortfall_sample=sf28)
CREDIT_RATIO[0] = ratio
ind_e = indicators_econ(rng)
thr_m = lambda n: rng.uniform(*MARKET_CRASH_SHORTFALL, size=n)
thr_e = lambda n: np.full(n, BUST_SHORTFALL)
weights_e = {"fundamentals": 1, "history": 1, "market": 1, "indicators": 1}
weights_m = {"fundamentals": 1, "history": 0.5, "market": 1, "indicators": 0.5} # history+indicators share GSY
res["economic_bust"], res["market_crash"] = {}, {}
for name, h in HORIZONS.items():
ge = g_draws(sim, thr_e, h)
gm = g_draws(sim, thr_m, h)
draws_e = {"fundamentals": ge, "history": hist_e[name], "market": cred_e[name], "indicators": ind_e[name]}
draws_m = {"fundamentals": gm, "history": hist_m[name], "market": mkt_m[name], "indicators": ind_m[name]}
res["economic_bust"][name] = {
"methods": {"fundamentals": {"median": fund["econ"][name], **{k: v for k, v in summarize_method(ge).items() if k != "median"}},
**{k: summarize_method(v) for k, v in draws_e.items() if k != "fundamentals"}},
"combined": pool(draws_e, weights_e, rng)}
res["market_crash"][name] = {
"methods": {"fundamentals": {"median": fund["market"][name], **{k: v for k, v in summarize_method(gm).items() if k != "median"}},
**{k: summarize_method(v) for k, v in draws_m.items() if k != "fundamentals"}},
"combined": pool(draws_m, weights_m, rng)}
# leave-one-out
for key, draws, wts in (("economic_bust", draws_e, weights_e), ("market_crash", draws_m, weights_m)):
loo = {}
for drop in draws:
d2 = {k: v for k, v in draws.items() if k != drop}
w2 = {k: wts[k] for k in d2}
loo[drop] = pool(d2, w2, rng)["median"]
res[key][name]["leave_one_out"] = loo
res["cross_checks"] = {
"neocloud_default_credit_implied": {k: summarize_method(v) for k, v in cred_nc.items()},
"credit_to_bust_ratio": ratio,
"oracle_default_cds_implied": oracle_implied(),
"polymarket": {"industry downturn by 31 Dec 2026": 0.071, "by 30 Jun 2027 (thin market)": 0.18},
"ny_fed_recession_12m": 0.139,
}
res["history_info"] = hist_info
# Model G detail
times = sim["times"]
res["model_g"] = {
"plan_growth": {str(int(y) - 1): float(np.median(sim["plan"][:, int(np.where(np.isclose(times, y))[0][0])]
/ sim["plan"][:, int(np.where(np.isclose(times, y - 1))[0][0])]) - 1)
for y in (2028.0, 2029.0, 2030.0)},
"fan": [{"t": float(t), "p10": float(np.percentile(sim["ratio"][:, j], 10)),
"p25": float(np.percentile(sim["ratio"][:, j], 25)),
"p50": float(np.median(sim["ratio"][:, j])),
"p75": float(np.percentile(sim["ratio"][:, j], 75)),
"p90": float(np.percentile(sim["ratio"][:, j], 90))} for j, t in enumerate(times)],
"p_stall_by_end_2029": float(np.mean(sim["stalled"])),
"shortfall_2028_pct": {q: float(np.percentile(sf28, q)) for q in (10, 25, 50, 75, 90)},
"econ": fund["econ"], "market": fund["market"],
"median_revenue_growth": {
str(int(y) - 1): float(np.median(
sim["R"][:, int(np.where(np.isclose(times, y))[0][0])]
/ sim["R"][:, int(np.where(np.isclose(times, y - 1))[0][0])]) - 1)
for y in (2028.0, 2029.0, 2030.0)},
}
res["tracker_g"], res["median_h1_2027_growth"] = tracker_from_g(sim)
h = "end-2028"
ge28 = g_draws(sim, thr_e, HORIZONS[h]); gm28 = g_draws(sim, thr_m, HORIZONS[h])
de = {"fundamentals": ge28, "history": hist_e[h], "market": cred_e[h], "indicators": ind_e[h]}
dm = {"fundamentals": gm28, "history": hist_m[h], "market": mkt_m[h], "indicators": ind_m[h]}
res["tracker_combined"] = tracker_combined(de, weights_e, dm, weights_m, res["tracker_g"], rng)
late, _ = history_econ(rng, o_lo=2025.25, o_hi=2026.25)
res["history_onset_sensitivity"] = {"onset 2025-26 (AI capex >1% of GDP)": {k: float(np.median(v)) for k, v in late.items()}}
res["combined_fixed_weights"] = {
"economic_bust_end-2028": pool(de, weights_e, rng, random_weights=False),
"market_crash_end-2028": pool(dm, weights_m, rng, random_weights=False)}
res["tornado_g"] = g_tornado()
res["tier_diagnostics"] = tier_diagnostics(sim)
res["g_upgrade_ablation"] = upgrade_ablation()
res["softer_capex_plan"] = softer_capex_plan()
res["backtest"] = BACKTEST
res["scorecard"] = SCORECARD
res["booms"] = [{"boom": b[0], "onset": b[1], "bust": b[2], "gap": b[2] - b[1], "note": b[3]} for b in BOOMS]
# Part III
p_bust_28 = res["economic_bust"]["end-2028"]["combined"]["median"]
res["expected_damage"] = expected_damage(sf28, p_bust_28, run_f=not fast)
res["expected_damage"]["p_bust_used"] = p_bust_28
# Odds of deeper shortfalls by end-2028, on the same reweighted distribution (no random draws consumed)
hi = sf28[sf28 >= BUST_SHORTFALL]
res["expected_damage"]["p_shortfall_at_least"] = {
f"{int(t * 100)}%": float(p_bust_28 * np.mean(hi >= t)) for t in (0.15, 0.20, 0.25, 0.30, 0.40)}
# Extra sensitivities (own random streams)
res["mapping_sensitivity"] = mapping_sensitivity(sf28, p_bust_28, de, weights_e, dm, weights_m, sim, run_f=not fast)
res["stall_decomposition"] = stall_decomposition()
res["mid_2027"] = mid_2027_odds(sim, weights_e)
return res
if __name__ == "__main__":
ap = argparse.ArgumentParser(description="Part I odds and Part III expected damage for the AI-bust paper.")
ap.add_argument("--fast", action="store_true", help="skip the Model F rerun")
ap.add_argument("--out", default=None,
help="output file (default probability_results.json, or probability_results_fast.json with --fast)")
ap.add_argument("--calibrate-mapping", action="store_true",
help="only rebuild the Model G -> Model F shortfall mapping and write abm_shock_to_realised.json")
a = ap.parse_args()
if a.calibrate_mapping:
with open("abm_shock_to_realised.json", "w") as f:
json.dump(calibrate_mapping(), f, indent=1)
raise SystemExit(0)
out = a.out or ("probability_results_fast.json" if a.fast else "probability_results.json")
r = run(fast=a.fast)
with open(out, "w") as f:
json.dump(r, f, indent=1, default=float)
for key in ("economic_bust", "market_crash"):
print(key)
for h, v in r[key].items():
print(" ", h, {m: round(x["median"], 3) for m, x in v["methods"].items()},
"-> combined", round(v["combined"]["median"], 3), f"({v['combined']['p10']:.2f}-{v['combined']['p90']:.2f})")