A corrugator engineer values International Paper

"""My bet: the value distribution mapped onto the instruments I could actually buy.

The Monte Carlo gives a distribution of what the shares are WORTH. A position pays on
what they TRADE at, and the market can disagree with value for a long time. So the one
extra assumption here is convergence: over the horizon the price moves a fraction c of
the way from $32 toward the simulated value (c = 1 is full convergence, 0 is none).

Option premiums are ILLUSTRATIVE placeholders — check the live chain before trading.
Run: python3 bet.py  (writes bet.json for the post's chart)
"""
import io, contextlib, json
import numpy as np

# re-run the simulation to get the value draws (same seed as montecarlo.py)
src = open("montecarlo.py").read().replace('json.dump(out, open("mc.json", "w"), indent=1)', '')
ns = {}
with contextlib.redirect_stdout(io.StringIO()):
    exec(src, ns)
value, down, price = ns["value"], ns["down"], ns["price"]
rng = np.random.default_rng(7)

# ---- horizon price: partial convergence toward value, plus noise the model can't see
c = 0.6                                              # 60% of the gap closes over ~12–15 months
price_T = price + c * (value - price)
MKT_NOISE = 0.12                                     # lognormal market noise the model can't see
price_T = price_T * np.exp(rng.normal(0, MKT_NOISE, len(value)))
price_T = np.maximum(price_T, 3)

# ---- dividend over the horizon: $1.85 held in the base regime, halved in a down-cycle (cut)
div = np.where(down, 0.925, 1.85)

# ---- instruments (per 100 shares / one contract). Premiums are placeholders.
CALL_K, CALL_PREM = 35.0, 3.50          # Jan-2028 $35 call, illustrative
PUT_K, PUT_PREM = 28.0, 1.75            # 6-month $28 cash-secured put, illustrative
SPR_K1, SPR_K2, SPR_NET = 35.0, 45.0, 2.25   # 35/45 call spread, illustrative net debit

instruments = {
    "100 shares at $32": dict(pnl=100 * (price_T - price) + 100 * div, capital=100 * price),
    "1 call, Jan-2028 $35 (~$3.50)": dict(pnl=100 * (np.maximum(0, price_T - CALL_K) - CALL_PREM), capital=100 * CALL_PREM),
    "1 cash-secured put, $28 6-mo (~$1.75)": dict(pnl=100 * (PUT_PREM - np.maximum(0, PUT_K - price_T)), capital=100 * PUT_K),
    "1 call spread 35/45 (~$2.25)": dict(pnl=100 * (np.minimum(SPR_K2 - SPR_K1, np.maximum(0, price_T - SPR_K1)) - SPR_NET), capital=100 * SPR_NET),
}
rows = []
for name, d in instruments.items():
    p = d["pnl"]
    rows.append(dict(name=name, capital=float(d["capital"]), mean=float(p.mean()), p_profit=float((p > 0).mean()),
                     p10=float(np.percentile(p, 10)), p50=float(np.percentile(p, 50)), p90=float(np.percentile(p, 90)),
                     worst=float(p.min()), roc=float(p.mean() / d["capital"]),
                     mean_down=float(p[down].mean()), mean_base=float(p[~down].mean())))
    print(f"{name:40s} capital ${d['capital']:>6,.0f}  E[P&L] ${p.mean():>7,.0f}  P(profit) {(p>0).mean():5.1%}  10th ${np.percentile(p,10):>7,.0f}  90th ${np.percentile(p,90):>7,.0f}  worst ${p.min():>7,.0f}  on capital {p.mean()/d['capital']:+.1%}")

