2026: Atualização Anual da Taxa Segura de Retirada no Brasil

#instalar modulos: pip install numpy, matplotlib
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.ticker as mtick

# ============================
# 1. DADOS HISTÓRICOS (1995–2025)
# ============================

# Valores em DECIMAIS (ex.: -1.24% = -0.0124)

r_ibov = np.array([
    -0.0124,  # 1995
    0.6374,   # 1996
    0.4485,   # 1997
    -0.3346,  # 1998
    1.5193,   # 1999
    -0.1072,  # 2000
    -0.1102,  # 2001
    -0.1701,  # 2002
    0.9734,   # 2003
    0.1781,   # 2004
    0.2771,   # 2005
    0.3293,   # 2006
    0.4311,   # 2007
    -0.41,    # 2008
    0.8266,   # 2009
    0.0103,   # 2010
    -0.181,   # 2011
    0.074,    # 2012
    -0.155,   # 2013
    -0.0291,  # 2014
    -0.1331,  # 2015
    0.3893,   # 2016
    0.2686,   # 2017
    0.1503,   # 2018
    0.3158,   # 2019
    0.0292,   # 2020
    -0.1193,  # 2021
    0.0469,   # 2022
    0.2228,   # 2023
    -0.1036,  # 2024
    0.3395    # 2025
])

r_selic = np.array([
    0.426,  # 1995
    0.242,  # 1996
    0.222,  # 1997
    0.253,  # 1998
    0.228,  # 1999
    0.161,  # 2000
    0.16,   # 2001
    0.175,  # 2002
    0.21,   # 2003
    0.151,  # 2004
    0.174,  # 2005
    0.14,   # 2006
    0.112,  # 2007 
    0.118,  # 2008
    0.095,  # 2009
    0.093,  # 2010
    0.11,   # 2011
    0.082,  # 2012
    0.079,  # 2013
    0.104,  # 2014
    0.125,  # 2015
    0.131,  # 2016 
    0.095,  # 2017
    0.062,  # 2018
    0.058,  # 2019
    0.027,  # 2020
    0.043,  # 2021
    0.117,  # 2022
    0.123,  # 2023
    0.108,  # 2024
    0.145   # 2025
])

inflacao = np.array([
    0.2241,  # 1995
    0.0956,  # 1996
    0.0522,  # 1997
    0.0165,  # 1998
    0.0894,  # 1999
    0.0597,  # 2000
    0.0767,  # 2001
    0.1253,  # 2002
    0.093,   # 2003
    0.076,   # 2004
    0.0569,  # 2005
    0.0314,  # 2006
    0.0446,  # 2007
    0.059,   # 2008
    0.0431,  # 2009
    0.0591,  # 2010
    0.065,   # 2011
    0.0584,  # 2012
    0.0591,  # 2013
    0.0641,  # 2014
    0.1067,  # 2015
    0.0629,  # 2016
    0.0295,  # 2017
    0.0375,  # 2018
    0.0431,  # 2019
    0.0452,  # 2020
    0.1006,  # 2021
    0.0579,  # 2022
    0.0462,  # 2023
    0.0483,  # 2024
    0.0426   # 2025
])

# ============================
# 2. CONFIGURAÇÃO DA CARTEIRA
# ============================
print("|||| AA40- Simulação de Monte Carlo para TSR Máxima |||")

# Parametrizacao pelo usuário
print("Configuração personalizada para a simulação de Monte Carlo")

# carteira
ibov_percentage = float(input("Digite a porcentagem do IBOV na carteira (0-100): ")) / 100
selic_percentage = float(input("Digite a porcentagem do SELIC na carteira (0-100): ")) / 100

# modo de simulacao
mode = input("Escolha o modo de cálculo (TSR/PWR): ").strip().upper() or "TSR"
if mode not in {"TSR", "PWR"}:
    raise ValueError("Modo inválido. Use TSR ou PWR.")

# simulacao
n_years = int(input("Digite o número de anos de aposentadoria (ex.: 30): "))
initial_wealth = float(input("Digite o valor inicial do patrimônio (ex.: 1000000): "))
target_success = float(input("Digite a taxa de sucesso alvo (0-100, ex.: 95): ")) / 100

# repeticoes
n_sims = 10000

# validacao
if ibov_percentage + selic_percentage != 1.0:
    raise ValueError("A soma das porcentagens do IBOV e SELIC deve ser 100%.")

# atualiza carteira
r_5050 = ibov_percentage * r_ibov + selic_percentage * r_selic

print("Aguarde... rodando resultados...")
# ============================
# 3. SIMULAÇÃO COM PERMUTAÇÃO COMPLETA (SORR)
# ============================

np.random.seed(42)  # Reprodutibilidade

def simulate_batch(tsr, n_sims, return_inflation=False):
    """
    Simula n_sims trajetórias com TSR fixa para qualquer período (20-60 anos),
    amostrando com reposição da sequência histórica (SORR).
    Retorna array com valores finais (0.0 se ruína) e opcionalmente os fatores de inflação.
    """
    final_values = np.zeros(n_sims)
    inflation_factors = np.ones(n_sims) if return_inflation else None

    for i in range(n_sims):
        wealth = initial_wealth
        withdrawal = tsr * initial_wealth

        # Amostragem com reposição para permitir qualquer período
        seq_r = np.random.choice(r_5050, size=n_years, replace=True)
        seq_inf = np.random.choice(inflacao, size=n_years, replace=True)
        cumulative_inflation = np.prod(1 + seq_inf)
        if return_inflation:
            inflation_factors[i] = cumulative_inflation

        for t in range(n_years):
            wealth -= withdrawal
            if wealth <= 0:
                final_values[i] = 0.0
                break
            wealth *= (1 + seq_r[t])
            withdrawal *= (1 + seq_inf[t])
        else:
            final_values[i] = wealth

    if return_inflation:
        return final_values, inflation_factors
    return final_values

