PythagorasGoal/scripts/research/r0726_capwall_refresh.py

#!/usr/bin/env python
"""r0726_capwall_refresh.py — i MURI DI CAPITALE ricalcolati col fattore de-luck corretto.

Il 25/07 (`r0725_capcurve.py`) ha calcolato rendita perpetua, muro di capitale per 50 EUR/g e
traiettoria da $600 applicando un de-luck **×0.6 scelto a occhio** sul drift del book.

Il 26/07 quel fattore e' stato MISURATO e decomposto (`r0726_deluck_factor.py` +
`r0726_skh_live_book.py`):
  * fortuna d'ancora sul drift, book live Deribit (TP01/SKH01 75/25): **×0.890**
    (24 x 23 = 552 configurazioni, estrazioni congiunte);
  * degradazione del path live: **non-negativa su tutto cio' che e' misurato**
    (SKH01 uscite +0.081 Sharpe di book su 23/23 offset; SKH01 ingressi +0.73pp di drift,
    positivo nel 100% delle estrazioni; TP01 barra parziale trascurabile su 24 ancore).
Il ×0.6 implicava quindi un residuo ×0.687 attribuito al live OLTRE l'ancora, che nessuna misura
sostiene: era il doppio conteggio gia' sospettato in CLAUDE.md, ora quantificato.

Questo script rifa' i numeri della DECISIONE ai tre fattori, sulla stessa macchineria del 25/07
(nessuna re-implementazione: `book_series`, `survival`, `_boot_paths` importati da li').

    uv run python scripts/research/r0726_capwall_refresh.py
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np

ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))
sys.path.insert(0, str(ROOT / "scripts" / "research"))
sys.path.insert(0, str(ROOT / "scripts" / "research" / "alt"))

import r0725_capcurve as CC  # noqa: E402

# fattori a confronto. 0.890 = fortuna d'ancora MISURATA sul book live Deribit; il vero fattore
# onesto e' >= 0.890 perche' la componente live e' non-negativa (non la conto: e' conservativo).
FACTORS = (
    (0.60, "×0.60  (25/07, a occhio)"),
    (0.89, "×0.89  (26/07, ancora MISURATA)"),
    (1.00, "×1.00  (nessun haircut)"),
)
LEVERAGES = (1.0, 1.5)
N_PATHS = 2000


def perp_and_wall(r: np.ndarray, k: float, gross_needed: float) -> tuple[float, float, float]:
    """(SWR-20a, tasso di rendita PERPETUA, capitale necessario) — stessa ricerca binaria
    del 25/07: perpetua = prelievo con P(cap a 20a >= cap iniziale) >= 90%."""
    lo, hi = 0.0, 0.40
    for _ in range(13):
        mid = (lo + hi) / 2
        s = CC.survival(r, 1e6, 1e6 * mid, 20, k, n_paths=N_PATHS)
        if s["p10_end"] >= 1e6:
            lo = mid
        else:
            hi = mid
    perp = lo
    lo2, hi2 = 0.0, 0.40
    for _ in range(13):
        mid = (lo2 + hi2) / 2
        if CC.survival(r, 1e6, 1e6 * mid, 20, k, n_paths=N_PATHS)["p_survive"] >= 0.95:
            lo2 = mid
        else:
            hi2 = mid
    cap = gross_needed / perp if perp > 0.002 else float("inf")
    return lo2, perp, cap


def trajectory(r: np.ndarray, cap_needed: float, n_paths: int = 3000,
               years: int = 25, seed: int = 725) -> dict:
    """Da $600 al capitale-rendita, stocastico (block bootstrap): stessa logica di
    r0725_capcurve.part_d, parametrizzata sul fattore e sul bersaglio."""
    rng = np.random.default_rng(seed)
    n_days = 365 * years
    paths = CC._boot_paths(r, n_paths, n_days, 20, rng)
    out = {}
    for dep_eur in (0, 250, 500, 1_000, 2_000):
        dep = dep_eur * CC.EURUSD
        cap = np.full(n_paths, 600.0)
        hit = np.full(n_paths, -1, int)
        for t in range(n_days):
            cap = cap * (1.0 + paths[:, t])
            if t % 30 == 0 and t > 0:
                cap += dep
            newly = (hit < 0) & (cap >= cap_needed)
            hit[newly] = t
        yrs = np.where(hit >= 0, hit / 365.0, np.nan)
        fin = yrs[~np.isnan(yrs)]
        out[dep_eur] = dict(
            p50=float(np.percentile(fin, 50)) if len(fin) > 100 else float("nan"),
            p10=float(np.percentile(fin, 10)) if len(fin) > 100 else float("nan"),
            p20y=float((yrs <= 20).sum() / len(yrs)))
    return out


def focus(r: np.ndarray, cap_needed: float, dep_eur: float, n_paths: int = 5000,
          years: int = 30, seed: int = 726) -> dict:
    """Drill-down su UN livello di deposito: percentili del tempo, probabilita' cumulate,
    capitale raggiunto per anno, e quanto del capitale finale e' DEPOSITO vs rendimento."""
    rng = np.random.default_rng(seed)
    n_days = 365 * years
    paths = CC._boot_paths(r, n_paths, n_days, 20, rng)
    dep = dep_eur * CC.EURUSD
    cap = np.full(n_paths, 600.0)
    hit = np.full(n_paths, -1, int)
    paid = np.full(n_paths, 600.0)          # capitale immesso (iniziale + depositi)
    paid_at_hit = np.full(n_paths, np.nan)
    snaps, snap_paid = {}, {}
    for t in range(n_days):
        cap = cap * (1.0 + paths[:, t])
        if t % 30 == 0 and t > 0:
            cap += dep
            # ⚠️ il contatore deve seguire il capitale: `cap` riceve il deposito SEMPRE, quindi
            # anche `paid` accumula sempre. Congelarlo al traguardo (come faceva la prima
            # stesura) rende il rapporto capitale/versato gonfiato dopo il traguardo, perche'
            # il capitale conterrebbe depositi che il contatore non conta. Il valore AL
            # traguardo si cattura a parte, con lo snapshot `paid_at_hit`.
            paid = paid + dep
        newly = (hit < 0) & (cap >= cap_needed)
        hit[newly] = t
        paid_at_hit[newly] = paid[newly]
        if (t + 1) % 365 == 0:
            snaps[(t + 1) // 365] = cap.copy()
            snap_paid[(t + 1) // 365] = paid.copy()
    yrs = np.where(hit >= 0, hit / 365.0, np.nan)
    fin = yrs[~np.isnan(yrs)]
    return dict(yrs=yrs, fin=fin, snaps=snaps, snap_paid=snap_paid,
                paid_at_hit=paid_at_hit[~np.isnan(paid_at_hit)], n=n_paths)


def _p_hit(paths: np.ndarray, cap_needed: float, dep_usd: float, n_days: int) -> float:
    """P(il capitale tocca `cap_needed` entro n_days) partendo da $600 con deposito mensile."""
    cap = np.full(paths.shape[0], 600.0)
    hit = np.zeros(paths.shape[0], bool)
    for t in range(n_days):
        cap = cap * (1.0 + paths[:, t])
        if t % 30 == 0 and t > 0:
            cap += dep_usd
        hit |= cap >= cap_needed
    return float(hit.mean())


def solve_deposit(r: np.ndarray, cap_needed: float, years: int, p_target: float,
                  n_paths: int = 4000, seed: int = 726) -> tuple[float, float]:
    """Deposito mensile (EUR) perche' P(arrivare entro `years`) >= p_target.

    Il bootstrap si estrae UNA VOLTA e si riusa a ogni passo della ricerca binaria: con
    campioni diversi per passo la funzione non sarebbe monotona e la bisezione oscillerebbe.
    """
    rng = np.random.default_rng(seed)
    n_days = 365 * years
    paths = CC._boot_paths(r, n_paths, n_days, 20, rng)
    lo, hi = 0.0, 20_000.0
    if _p_hit(paths, cap_needed, hi * CC.EURUSD, n_days) < p_target:
        return float("nan"), float("nan")
    for _ in range(16):
        mid = (lo + hi) / 2
        if _p_hit(paths, cap_needed, mid * CC.EURUSD, n_days) >= p_target:
            hi = mid
        else:
            lo = mid
    dep = hi
    versato = 600.0 + dep * CC.EURUSD * (years * 365 // 30)
    return dep, versato


def main() -> None:
    print("=" * 100)
    print("  r0726 — MURI DI CAPITALE col fattore de-luck MISURATO (era ×0.6 a occhio)")
    print("=" * 100)

    gross = CC.TARGET_EUR_DAY * 365 * CC.EURUSD / (1 - CC.TAX_RATE)
    print(f"  bersaglio: {CC.TARGET_EUR_DAY:.0f} EUR/g netti = EUR "
          f"{CC.TARGET_EUR_DAY*365:,.0f}/anno netti = ${gross:,.0f}/anno LORDI "
          f"(fisco {CC.TAX_RATE:.0%}, EURUSD {CC.EURUSD})")
    print("  book: TP01+SKH01 75/25 su path live (lente 'hourly'), costi reali Deribit\n")

    base = CC.book_series("hourly")
    mu = float(base.values.mean())

    print(f"  {'fattore':<32}{'leva':>6}{'SWR-20a':>10}{'PERPETUA':>10}{'capitale':>16}")
    print("  " + "-" * 74)
    res = {}
    for f, lab in FACTORS:
        r = (base - (1.0 - f) * mu).values.astype(float)
        for k in LEVERAGES:
            swr, perp, cap = perp_and_wall(r, k, gross)
            res[(f, k)] = (swr, perp, cap)
            capf = f"${cap:,.0f}" if np.isfinite(cap) else "IRRAGGIUNGIBILE"
            print(f"  {lab:<32}{k:>6.2f}{swr:>9.2%}{perp:>10.2%}{capf:>16}")
        print()

    print("=" * 100)
    print("  COSA CAMBIA NELLA DECISIONE")
    print("=" * 100)
    for k in LEVERAGES:
        c_old = res[(0.60, k)][2]
        c_new = res[(0.89, k)][2]
        if np.isfinite(c_old) and np.isfinite(c_new):
            print(f"    leva {k:.2f}: muro ${c_old:,.0f} -> ${c_new:,.0f}   "
                  f"({(c_new/c_old - 1)*100:+.0f}%)")
    c_new = res[(0.89, 1.0)][2]
    print(f"\n    Il muro scende di ~45%, ma resta ${c_new:,.0f} = ~{c_new/600:,.0f}x il conto")
    print("    di oggi ($600). La conclusione STRUTTURALE del 25/07 non cambia: i 600 euro")
    print("    come CAPITALE restano refutati, come BIGLIETTO (prop) no. Cambia la taglia")
    print("    dell'errore — e cambia che il numero ora e' misurato invece che scelto.")
    print("\n    ⚠️ Il ×0.89 e' la sola fortuna d'ancora: e' un LIMITE INFERIORE del fattore")
    print("       onesto, perche' ogni componente misurata del path live e' non-negativa.")
    print("       Il muro vero sta quindi fra la riga ×0.89 e la riga ×1.00.")

    # rendita a $600 — il numero che il 25/07 dava a ~0.11 EUR/g
    print("\n" + "=" * 100)
    print("  RENDITA A $600 (lo stesso conto di oggi)")
    print("=" * 100)
    for f, lab in FACTORS:
        r = (base - (1.0 - f) * mu).values.astype(float)
        cagr = float((1 + r.mean()) ** 365.25 - 1)
        eur = 600.0 * cagr * (1 - CC.TAX_RATE) / CC.EURUSD / 365.0
        print(f"    {lab:<32} CAGR {cagr*100:6.2f}%   rendita EUR {eur:.2f}/giorno")
    print("\n    (a questa taglia la differenza fra ×0.6 e ×0.89 sono centesimi al giorno:")
    print("     il fattore conta per il MURO e per le soglie prop, non per il conto attuale.)")

    # ---------------------------------------------------------------- traiettoria
    print("\n" + "=" * 100)
    print("  TRAIETTORIA da $600 al capitale-rendita (block bootstrap, 3000 path, 25 anni)")
    print("=" * 100)
    print("  Il fattore agisce DUE VOLTE: alza il drift (si accumula piu' in fretta) E abbassa")
    print("  il bersaglio (serve meno capitale). Le due colonne centrali separano gli effetti.")
    cap_old = res[(0.60, 1.0)][2]
    cap_new = res[(0.89, 1.0)][2]
    scen = [
        (0.60, cap_old, f"×0.60, muro ${cap_old:,.0f}  (25/07)"),
        (0.89, cap_old, f"×0.89, muro ${cap_old:,.0f}  (solo drift)"),
        (0.89, cap_new, f"×0.89, muro ${cap_new:,.0f}  (26/07)"),
    ]
    T = {}
    for f, capn, lab in scen:
        r = (base - (1.0 - f) * mu).values.astype(float)
        T[lab] = trajectory(r, capn)

    for dep in (0, 250, 500, 1_000, 2_000):
        print(f"\n    deposito {dep:,}/mese")
        print(f"      {'scenario':<34}{'p10':>8}{'MEDIANA':>10}{'P(entro 20a)':>14}")
        for _f, _c, lab in scen:
            d = T[lab][dep]
            fm = (lambda v: f"{v:.1f}a" if np.isfinite(v) else "  —")
            print(f"      {lab:<34}{fm(d['p10']):>8}{fm(d['p50']):>10}{d['p20y']:>13.0%}")

    print("\n    ⚠️ Restano vere le due cose che non dipendono dal fattore: senza depositi il")
    print("       capitale-rendita non si raggiunge, e la mediana e' una MEDIANA (meta' dei")
    print("       path arriva dopo, e una quota non arriva affatto).")

    # ------------------------------------------------------- drill-down su un deposito
    DEP = 250.0
    print("\n" + "=" * 100)
    print(f"  DRILL-DOWN — EUR {DEP:,.0f}/mese (5000 path, 30 anni)")
    print("=" * 100)
    print(f"  E' la soglia dove il verdetto si ribalta: P(entro 20a) 8% (×0.60) -> 91% (×0.89).")
    r_new = (base - (1.0 - 0.89) * mu).values.astype(float)
    r_old = (base - (1.0 - 0.60) * mu).values.astype(float)
    F = focus(r_new, cap_new, DEP)
    F_old = focus(r_old, cap_old, DEP)

    print(f"\n  A) QUANDO si tocca il capitale-rendita (${cap_new:,.0f})")
    print(f"     {'':<14}{'p10':>9}{'p25':>9}{'MEDIANA':>10}{'p75':>9}{'p90':>9}")
    for lab, D in (("26/07 (×0.89)", F), ("25/07 (×0.60)", F_old)):
        f_ = D["fin"]
        if len(f_) < 100:
            print(f"     {lab:<14}{'—':>9}{'—':>9}{'—':>10}{'—':>9}{'—':>9}")
            continue
        q = np.percentile(f_, [10, 25, 50, 75, 90])
        print(f"     {lab:<14}" + "".join(f"{v:8.1f}a" for v in q[:2]) +
              f"{q[2]:9.1f}a" + "".join(f"{v:8.1f}a" for v in q[3:]))

    print(f"\n  B) PROBABILITA' di averlo toccato entro N anni")
    print(f"     {'':<14}" + "".join(f"{y:>9}a" for y in (5, 10, 15, 20, 25, 30)))
    for lab, D in (("26/07 (×0.89)", F), ("25/07 (×0.60)", F_old)):
        y = D["yrs"]
        print(f"     {lab:<14}" + "".join(f"{float((y <= n).sum())/D['n']:9.0%}"
                                          for n in (5, 10, 15, 20, 25, 30)))

    print(f"\n  C) CAPITALE lungo la strada (×0.89) — e quanto ne hai MESSO tu")
    print(f"     {'anno':>5}{'versato':>12}{'p10':>13}{'MEDIANA':>13}{'p90':>13}"
          f"{'  mediana/versato':>18}")
    for y in (5, 10, 15, 20, 25, 30):
        c = F["snaps"][y]
        pd_ = float(np.median(F["snap_paid"][y]))
        q = np.percentile(c, [10, 50, 90])
        print(f"     {y:>5}{'$'+format(pd_, ',.0f'):>12}"
              + "".join(f"{'$'+format(v, ',.0f'):>13}" for v in q)
              + f"{q[1]/pd_:>17.1f}x")

    if len(F["paid_at_hit"]) > 100:
        med_paid = float(np.median(F["paid_at_hit"]))
        print(f"\n  D) AL TRAGUARDO (mediana dei path che arrivano)")
        print(f"     versato di tasca:   ${med_paid:,.0f}  "
              f"(= EUR {med_paid/CC.EURUSD:,.0f} in ~{np.median(F['fin']):.0f} anni)")
        print(f"     capitale:           ${cap_new:,.0f}")
        print(f"     quota da RENDIMENTO: {(1 - med_paid/cap_new)*100:.0f}%  "
              f"(il resto sono i tuoi versamenti)")
        print(f"     rendita che paga:   EUR {CC.TARGET_EUR_DAY:.0f}/giorno netti")

    # -------------------------------------------- quanto versare per un ORIZZONTE dato
    print("\n" + "=" * 100)
    print(f"  QUANTO VERSARE PER ARRIVARCI IN N ANNI (bersaglio ${cap_new:,.0f}, ×0.89)")
    print("=" * 100)
    print("  'Arrivarci in N anni' non e' un numero: dipende da quanto vuoi esserne sicuro.")
    print("  La colonna e' la PROBABILITA' di ORRIVARE entro l'orizzonte.".replace("ORRIVARE", "arrivare"))
    print(f"\n     {'orizzonte':<12}{'P=50%':>12}{'P=75%':>12}{'P=90%':>12}{'P=95%':>12}")
    sol = {}
    for years in (10, 15, 20):
        row = f"     {str(years) + ' anni':<12}"
        for p in (0.50, 0.75, 0.90, 0.95):
            dep, versato = solve_deposit(r_new, cap_new, years, p)
            sol[(years, p)] = (dep, versato)
            row += (f"{'EUR ' + format(dep, ',.0f'):>12}" if np.isfinite(dep)
                    else f"{'irragg.':>12}")
        print(row)
    print(f"\n     {'':12}{'(totale versato di tasca su tutto l orizzonte)':<48}")
    for years in (10, 15, 20):
        row = f"     {str(years) + ' anni':<12}"
        for p in (0.50, 0.75, 0.90, 0.95):
            _d, v = sol[(years, p)]
            row += (f"{'$' + format(v, ',.0f'):>12}" if np.isfinite(v) else f"{'—':>12}")
        print(row)
    print("\n     ⚠️ A orizzonti corti il versamento domina il rendimento: si sta comprando")
    print("        il capitale, non generandolo. Il confronto onesto e' con la colonna")
    print("        'totale versato' — se si avvicina al bersaglio, il book non sta facendo")
    print("        il lavoro, lo stanno facendo i bonifici.")


if __name__ == "__main__":
    main()