# payoff curves for the chart (per 100 shares, before dividends)
grid = np.arange(10, 70.5, 0.5)
curves = {
    "shares": (100 * (grid - price)).tolist(),
    "call": (100 * (np.maximum(0, grid - CALL_K) - CALL_PREM)).tolist(),
    "put": (100 * (PUT_PREM - np.maximum(0, PUT_K - grid))).tolist(),
    "spread": (100 * (np.minimum(SPR_K2 - SPR_K1, np.maximum(0, grid - SPR_K1)) - SPR_NET)).tolist(),
}
hist, edges = np.histogram(price_T, bins=np.arange(0, 100, 2.5))
json.dump(dict(convergence=c, noise=MKT_NOISE, spread_wipeout=float((price_T <= SPR_K1).mean()), price=price, grid=grid.tolist(), curves=curves, rows=rows,
               priceT=dict(edges=edges.tolist(), counts=hist.tolist(), p10=float(np.percentile(price_T, 10)), p50=float(np.percentile(price_T, 50)), p90=float(np.percentile(price_T, 90))),
               premiums=dict(call=CALL_PREM, put=PUT_PREM, spread=SPR_NET)), open("bet.json", "w"))
print(f"\nhorizon price (c={c}): 10th ${np.percentile(price_T,10):.0f}  median ${np.percentile(price_T,50):.0f}  90th ${np.percentile(price_T,90):.0f}")

PASS  market cap = price × shares: $16,944M vs $16,944M
PASS  net debt = gross debt − cash: $8,491M
PASS  10-K pension lines tie out: surplus plans − deficit plans = the filing's net, and the U.S. plan's table = its stated surplus: $486M − $342M = $144M (10-K: $144M); U.S. qualified $8,482M − $8,116M = $366M
PASS  debt-like items rebuild from the 10-K lines and come off the equity, not the WACC: $957M after tax = $957M from the lines; equity = EV − $8,491M − $957M
PASS  2026E EBITDA = guidance midpoint: $3,300M vs $3,300M guide midpoint
PASS  2026E levered FCF inside the company's Jan guide ($300–500M): $303M (guide given before the April cut; the low end is the honest reading)
PASS  dividend cost ≈ $980M: $980M = $1.85 × 529.5M
PASS  2026E levered FCF does NOT cover the dividend: 0.31× coverage
PASS  Gordon TV implies an exit multiple between SW and PKG: 8.2x (SW 7.2x, PKG 12.2x)
PASS  Exit-multiple TV implies a sane perpetual growth (0–4%): 1.89%
PASS  WACC > terminal growth: 7.36% > 2.0%
PASS  terminal value share of EV disclosed and in the mature-company band (60–90%): 82%
PASS  DCF sheet ≈ Scenarios mid-cycle column (same chain, two layouts): $38.44 vs $38.69
PASS  Sensitivity grid (17% × 8.5x) ≈ DCF within the interpolation error: $40.11 vs $38.44
PASS  Bear < Mid < Bull: $11 < $39 < $69
PASS  25% of $170/ton moves 2027E EBITDA by 42.5 × $9M: $382.5M vs $382.5M (CFO: $10/ton ≈ $90M)
PASS  Python chain reproduces Scenarios!C21 (2027E EBITDA): $4,136M vs $4,136M
PASS  P(value > price) falls monotonically as P(down-cycle) rises: 79% → 71% → 62% → 53% → 45% → 37% → 28%
PASS  breakeven P(down-cycle) ≈ 50% (market price = expected value): mean at p=0.5: $30.9
PASS  MC mean sits between the Bear and Bull point estimates: $35.9
PASS  regime split is real: P(>price | down) ≪ P(>price | base): 3% vs 87%

0 failure(s)

"""Validation: does the model tie out?  Run: python3 validate.py

Each check is a thing a reviewer would ask. A check prints PASS/FAIL with the numbers,
and the script exits non-zero if anything fails, so it can gate a republish.
"""
import json, sys, subprocess
from openpyxl import load_workbook

