import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def run_simulation_final(n_sim=10000):
    # --- 1. PARAMETRI DI CONFIGURAZIONE ---
    ANNI = 10
    RAL = 35000

    # Percentuali Contributive
    PERC_TFR = 0.0691
    PERC_DIP = 0.012
    PERC_DAT = 0.022

    # Fiscalità
    TAX_TFR_REND = 0.17
    TAX_FONDO_REND = 0.20
    TAX_ETF_GAIN = 0.26
    TAX_BOND_GAIN = 0.125
    TAX_FONDO_EXIT = 0.23
    ALIQUOTA_MEDIA_TFR = 0.25
    ALIQUOTA_MARGINALE = 0.35
    BOLLO = 0.002

    # Costi Strumenti (TER / ISC)
    COSTO_COMETA = 0.001       # ~0.10% (Negoziale)
    COSTO_FPA = 0.012          # ~1.20% (Fondo Aperto Azionario)

    # Scenari di Mercato (Media, Deviazione Standard)
    INFLAZIONE_PARAMS = (0.02, 0.015)
    COMETA_PARAMS = (0.040, 0.06)    # Profilo Bilanciato
    AZIONARIO_PARAMS = (0.080, 0.14) # Azionario Globale
    BOND_PARAMS = (0.030, 0.04)      # Titoli di Stato Euro

    print(f"Esecuzione di {n_sim} simulazioni Monte Carlo... (Orizzonte {ANNI} anni)")

    # --- 2. MOTORE MONTE CARLO ---
    inflazione = np.random.normal(*INFLAZIONE_PARAMS, (n_sim, ANNI))
    r_cometa = np.random.normal(*COMETA_PARAMS, (n_sim, ANNI))
    r_azionario = np.random.normal(*AZIONARIO_PARAMS, (n_sim, ANNI))
    r_bond = np.random.normal(*BOND_PARAMS, (n_sim, ANNI))

    # Inizializzazione vettori
    cometa_m = np.zeros(n_sim)
    cometa_vers = np.zeros(n_sim)
    fpa_m = np.zeros(n_sim)
    fpa_vers = np.zeros(n_sim)
    tfr_m = np.zeros(n_sim)
    tfr_cap = np.zeros(n_sim)
    etf_m = np.zeros(n_sim)
    etf_inv = np.zeros(n_sim)
    bond_m = np.zeros(n_sim)
    bond_inv = np.zeros(n_sim)

    # Ciclo annuale
    for t in range(ANNI):
        q_tfr = RAL * PERC_TFR
        q_dip_lorda = RAL * PERC_DIP
        q_dat = RAL * PERC_DAT

        # A. COMETA
        ing_cometa = q_tfr + q_dip_lorda + q_dat
        cometa_vers += ing_cometa
        gain_cometa = cometa_m * r_cometa[:, t]
        tax_cometa = np.maximum(0, gain_cometa * TAX_FONDO_REND)
        cometa_m = cometa_m + gain_cometa - tax_cometa - (cometa_m * COSTO_COMETA) + ing_cometa

        # B. FPA AZIONARIO
        ing_fpa = q_tfr + q_dip_lorda
        fpa_vers += ing_fpa
        gain_fpa = fpa_m * r_azionario[:, t]
        tax_fpa = np.maximum(0, gain_fpa * TAX_FONDO_REND)
        fpa_m = fpa_m + gain_fpa - tax_fpa - (fpa_m * COSTO_FPA) + ing_fpa

        # C. TFR AZIENDA
        coeff_riv = 0.015 + (0.75 * np.maximum(inflazione[:, t], 0))
        riv_netta = (tfr_m * coeff_riv) * (1 - TAX_TFR_REND)
        tfr_m += riv_netta + q_tfr
        tfr_cap += q_tfr

        # D. INVESTIMENTI PRIVATI
        netto_inv = q_dip_lorda * (1 - ALIQUOTA_MARGINALE)
        etf_inv += netto_inv
        etf_m = etf_m * (1 + r_azionario[:, t]) * (1 - BOLLO) + netto_inv
        bond_inv += netto_inv
        bond_m = bond_m * (1 + r_bond[:, t]) * (1 - BOLLO) + netto_inv

    # --- 3. CALCOLO VALORI NETTI DI LIQUIDAZIONE ---
    def get_netto_fondo(montante, versato):
        perc_rend = np.maximum(0, (montante - versato) / montante)
        imponibile = montante * (1 - perc_rend)
        return montante - (imponibile * TAX_FONDO_EXIT)

    netto_totale_cometa = get_netto_fondo(cometa_m, cometa_vers)
    netto_totale_fpa = get_netto_fondo(fpa_m, fpa_vers)
    netto_totale_tfr = tfr_m - (tfr_cap * ALIQUOTA_MEDIA_TFR)
    netto_totale_etf = etf_m - (np.maximum(0, etf_m - etf_inv) * TAX_ETF_GAIN)
    netto_totale_bond = bond_m - (np.maximum(0, bond_m - bond_inv) * TAX_BOND_GAIN)

    # --- 4. COSTRUZIONE DATAFRAME DETTAGLIATI ---

    # Rinomino le chiavi per renderle univoche nella mega tabella finale
    df_parziale = pd.DataFrame({
        'Cometa 75% (Parz)': netto_totale_cometa * 0.75,
        'FPA 75% (Parz)': netto_totale_fpa * 0.75,
        'No Fondo ETF (Parz)': (netto_totale_tfr * 0.70) + netto_totale_etf,
        'No Fondo Bond (Parz)': (netto_totale_tfr * 0.70) + netto_totale_bond
    })

    df_totale = pd.DataFrame({
        'Cometa 100% (Tot)': netto_totale_cometa,
        'No Fondo ETF (Tot)': netto_totale_tfr + netto_totale_etf,
        'No Fondo Bond (Tot)': netto_totale_tfr + netto_totale_bond
    })

    # --- 5. OUTPUT TABELLARE CON PROBABILITÀ ---

    def print_stats_with_prob(df, benchmark_col, title):
        stats = df.describe(percentiles=[0.05, 0.5, 0.95]).T
        stats = stats[['mean', '5%', '50%', '95%']]
        stats.columns = ['Media (€)', 'Pessimo 5% (€)', 'Mediana (€)', 'Ottimo 95% (€)']

        probs = []
        benchmark_vec = df[benchmark_col].values
        for col in df.columns:
            if col == benchmark_col:
                probs.append("-")
            else:
                win_rate = np.mean(benchmark_vec > df[col].values) * 100
                probs.append(f"{win_rate:.1f}%")
        stats['Prob. Cometa Vince'] = probs

        print("\n" + "="*90)
        print(title)
        print("="*90)
        pd.set_option('display.float_format', '{:,.0f}'.format)
        print(stats.sort_values(by='Media (€)', ascending=False))

    print_stats_with_prob(df_parziale, 'Cometa 75% (Parz)',
                          "SCENARIO A: ANTICIPO CASA (Mantengo Lavoro - Liquidità Parziale)")

    print_stats_with_prob(df_totale, 'Cometa 100% (Tot)',
                          "SCENARIO B: CAMBIO LAVORO (Riscatto Totale - Liquidità Totale)")

    # --- 6. PLOTTING (GRAFICI) ---
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Plot 1: Parziale
    df_parziale.boxplot(ax=axes[0], showfliers=False, patch_artist=True,
                        boxprops=dict(facecolor="lightblue"))
    axes[0].set_title('Scenario A: Anticipo Casa (Parziale)')
    axes[0].set_ylabel('Liquidità Netta (€)')
    axes[0].grid(True, linestyle='--', alpha=0.7)
    axes[0].tick_params(axis='x', rotation=45)

    # Plot 2: Totale
    df_totale.boxplot(ax=axes[1], showfliers=False, patch_artist=True,
                      boxprops=dict(facecolor="lightgreen"))
    axes[1].set_title('Scenario B: Cambio Lavoro (Totale)')
    axes[1].grid(True, linestyle='--', alpha=0.7)
    axes[1].tick_params(axis='x', rotation=45)

    # Nota rossa in alto
    axes[1].text(0.5, 0.95, "Nota: FPA escluso (non riscattabile)",
                 transform=axes[1].transAxes, ha='center', va='top',
                 color='red', fontsize=10, weight='bold',
                 bbox=dict(facecolor='white', alpha=0.8, edgecolor='none'))

    plt.tight_layout()
    plt.show()

    # --- 7. MEGA TABELLA RIASSUNTIVA (RANKING ASSOLUTO) ---

    # Unione dei dataframe per classifica globale
    df_mega = pd.concat([df_parziale, df_totale], axis=1)

    stats_mega = df_mega.describe(percentiles=[0.05, 0.5, 0.95]).T
    stats_mega = stats_mega[['mean', '5%', '50%', '95%']]
    stats_mega.columns = ['Media (€)', 'Pessimo 5% (€)', 'Mediana (€)', 'Ottimo 95% (€)']

    # Aggiungo colonna "Tipo Ritiro" per chiarezza
    stats_mega['Tipo'] = ['Parziale']*len(df_parziale.columns) + ['Totale']*len(df_totale.columns)

    print("\n" + "="*90)
    print("🏆 CLASSIFICA GENERALE LIQUIDITÀ (RANKING ASSOLUTO)")
    print("Confronto unificato di tutti gli scenari (Parziali e Totali)")
    print("="*90)
    print(stats_mega.sort_values(by='Media (€)', ascending=False))

if __name__ == "__main__":
    run_simulation_final()