PythagorasGoal/scripts/research/r0807_piano_netto.py

#!/usr/bin/env python
"""r0807_piano_netto.py — il piano di accumulo rifatto AL NETTO: muro, traiettorie, versamenti.

IL DEBITO. Il 07/08 `r0727_tasse.py` / `r0807_growth_yearly.py` hanno misurato che l'accumulo
composto al lordo sovrastima il capitale del **15.7% a 5 anni, 29.8% a 10, 42.6% a 15**. Il fatto
e' stato scritto in memoria ma i NUMERI del piano non sono stati rifatti — e sono numeri che sono
serviti a decidere: la traiettoria da $600 (\"€250/mese -> 16.2 anni, P(entro 20a) 92%\"), la
tabella \"quanto versare per un orizzonte dato\", la rendita €/g per livello di versamento.

⚠️ E IL MURO STESSO NON E' NEUTRO. `$272.061` viene da una convenzione ASIMMETRICA: il prelievo
viene lordizzato (`€50/g netti -> $27.375/anno lordi` a fisco 33%) ma il capitale che resta
investito compone **senza mai pagare imposte**, in accumulo come in prelievo. Le due meta' del
conto non usano lo stesso fisco. Qui si misura la convenzione COERENTE — le imposte si pagano
dentro il portafoglio, ogni anno, e cio' che si preleva e' gia' netto — e si guarda di quanto il
muro si sposta. I due errori vanno in versi OPPOSTI (la lordizzazione alza il muro, l'assenza di
imposte in prelievo lo abbassa), quindi il segno del netto non e' deducibile: va misurato.

COSA C'E' QUI
  (0) CONTROLLO DI REPLICA — a fisco spento tutto deve riprodurre i numeri pubblicati. Senza
      questo, un numero diverso non si distingue da un bug.
  (1) IL MURO, nelle due convenzioni.
  (2) TRAIETTORIE da $600 per livello di versamento, lorde vs nette, sugli STESSI path.
  (3) QUANTO VERSARE per un orizzonte dato, al netto.
  (4) RENDITA €/g mediana per livello di versamento, al netto.

⚠️ Le assunzioni fiscali sono quelle DICHIARATE il 27/07 e NON sono un parere fiscale: 33%
(sensibilita' a 26%), minusvalenze in carry 4 anni, patrimoniale 0.2%/anno, e il modello tassa la
variazione ANNUA di valore — cioe' un LIMITE SUPERIORE rispetto alla pura realizzazione, vicino
al vero perche' TP01 ribilancia ogni giorno e SKH01 chiude round-trip discreti.

    uv run python scripts/research/r0807_piano_netto.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"))

import r0725_capcurve as CC              # noqa: E402
import r0726_capwall_refresh as WR       # noqa: E402
import r0726_deposits as DP              # noqa: E402
import r0727_lumpsum_split as LS         # noqa: E402
import r0727_tasse as TX                 # noqa: E402

N_PATHS = 4000
BLOCK = 20
SEED = 20260807
SEED_MURO = 725                      # lo stesso di CC.survival: consente la replica esatta
ALIQUOTA = CC.TAX_RATE               # 0.33
VERSAMENTI = (0, 250, 500, 800, 1_000, 2_000)
ORIZZONTI = (10, 15, 20)
CONFIDENZE = (0.50, 0.75, 0.90, 0.95)


# ===========================================================================
# il passo fiscale, uno solo, riusato ovunque
# ===========================================================================
def passo_fiscale(cap: np.ndarray, anno_start: np.ndarray, flusso_anno: np.ndarray,
                  carry: np.ndarray, aliquota: float, patrimoniale: float) -> np.ndarray:
    """Imposte di fine anno, in loco su `cap`. Ritorna l'imposta pagata.

    `flusso_anno` = versamenti MENO prelievi dell'anno: la plusvalenza e' la variazione di valore
    depurata dai movimenti di cassa, quindi la stessa formula serve l'accumulo e il prelievo.
    Identica meccanica a `r0807_growth_yearly.traiettoria` (carry a 4 anni, piu' vecchie prima).
    """
    imposta = np.zeros(len(cap))
    plus = cap - anno_start - flusso_anno
    if aliquota > 0:
        usate = np.minimum(np.maximum(plus, 0.0), carry.sum(axis=1))
        imponibile = np.maximum(plus, 0.0) - usate
        resid = usate.copy()
        for k in range(carry.shape[1] - 1, -1, -1):
            presa = np.minimum(carry[:, k], resid)
            carry[:, k] -= presa
            resid -= presa
        dovuta = imponibile * aliquota
        cap -= np.minimum(dovuta, cap)
        imposta += dovuta
        carry[:] = np.roll(carry, 1, axis=1)
        carry[:, 0] = np.maximum(-plus, 0.0)
    if patrimoniale > 0:
        p = cap * patrimoniale
        cap -= np.minimum(p, cap)
        imposta += p
    return imposta


# ===========================================================================
# (1) il muro: sopravvivenza in prelievo, con e senza fisco dentro
# ===========================================================================
def sopravvivenza(r: np.ndarray, cap0: float, prelievo_anno: float, anni: int,
                  aliquota: float, patrimoniale: float, n_paths: int = N_PATHS,
                  seed: int = SEED_MURO) -> dict:
    """Gemella di `CC.survival` con il fisco DENTRO il portafoglio.

    Con aliquota=0 e patrimoniale=0 deve dare esattamente `CC.survival` — e' il controllo che
    rende credibile la riga col fisco (test `test_piano_netto.py`).
    """
    rng = np.random.default_rng(seed)
    n_days = int(anni * 365)
    paths = CC._boot_paths(r, n_paths, n_days, BLOCK, rng)
    cap = np.full(n_paths, cap0, float)
    vivo = np.ones(n_paths, bool)
    giorno = prelievo_anno / 365.0
    carry = np.zeros((n_paths, TX.CARRY_ANNI))
    anno_start = cap.copy()
    flusso = np.zeros(n_paths)
    for t in range(n_days):
        cap = np.where(vivo, cap * (1.0 + paths[:, t]) - giorno, cap)
        flusso = np.where(vivo, flusso - giorno, flusso)
        vivo &= cap > 0
        if (t + 1) % 365 == 0 and (aliquota > 0 or patrimoniale > 0):
            vivi = vivo.copy()
            passo_fiscale(cap, anno_start, flusso, carry, aliquota, patrimoniale)
            cap = np.where(vivi, cap, 0.0)
            vivo &= cap > 0
            anno_start = cap.copy()
            flusso[:] = 0.0
    fine = np.where(vivo, cap, 0.0)
    return dict(p_survive=float(vivo.mean()), median_end=float(np.median(fine)),
                p10_end=float(np.percentile(fine, 10)))


def muro(r: np.ndarray, prelievo_netto_anno: float, aliquota: float,
         patrimoniale: float) -> tuple[float, float]:
    """(tasso di rendita perpetua, capitale necessario). Perpetua = il prelievo piu' alto con
    P(capitale a 20 anni >= capitale iniziale) >= 90% — stessa definizione del 25/07."""
    lo, hi = 0.0, 0.40
    for _ in range(13):
        mid = (lo + hi) / 2
        s = sopravvivenza(r, 1e6, 1e6 * mid, 20, aliquota, patrimoniale)
        if s["p10_end"] >= 1e6:
            lo = mid
        else:
            hi = mid
    return lo, (prelievo_netto_anno / lo if lo > 0.002 else float("inf"))


# ===========================================================================
# (2-4) accumulo
# ===========================================================================
def accumula(paths: np.ndarray, dep_eur: float, bersaglio: float, aliquota: float,
             patrimoniale: float, start: float = LS.START,
             lump_eur: float = 0.0) -> dict:
    """Accumulo su path GIA' ESTRATTI (confronto appaiato fra regimi fiscali e versamenti).

    Ritorna il capitale finale, il giorno di primo passaggio al bersaglio e il totale versato.
    """
    n_paths, n_days = paths.shape
    dep = dep_eur * CC.EURUSD
    cap = np.full(n_paths, start + lump_eur * CC.EURUSD)
    versato = start + lump_eur * CC.EURUSD
    carry = np.zeros((n_paths, TX.CARRY_ANNI))
    anno_start = cap.copy()
    flusso = np.zeros(n_paths)
    colpito = np.full(n_paths, -1, int)
    for t in range(n_days):
        cap *= (1.0 + paths[:, t])
        if t % 30 == 0 and t > 0 and dep > 0:
            cap += dep
            flusso += dep
            versato += dep
        nuovi = (colpito < 0) & (cap >= bersaglio)
        colpito[nuovi] = t
        if (t + 1) % 365 == 0 and (aliquota > 0 or patrimoniale > 0):
            passo_fiscale(cap, anno_start, flusso, carry, aliquota, patrimoniale)
            anno_start = cap.copy()
            flusso[:] = 0.0
    return dict(cap=cap, colpito=colpito, versato=versato)


def anni_al_traguardo(colpito: np.ndarray) -> tuple[float, float]:
    """(anni mediani al traguardo CONDIZIONATI all'arrivo, P(traguardo entro 20 anni)).

    ⚠️ La mediana e' condizionata — e' la convenzione di `r0726_capwall_refresh.trajectory`, quindi
    i due numeri sono confrontabili con la tabella pubblicata. Va letta insieme alla probabilita':
    «15 anni» su meta' dei path che arrivano non e' «15 anni». Prendere la mediana sull'INTERO
    vettore (con i non-arrivi a -1) darebbe un numero piu' BASSO quanto peggio va la colonna, cioe'
    esattamente il contrario di cio' che si vuole leggere.
    """
    arrivati = colpito[colpito >= 0]
    med = float(np.percentile(arrivati, 50)) / 365.0 if len(arrivati) > 100 else float("nan")
    entro20 = float(((colpito >= 0) & (colpito <= 20 * 365)).mean())
    return med, entro20


def versamento_per_confidenza(paths: np.ndarray, bersaglio: float, conf: float,
                              aliquota: float, patrimoniale: float) -> float:
    """Il versamento mensile minimo con P(traguardo entro l'orizzonte) >= `conf`.
    Bisezione: P e' monotona nel versamento."""
    lo, hi = 0.0, 6_000.0
    if accumula(paths, hi, bersaglio, aliquota, patrimoniale)["colpito"].__ge__(0).mean() < conf:
        return float("inf")
    for _ in range(12):
        mid = (lo + hi) / 2
        p = float((accumula(paths, mid, bersaglio, aliquota, patrimoniale)["colpito"] >= 0).mean())
        if p >= conf:
            hi = mid
        else:
            lo = mid
    return hi


def sezione(t: str) -> None:
    print("\n" + "=" * 100)
    print(f"  {t}")
    print("=" * 100)


def main() -> None:
    print("=" * 100)
    print("  r0807 — IL PIANO AL NETTO: muro, traiettorie e versamenti col fisco d'accumulo dentro")
    print("=" * 100)
    r = DP.deluck_returns()
    netto_anno = CC.TARGET_EUR_DAY * 365 * CC.EURUSD          # €50/g NETTI, senza lordizzare
    lordo_anno = netto_anno / (1 - ALIQUOTA)
    print(f"\n  book live de-luckato ×{DP.DELUCK} · {len(r)} giorni · drift {r.mean() * 365:.1%} · "
          f"vol {r.std() * 365 ** 0.5:.1%}")
    print(f"  bersaglio di rendita: €{CC.TARGET_EUR_DAY:.0f}/g netti = ${netto_anno:,.0f}/anno "
          f"(${lordo_anno:,.0f} se lordizzati al {ALIQUOTA:.0%})")

    # ------------------------------------------------------------------ (0)
    sezione("(0) CONTROLLO DI REPLICA — a fisco spento deve uscire il numero pubblicato")
    perp_mio, muro_mio = muro(r, lordo_anno, 0.0, 0.0)
    n_orig = WR.N_PATHS
    _, perp_2k, muro_2k = WR.perp_and_wall(r, 1.0, lordo_anno)          # default dello script, 2000
    try:
        WR.N_PATHS = N_PATHS
        _, perp_4k, muro_4k = WR.perp_and_wall(r, 1.0, lordo_anno)      # alla taglia pubblicata
    finally:
        WR.N_PATHS = n_orig
    print(f"\n    costante pubblicata (VR.TARGET, citata ovunque) :"
          f"{'':>13}muro ${LS.TARGET:,.0f}")
    print(f"    questa implementazione, fisco OFF, {N_PATHS} path : perpetua {perp_mio:.4%} · "
          f"muro ${muro_mio:,.0f}")
    print(f"    r0726_capwall_refresh ai suoi {n_orig} path        : perpetua {perp_2k:.4%} · "
          f"muro ${muro_2k:,.0f}")
    print(f"    r0726_capwall_refresh portato a {N_PATHS} path      : perpetua {perp_4k:.4%} · "
          f"muro ${muro_4k:,.0f}")
    ok = abs(muro_mio - LS.TARGET) < 1.0 and abs(perp_mio - perp_4k) < 1e-9
    print(f"    -> {'REPLICA ESATTA' if ok else 'NON REPLICA — fermarsi qui'}: due implementazioni "
          "separate, stesso dollaro")
    print(f"\n    ⚠️ Trovato per strada: `perp_and_wall` gira a {n_orig} path e a quella taglia da'")
    print(f"       ${muro_2k:,.0f} invece di ${LS.TARGET:,.0f} — {abs(muro_2k / LS.TARGET - 1):.1%} di rumore Monte Carlo.")
    print("       Il muro pubblicato e' corretto (fu calcolato a 4000), ma la sua TERZA cifra")
    print("       significativa non e' un'informazione: si cita $272k, non $272.061.")
    if not ok:
        raise SystemExit("la replica non torna: qualunque numero col fisco sarebbe indistinguibile da un bug")
    muro_pub, perp_pub = LS.TARGET, perp_mio

    # ------------------------------------------------------------------ (1)
    sezione("(1) IL MURO — la convenzione pubblicata e quella coerente")
    print("\n    (a) PUBBLICATA: prelievo lordizzato, capitale che compone senza mai pagare imposte.")
    print("    (b) COERENTE  : imposte pagate ogni anno dentro il portafoglio, prelievo gia' netto.")
    print(f"\n  {'convenzione':<34}{'prelievo/anno':>15}{'perpetua':>11}{'muro':>14}")
    print(f"  {'(a) pubblicata (fisco solo sul prelievo)':<34}{lordo_anno:>15,.0f}"
          f"{perp_pub:>10.2%}{muro_pub:>14,.0f}")
    righe = {}
    for aliq in (ALIQUOTA, 0.26):
        p, w = muro(r, netto_anno, aliq, TX.PATRIMONIALE)
        righe[aliq] = (p, w)
        print(f"  {f'(b) coerente, aliquota {aliq:.0%}':<34}{netto_anno:>15,.0f}{p:>10.2%}{w:>14,.0f}")
    perp_n, muro_n = righe[ALIQUOTA]
    print(f"\n    Il muro coerente e' ${muro_n:,.0f} contro ${muro_pub:,.0f} pubblicati "
          f"({muro_n / muro_pub - 1:+.1%}).")
    print("    I due errori della convenzione (a) si compensano quasi: lordizzare il prelievo")
    print("    ALZA il muro, non tassare il capitale investito lo ABBASSA. Restano scoperti la")
    print(f"    patrimoniale ({TX.PATRIMONIALE:.1%}/anno) e la non-linearita' della rendita perpetua")
    print("    nel drift, che e' cio' che si vede nella differenza.")

    # ------------------------------------------------------------------ (2)
    sezione(f"(2) TRAIETTORIE DA ${LS.START:.0f} — anni al traguardo, lordo vs netto, sugli STESSI path")
    anni = 25                      # come `r0726_capwall_refresh.trajectory`: 20 troncherebbe la coda
    rng = np.random.default_rng(725)                                  # stesso seed della tabella
    paths = CC._boot_paths(r, 3_000, anni * 365, BLOCK, rng)          # e stessa taglia
    print(f"\n    orizzonte {anni} anni · 3000 path · seed 725 — cioe' ESATTAMENTE la macchina che")
    print("    ha prodotto la tabella pubblicata: la colonna LORDA e' una replica di controllo.")
    print(f"    (LORDO usa il muro pubblicato ${muro_pub:,.0f}, NETTO quello coerente ${muro_n:,.0f})")
    print(f"\n  {'€/mese':>8}{'versato in 20a':>16}   |{'LORDO: anni':>13}{'P(20a)':>9}"
          f"   |{'NETTO: anni':>13}{'P(20a)':>9}{'ritardo':>10}")
    tab2 = []
    for dep in VERSAMENTI:
        L = accumula(paths, dep, muro_pub, 0.0, 0.0)
        N = accumula(paths, dep, muro_n, ALIQUOTA, TX.PATRIMONIALE)
        aL, pL = anni_al_traguardo(L["colpito"])
        aN, pN = anni_al_traguardo(N["colpito"])
        rit = f"{aN - aL:+.1f}a" if not (np.isnan(aL) or np.isnan(aN)) else "—"
        f = lambda a: f"{a:.1f}a" if not np.isnan(a) else "mai"      # noqa: E731
        versato_20a = LS.START + dep * CC.EURUSD * (20 * 365 // 30)
        print(f"  {dep:>8}{versato_20a:>16,.0f}   |{f(aL):>13}{pL:>9.0%}   "
              f"|{f(aN):>13}{pN:>9.0%}{rit:>10}")
        tab2.append(dict(dep=dep, anni_lordo=aL, p_lordo=pL, anni_netto=aN, p_netto=pN))
    print("\n    ⚠️ La colonna «anni» e' la mediana CONDIZIONATA all'arrivo: va letta insieme alla")
    print("       probabilita' accanto. Dove P scende molto (€250/mese al netto) l'anno mediano")
    print("       resta basso perche' descrive solo chi arriva — non e' un miglioramento.")

    # ------------------------------------------------------------------ (3)
    sezione("(3) QUANTO VERSARE PER UN ORIZZONTE DATO — al netto")
    print(f"\n    bersaglio ${muro_n:,.0f} · versamento mensile minimo per la confidenza indicata")
    print(f"\n  {'orizzonte':>10}" + "".join(f"{f'P={c:.0%}':>12}" for c in CONFIDENZE)
          + f"{'tot. versato @P=90%':>22}")
    for anni_h in ORIZZONTI:
        rng_h = np.random.default_rng(SEED)
        ph = CC._boot_paths(r, N_PATHS, anni_h * 365, BLOCK, rng_h)
        celle, tot90 = [], 0.0
        for c in CONFIDENZE:
            v = versamento_per_confidenza(ph, muro_n, c, ALIQUOTA, TX.PATRIMONIALE)
            celle.append(f"€{v:,.0f}/m" if np.isfinite(v) else "irraggiung.")
            if c == 0.90 and np.isfinite(v):
                tot90 = accumula(ph, v, muro_n, ALIQUOTA, TX.PATRIMONIALE)["versato"]
        print(f"  {f'{anni_h} anni':>10}" + "".join(f"{x:>12}" for x in celle)
              + f"{f'${tot90:,.0f}':>22}")

    # ------------------------------------------------------------------ (4)
    sezione("(4) RENDITA €/g MEDIANA — al netto, per livello di versamento")
    print(f"\n    rendita = capitale × perpetua({perp_n:.2%}); le imposte sono gia' dentro la")
    print("    perpetua, quindi il numero e' NETTO e non va lordizzato una seconda volta.")
    print(f"\n  {'€/mese':>8}" + "".join(f"{f'{y}a':>12}" for y in (5, 10, 15, 20))
          + f"{'P(€50/g a 20a)':>17}")
    for dep in VERSAMENTI:
        celle, p50 = [], 0.0
        for y in (5, 10, 15, 20):
            rng_y = np.random.default_rng(SEED)
            py = CC._boot_paths(r, N_PATHS, y * 365, BLOCK, rng_y)
            cap = accumula(py, dep, muro_n, ALIQUOTA, TX.PATRIMONIALE)["cap"]
            rend = cap * perp_n / 365.0 / CC.EURUSD
            celle.append(f"{np.median(rend):,.2f}")
            if y == 20:
                p50 = float((rend >= CC.TARGET_EUR_DAY).mean())
        print(f"  {dep:>8}" + "".join(f"{x:>12}" for x in celle) + f"{p50:>17.1%}")

    # ------------------------------------------------------------------ verdetto
    sezione("VERDETTO — cosa cambia nelle conclusioni gia' usate per decidere")
    d250 = next(x for x in tab2 if x["dep"] == 250)
    d500 = next(x for x in tab2 if x["dep"] == 500)
    print(f"""
    1. Il MURO si sposta poco: ${muro_n:,.0f} contro ${muro_pub:,.0f} ({muro_n / muro_pub - 1:+.1%}).
       La convenzione pubblicata sbagliava DUE volte in versi opposti e le due meta' quasi si
       cancellano — ma per caso, non per costruzione: e' il tipo di errore che si scopre solo
       rifacendo il conto in modo coerente.

    2. E' sulle TRAIETTORIE che il fisco morde, e si legge nella PROBABILITA' piu' che negli anni:
       €250/mese passa da P(entro 20a) {d250['p_lordo']:.0%} a {d250['p_netto']:.0%} (mediana condizionata
       {d250['anni_lordo']:.1f} -> {d250['anni_netto']:.1f} anni); €500/mese da {d500['p_lordo']:.0%} a {d500['p_netto']:.0%}
       ({d500['anni_lordo']:.1f} -> {d500['anni_netto']:.1f} anni). Il livello di versamento che il 26/07 dava come
       «P 92%, il piano funziona» e' quello che al netto diventa una moneta.

    3. Cio' che NON cambia: senza versamenti il capitale-rendita non si raggiunge mai, a
       qualunque lente fiscale. L'ordine di importanza delle leve (versare > quando > quanto
       presto si smette > piatto vs crescente > frequenza) e' invariato, perche' il fisco
       colpisce tutte le colonne allo stesso modo.
    """)


if __name__ == "__main__":
    main()