wb = load_workbook("IP_valuation_model.xlsx", data_only=True)
inp = {r[0].value: r[1].value for r in wb["Inputs"].iter_rows(min_row=5) if r[0].value and r[1].value is not None}
M, D, S, V = wb["Model"], wb["DCF"], wb["Scenarios"], wb["Sensitivity"]
DCF = {}  # DCF sheet by label, so a new row can't silently shift what a check reads
for row in D.iter_rows(min_row=1):
    if isinstance(row[0].value, str) and row[1].value is not None:
        DCF.setdefault(row[0].value.strip(), row[1].value)
mc = json.load(open("mc.json")); pdown = json.load(open("mc_pdown.json"))
fails = 0
def check(name, ok, detail):
    global fails
    print(f"{'PASS' if ok else 'FAIL'}  {name}: {detail}")
    fails += 0 if ok else 1

# --- identities and tie-outs to what the company said
price, shares = inp["Share price ($)"], inp["Diluted shares (millions)"]
mcap, netdebt, debtlike = DCF["Market capitalization"], DCF["Net debt (gross debt − cash)"], DCF["Debt-like items from the 10-K, after tax (Inputs)"]
check("market cap = price × shares", abs(mcap - price * shares) < 1, f"${mcap:,.0f}M vs ${price*shares:,.0f}M")
check("net debt = gross debt − cash", abs(netdebt - (inp["Gross debt"] - inp["Cash and temporary investments"])) < 1, f"${netdebt:,.0f}M")
L = lambda prefix: next(v for k, v in inp.items() if k.startswith(prefix))
surplus, deficit = L("Pension plans in surplus"), L("Pension plans in deficit")
check("10-K pension lines tie out: surplus plans − deficit plans = the filing's net, and the U.S. plan's table = its stated surplus",
      abs((surplus - deficit) - L("Pensions, all plans")) < 1 and abs(L("Pension, U.S. qualified plan: assets") - L("Pension, U.S. qualified plan: obligation") - 366) < 1,
      f"${surplus:,.0f}M − ${deficit:,.0f}M = ${surplus - deficit:,.0f}M (10-K: ${L('Pensions, all plans'):,.0f}M); U.S. qualified ${L('Pension, U.S. qualified plan: assets'):,.0f}M − ${L('Pension, U.S. qualified plan: obligation'):,.0f}M = $366M")
dl = ((deficit + L("Retiree medical") + L("Environmental remediation") + L("Asbestos claims")) * (1 - inp["Cash tax rate"])
      - surplus * (1 - inp["Cash tax rate"]) * L("Count the pension surplus") + L("Brazil tax guarantee")
      + L("Deferred tax on the Temple-Inland") - (L("Timber notes held") - L("Nonrecourse loans")))
check("debt-like items rebuild from the 10-K lines and come off the equity, not the WACC", abs(dl - debtlike) < 1 and abs(
      DCF["Equity value = EV − net debt − debt-like"] - (DCF["Enterprise value"] - netdebt - debtlike)) < 1,
      f"${debtlike:,.0f}M after tax = ${dl:,.0f}M from the lines; equity = EV − ${netdebt:,.0f}M − ${debtlike:,.0f}M")
e26 = M["C14"].value
check("2026E EBITDA = guidance midpoint", abs(e26 - (inp["2026 adj. EBITDA guidance — low"] + inp["2026 adj. EBITDA guidance — high"]) / 2) < 1, f"${e26:,.0f}M vs $3,300M guide midpoint")
lfcf26 = M["C30"].value
check("2026E levered FCF inside the company's Jan guide ($300–500M)", 250 <= lfcf26 <= 550, f"${lfcf26:,.0f}M (guide given before the April cut; the low end is the honest reading)")
div = -M["C31"].value
check("dividend cost ≈ $980M", abs(div - 980) < 20, f"${div:,.0f}M = $1.85 × 529.5M")
check("2026E levered FCF does NOT cover the dividend", lfcf26 < div, f"{lfcf26/div:.2f}× coverage")

