#!/usr/bin/env python3
"""The model behind register entry 2026-09-003, as published. Run it; it prints what the entry says.

    python3 scripts/forecasts/2026-09-003_mc.py

    P(YES) = R x K x L + (1 - R x K x L) x O

    R  Starship is credible enough by mid-to-late 2027 that SpaceX sells 2029 seats to outside customers
    K  SpaceX keeps rideshare and moves it to Starship
    L  given R and K, a per-kg price is listed on spacex.com before 2028-01-01
    O  any other route to a listing

Each factor is a judgement, not a measurement, so each is drawn from a Beta distribution centred
on the judged value (mean m, concentration k: alpha = m k, beta = (1 - m) k). The factors move
together: if Starship works, the rest become more likely. A Gaussian copula carries that, with
correlation 0.3 between R and K and 0.15 between each of them and L. O is independent.

Needs numpy and scipy. Seed 7, 200,000 draws: the output is the same on every run.
"""
import numpy as np
from scipy import stats

ENTRY = "2026-09-003"
SEED, DRAWS = 7, 200_000
FACTORS = [("R", 0.45, 12), ("K", 0.50, 10), ("L", 0.65, 12), ("O", 0.04, 25)]     # name, mean, concentration
RHO = {("R", "K"): 0.30, ("R", "L"): 0.15, ("K", "L"): 0.15}
STEP = 0.15                                                                       # sensitivity: each factor up and down by this


def point(v):
    rkl = v["R"] * v["K"] * v["L"]
    return rkl + (1 - rkl) * v["O"]


def draws(rng, corr):
    z = rng.multivariate_normal(np.zeros(len(FACTORS)), corr, size=DRAWS)
    u = stats.norm.cdf(z)
    d = {name: stats.beta.ppf(u[:, i], m * k, (1 - m) * k) for i, (name, m, k) in enumerate(FACTORS)}
    rkl = d["R"] * d["K"] * d["L"]
    return rkl + (1 - rkl) * d["O"]


def main():
    names = [f[0] for f in FACTORS]
    corr = np.eye(len(FACTORS))
    for (a, b), r in RHO.items():
        corr[names.index(a), names.index(b)] = corr[names.index(b), names.index(a)] = r
    rng = np.random.default_rng(SEED)
    independent = draws(rng, np.eye(len(FACTORS)))
    correlated = draws(rng, corr)
    means = {name: m for name, m, _ in FACTORS}
    lo, hi = np.quantile(correlated, [0.10, 0.90])
    print("entry %s · seed %d · %d draws" % (ENTRY, SEED, DRAWS))
    print("point estimate, factors at their means   %.3f" % point(means))
    print("independent mean                         %.3f" % independent.mean())
    print("correlated mean                          %.3f   -> %d%%" % (correlated.mean(), round(correlated.mean() * 100)))
    print("80%% range                                %d%% to %d%%" % (round(lo * 100), round(hi * 100)))
    print("sensitivity (each factor +/- %.2f, others at their means):" % STEP)
    for name in ("R", "K", "L"):
        ends = [point(dict(means, **{name: means[name] + s})) for s in (-STEP, STEP)]
        print("  %s  %d%% to %d%%" % (name, round(ends[0] * 100), round(ends[1] * 100)))
    print("correlated_mean=%.4f" % correlated.mean())      # the line scripts/check_forecast_models.py reads


if __name__ == "__main__":
    main()