def success_rate(tsr, n_sims):
    results = simulate_batch(tsr, n_sims)
    return np.sum(results > 0) / n_sims


def pwr_success_rate(tsr, n_sims):
    final_values, inflation_factors = simulate_batch(tsr, n_sims, return_inflation=True)
    required_final = initial_wealth * inflation_factors
    return np.sum(final_values >= required_final) / n_sims

if mode == "TSR":
    # ============================
    # 4. BUSCA BINÁRIA PARA TSR MÁXIMA (96% SUCESSO)
    # ============================

    low, high = 0.00, 0.12  # busca entre 0% e 12% ao ano

    for _ in range(20):     # iterações reduzidas (convergência rápida)
        mid = (low + high) / 2
        sr = success_rate(mid, n_sims)
        if sr >= target_success:
            low = mid
        else:
            high = mid

    tsr_max = low
    print(f"TSR máxima segura ({target_success:.0%} de sucesso): {tsr_max:.4%}")

    # ============================
    # 5. DISTRIBUIÇÃO DO VALOR FINAL COM ESSA TSR
    # ============================

    final_values = simulate_batch(tsr_max, n_sims)

    p5  = np.percentile(final_values, 5)
    p25 = np.percentile(final_values, 25)
    p50 = np.percentile(final_values, 50)
    p75 = np.percentile(final_values, 75)
    p95 = np.percentile(final_values, 95)

    print("Percentis finais (R$):")
    print(f"P5  : R$ {p5:,.0f}")
    print(f"P25 : R$ {p25:,.0f}")
    print(f"P50 : R$ {p50:,.0f}")
    print(f"P75 : R$ {p75:,.0f}")
    print(f"P95 : R$ {p95:,.0f}")

    # ============================
    # 6. GRÁFICO (INÍCIO E FINAL)
    # ============================

    anos = np.array([0, n_years])  # Ano 0 (início) e ano final

    plt.figure(figsize=(12, 6))

    plt.plot(
        anos,
        [initial_wealth, p50],
        label=f"Mediana (50%) : R$ {p50:,.0f}",
        linewidth=2,
        color="blue"
    )
    plt.plot(
        anos,
        [initial_wealth, p25],
        label=f"25% : R$ {p25:,.0f}",
        linestyle="--",
        color="green"
    )
    plt.plot(
        anos,
        [initial_wealth, p75],
        label=f"75% : R$ {p75:,.0f}",
        linestyle="--",
        color="green"
    )
    plt.plot(
        anos,
        [initial_wealth, p5],
        label=f"5% : R$ {p5:,.0f}",
        linestyle=":",
        color="red"
    )
    plt.plot(
        anos,
        [initial_wealth, p95],
        label=f"95% : R$ {p95:,.0f}",
        linestyle=":",
        color="red"
    )

    # Update the title to reflect user-selected portfolio percentages
    plt.title(f"Evolução do Patrimônio – Carteira {ibov_percentage*100:.0f} IBOV /{selic_percentage*100:.0f} SELIC – TSR máx (MC@ {target_success:.0%}): {tsr_max:.2%}")


    plt.xlabel("Ano")
    plt.ylabel("Patrimônio (R$)")
    plt.grid(True)

    # Eixo Y em BRL, sem notação científica
    plt.gca().yaxis.set_major_formatter(mtick.StrMethodFormatter("R$ {x:,.0f}"))

    plt.legend()
    plt.tight_layout()

    # Add watermark
    plt.figtext(0.5, 0.01, "Copyright - AA40", wrap=True, horizontalalignment='center', fontsize=8, alpha=0.5)

    plt.show()
else:
    # ============================
    # BUSCA BINÁRIA PARA PWR (MANUTENÇÃO DO PRINCIPAL REAL)
    # ============================

    low, high = 0.00, 0.12
    for _ in range(20):
        mid = (low + high) / 2
        sr = pwr_success_rate(mid, n_sims)
        if sr >= target_success:
            low = mid
        else:
            high = mid

    pwr_rate = low
    print(f"PWR ({target_success:.0%} de manter o patrimônio real): {pwr_rate:.4%}")

    # Estatísticas do patrimônio corrigido pelo IPCA para a TSR encontrada
    final_values, inflation_factors = simulate_batch(pwr_rate, n_sims, return_inflation=True)
    real_final_values = final_values / inflation_factors

    median_real_final = np.percentile(real_final_values, 50)
    p5_real = np.percentile(real_final_values, 5)
    p95_real = np.percentile(real_final_values, 95)

    median_inflation = np.percentile(inflation_factors, 50)
    required_final_median = initial_wealth * median_inflation

    print("Patrimônio final corrigido pelo IPCA (valores reais de hoje):")
    print(f"P5 real : R$ {p5_real:,.0f}")
    print(f"P50 real: R$ {median_real_final:,.0f}")
    print(f"P95 real: R$ {p95_real:,.0f}")
    print(f"Requisito mediano de preservação real (Valor Inicial corrigido por IPCA): R$ {required_final_median:,.0f}")