# --- terminal value cross-checks (each method implies the other's parameter)
tv_mult = DCF["implied exit multiple of the Gordon TV"]; tv_g = DCF["implied perpetual growth of the exit-multiple TV"]
check("Gordon TV implies an exit multiple between SW and PKG", inp["Smurfit Westrock EV/EBITDA (LTM)"] <= tv_mult <= inp["PKG EV/EBITDA (LTM)"], f"{tv_mult:.1f}x (SW {inp['Smurfit Westrock EV/EBITDA (LTM)']}x, PKG {inp['PKG EV/EBITDA (LTM)']}x)")
check("Exit-multiple TV implies a sane perpetual growth (0–4%)", 0 <= tv_g <= 0.04, f"{tv_g*100:.2f}%")
check("WACC > terminal growth", DCF["WACC"] > inp["Terminal growth rate"], f"{DCF['WACC']*100:.2f}% > {inp['Terminal growth rate']*100:.1f}%")
tv_share = DCF["terminal value as % of EV"]
check("terminal value share of EV disclosed and in the mature-company band (60–90%)", 0.6 <= tv_share <= 0.9, f"{tv_share*100:.0f}%")

# --- three routes agree
dcf_ps, scen_mid_ps, sens_cell = DCF["Value per share ($)"], S["C30"].value, V["D8"].value  # sensitivity row 17% × col 8.5x
check("DCF sheet ≈ Scenarios mid-cycle column (same chain, two layouts)", abs(dcf_ps - scen_mid_ps) < 1.0, f"${dcf_ps:.2f} vs ${scen_mid_ps:.2f}")
check("Sensitivity grid (17% × 8.5x) ≈ DCF within the interpolation error", abs(sens_cell - dcf_ps) < 3.0, f"${sens_cell:.2f} vs ${dcf_ps:.2f}")
check("Bear < Mid < Bull", S["B30"].value < S["C30"].value < S["D30"].value, f"${S['B30'].value:.0f} < ${S['C30'].value:.0f} < ${S['D30'].value:.0f}")

# --- the CFO's sensitivity is reproduced
e0, e25 = V["C16"].value, V["C17"].value  # 2027E EBITDA at 0% and 25% stick
check("25% of $170/ton moves 2027E EBITDA by 42.5 × $9M", abs((e25 - e0) - 42.5 * 9) < 1, f"${e25-e0:,.1f}M vs ${42.5*9:,.1f}M (CFO: $10/ton ≈ $90M)")

# --- Python reproduces the spreadsheet chain (the Monte Carlo's point estimate equals the mid column)
import importlib.util, io, contextlib
src = open("montecarlo.py").read().replace('json.dump(out, open("mc.json", "w"), indent=1)', '')
ns = {}
with contextlib.redirect_stdout(io.StringIO()):
    exec(src, ns)
import numpy as np
# deterministic run of the same chain at the mid inputs
a = dict(stick=0.5, margin=0.17, mult=8.0, wacc=DCF["WACC"], vol=0.01, emea_m=0.11, emea_d=150.0, capex=0.065)
pd_ = a["stick"] * 170 * 9 - 175
rev_na_27 = ns["rev_na_26"] * (1 + a["vol"]) + pd_; rev_emea_27 = ns["rev_emea_26"] * 1.025
na27 = 2425 + pd_ + ns["rev_na_26"] * a["vol"] * 0.30 + (300 - 250) * ns["rev_na_26"] / ns["rev_26"]
emea27 = 950 + a["emea_d"] + (300 - 250) * ns["rev_emea_26"] / ns["rev_26"]
check("Python chain reproduces Scenarios!C21 (2027E EBITDA)", abs((na27 + emea27 - 75) - S["C21"].value) < 1, f"${na27+emea27-75:,.0f}M vs ${S['C21'].value:,.0f}M")

# --- Monte Carlo sanity
check("P(value > price) falls monotonically as P(down-cycle) rises", all(pdown[i]["p_above"] > pdown[i+1]["p_above"] for i in range(len(pdown)-1)), " → ".join(f"{r['p_above']:.0%}" for r in pdown))
check("breakeven P(down-cycle) ≈ 50% (market price = expected value)", any(abs(r["mean"] - 32) < 1.5 and abs(r["p_down"] - 0.5) < 1e-9 for r in pdown), f"mean at p=0.5: ${[r['mean'] for r in pdown if abs(r['p_down']-0.5)<1e-9][0]:.1f}")
check("MC mean sits between the Bear and Bull point estimates", S["B30"].value < mc["mean"] < S["D30"].value, f"${mc['mean']:.1f}")
check("regime split is real: P(>price | down) ≪ P(>price | base)", mc["p_above_price_given_down"] < 0.2 and mc["p_above_price_given_base"] > 0.8, f"{mc['p_above_price_given_down']:.0%} vs {mc['p_above_price_given_base']:.0%}")

print(f"\n{fails} failure(s)")
sys.exit(1 if fails else 0)

"""Monte Carlo over the IP valuation chain (same mechanics as the Scenarios sheet).

The deterministic model gives one number per story. This gives the distribution of
value per share when the four decisive assumptions are uncertain, and answers the
question a bet actually needs answered: P(value > price), the tails, and how big a
position the edge justifies (Kelly).

Everything here mirrors IP_valuation_model.xlsx — same inputs, same bridge, same
terminal mechanics. Run: python3 montecarlo.py  (writes mc.json for the post's d3 chart)
"""
import json
import numpy as np

rng = np.random.default_rng(2026)
N = 20000

# ---------------- fixed inputs, read from the workbook's Inputs sheet so the two can't drift ----------------
from openpyxl import load_workbook
I = {r[0].value: r[1].value for r in load_workbook("IP_valuation_model.xlsx", data_only=True)["Inputs"].iter_rows(min_row=5)
     if r[0].value and isinstance(r[1].value, (int, float))}
price, shares = float(I["Share price ($)"]), float(I["Diluted shares (millions)"])
net_debt = I["Gross debt"] - I["Cash and temporary investments"] + next(v for k, v in I.items() if "debt-like" in k.lower())
rev_2025_na, rev_2025_emea = float(I["2025 pro forma net sales — Packaging Solutions NA"]), float(I["2025 pro forma net sales — Packaging Solutions EMEA"])
rev_2026_h1 = float(I["H1 2026 net sales"])
na_26, emea_26 = float(I["2026 adj. EBITDA guidance — NA (midpoint)"]), float(I["2026 adj. EBITDA guidance — EMEA (midpoint)"])
corp = (I["2026 adj. EBITDA guidance — low"] + I["2026 adj. EBITDA guidance — high"]) / 2 - na_26 - emea_26
price_index_ytd = float(I["2026 containerboard index increase recognized YTD ($/ton, linerboard)"])
price_sens, price_in_2026 = float(I["EBITDA per $1/ton of price (annualized, NA)"]), float(I["Price benefit already inside 2026 NA guidance"])
contrib_margin = I["Contribution margin on incremental NA volume"]
costout_27, infl_27 = float(I["Incremental 80/20 cost-out 2027 (NA + EMEA)"]), float(I["Cost inflation 2027 (OCC, energy, labor, freight)"])
emea_27_delta = float(I["EMEA EBITDA change 2027 (box-price catch-up + cost-out, net)"])
rev_growth_lt, da_pct = I["Revenue growth 2028–2030 (nominal)"], I["D&A as % of revenue 2027+"]
nwc_pct, tax = I["Net working capital investment, % of incremental revenue"], I["Cash tax rate"]
capex_2027, sep_costs = float(I["Capex 2027"]), float(I["Separation / transaction costs expected 2026–27"])
g = I["Terminal growth rate"]

rev_na_26 = rev_2026_h1 * 2 * rev_2025_na / (rev_2025_na + rev_2025_emea)
rev_emea_26 = rev_2026_h1 * 2 * rev_2025_emea / (rev_2025_na + rev_2025_emea)
rev_26 = rev_na_26 + rev_emea_26

# ---------------- the uncertain assumptions, with a cycle ----------------
# Independent draws understate the tails: in a downturn price does NOT stick, demand falls,
# margins compress and the multiple shrinks together. So first draw the regime, then draw the
# assumptions conditional on it.  P(down-cycle over the horizon) is itself a judgment: 30%.
p_down = 0.30
down = rng.random(N) < p_down
# 1. Share of the $170/ton that sticks through 2027. Jason: "more like half, lag is months."
stick = np.where(down, rng.beta(2, 5, N), rng.beta(5, 4, N))          # down: mode ~0.2; base: mode ~0.57
# 2. Mid-cycle NA EBITDA margin: SW 13% floor ... PKG 21% ceiling.
na_margin = np.where(down, rng.triangular(0.11, 0.135, 0.16, N), rng.triangular(0.14, 0.175, 0.21, N))
# 3. Exit EV/EBITDA multiple.
mult = np.where(down, rng.triangular(6.0, 7.0, 8.5, N), rng.triangular(7.0, 8.5, 11.0, N))
# 4. WACC: normal around 7.4%, 50 bp sd, floored well above g.
wacc = np.clip(rng.normal(0.074, 0.005, N), 0.055, 0.10)
# 2027 volume: base +1% +/- 1.5%; down-cycle -2.5% +/- 1.5% (TTOBMA's -2% forecast lives here).
vol = np.where(down, rng.normal(-0.025, 0.015, N), rng.normal(0.01, 0.015, N))
# EMEA mid-cycle margin and the 2027 catch-up also depend on the regime.
emea_margin_mid = np.where(down, rng.triangular(0.08, 0.095, 0.11, N), rng.triangular(0.10, 0.115, 0.13, N))
emea_27_delta = np.where(down, rng.normal(0.0, 75.0, N), rng.normal(150.0, 75.0, N))
# Capex intensity 2028+: the 'assets in the ground' lever.
capex_pct = rng.triangular(0.060, 0.065, 0.078, N)
# Model risk: the structure itself is wrong by some amount we cannot see. Lognormal haircut/uplift on
# the terminal value with 15% sd — Damodaran's margin of safety, made explicit instead of applied at the end.
MODEL_NOISE = 0.15
model_noise = np.exp(rng.normal(0.0, MODEL_NOISE, N))

# ---------------- the chain (vectorized, same as Scenarios sheet) ----------------
price_delta = stick * price_index_ytd * price_sens - price_in_2026
rev_na_27 = rev_na_26 * (1 + vol) + price_delta
rev_emea_27 = rev_emea_26 * (1 + rev_growth_lt)
na_27 = na_26 + price_delta + rev_na_26 * vol * contrib_margin + (costout_27 - infl_27) * rev_na_26 / rev_26
emea_27 = emea_26 + emea_27_delta + (costout_27 - infl_27) * rev_emea_26 / rev_26
ebitda_27 = na_27 + emea_27 + corp
rev_27 = rev_na_27 + rev_emea_27
grow3 = (1 + rev_growth_lt) ** 3
rev_30 = rev_27 * grow3
ebitda_30 = rev_na_27 * grow3 * na_margin + rev_emea_27 * grow3 * emea_margin_mid + corp
ufcf_30 = (ebitda_30 - rev_30 * da_pct) * (1 - tax) + rev_30 * da_pct - rev_30 * capex_pct - rev_30 * rev_growth_lt / (1 + rev_growth_lt) * nwc_pct
ufcf_27 = (ebitda_27 - rev_27 * da_pct) * (1 - tax) + rev_27 * da_pct - capex_2027 - sep_costs * 0.5
pv_interim = (ufcf_27 / (1 + wacc) ** 0.7
              + (ufcf_27 * 2 / 3 + ufcf_30 / 3) / (1 + wacc) ** 1.7
              + (ufcf_27 / 3 + ufcf_30 * 2 / 3) / (1 + wacc) ** 2.7
              + ufcf_30 / (1 + wacc) ** 3.7)
tv = 0.5 * (ufcf_30 * (1 + g) / (wacc - g) + ebitda_30 * mult) * model_noise
ev = pv_interim + tv / (1 + wacc) ** 3.7
value = (ev - net_debt) / shares

# ---------------- what a bet needs to know ----------------
p_above = float((value > price).mean())
pct = {k: float(np.percentile(value, k)) for k in (5, 10, 25, 50, 75, 90, 95)}
ev_value = float(value.mean())
# Treat it as a bet that resolves over ~3 years: return r = value/price - 1 (ignoring the 5.8% dividend,
# which is upside if it survives and a cut if it doesn't). Growth-optimal (Kelly) fraction solves
# max_f E[log(1 + f r)] numerically; we report half of it, which is what practitioners actually size to.
r = value / price - 1
fs = np.linspace(0, 1, 201)
growth = np.array([np.mean(np.log1p(np.clip(f * r, -0.999, None))) for f in fs])
kelly = float(fs[growth.argmax()])
half_kelly = kelly / 2

# what drives it: rank correlation of each input with value
def spearman(a, b):
    ra, rb = a.argsort().argsort(), b.argsort().argsort()
    return float(np.corrcoef(ra, rb)[0, 1])
drivers = {"down-cycle regime": spearman(down.astype(float), value), "price sticks": spearman(stick, value), "NA mid-cycle margin": spearman(na_margin, value),
           "exit multiple": spearman(mult, value), "WACC": spearman(wacc, value),
           "2027 volume": spearman(vol, value), "capex intensity": spearman(capex_pct, value), "model noise": spearman(model_noise, value)}

hist, edges = np.histogram(value, bins=np.arange(0, 100, 2.5))
hist_down, _ = np.histogram(value[down], bins=edges)       # the same bins, split by regime: the two humps
hist_base, _ = np.histogram(value[~down], bins=edges)
out = {"n": N, "price": price, "p_down": p_down, "p_above_price_given_down": float((value[down] > price).mean()), "p_above_price_given_base": float((value[~down] > price).mean()),
       "mean_down": float(value[down].mean()), "mean_base": float(value[~down].mean()), "p_above_price": p_above, "mean": ev_value, "percentiles": pct,
       "kelly": kelly, "half_kelly": half_kelly, "drivers": drivers,
       "hist": {"edges": edges.tolist(), "counts": hist.tolist(), "down": hist_down.tolist(), "base": hist_base.tolist()},
       "median_down": float(np.median(value[down])), "median_base": float(np.median(value[~down])), "share_down": float(down.mean()),
       "model_noise_sd": MODEL_NOISE, "p_loss_half": float((value < price * 0.5).mean()), "p_double": float((value > price * 2).mean())}
json.dump(out, open("mc.json", "w"), indent=1)
print(f"P(value > ${price}) = {p_above:.1%}   mean ${ev_value:.1f}   median ${pct[50]:.1f}")
print("percentiles:", {k: round(v, 1) for k, v in pct.items()})
print(f"P(value < half the price) = {out['p_loss_half']:.1%}   P(value > 2x) = {out['p_double']:.1%}")
print(f"Kelly {kelly:.2f}  half-Kelly {half_kelly:.2f}")
print(f"given down-cycle ({p_down:.0%}): P(>price) {out['p_above_price_given_down']:.1%}, mean ${out['mean_down']:.1f};  given base: P(>price) {out['p_above_price_given_base']:.1%}, mean ${out['mean_base']:.1f}")
print("drivers (rank corr):", {k: round(v, 2) for k, v in drivers.items()})