用Python手撸蒙特卡洛模拟：别再把退休金押在“7%年化”上了

你用过任何一种退休计算器吗？它们都会问同一个问题：“你预期年化收益率是多少？” 这是个陷阱。假如你输入7%，它会直接算出一个精确到小数点后的数字——30年后你有多少钱。但真实市场不是这样的。市场今年涨24%，明年跌18%，后年又涨31%。收益出现的先后顺序比平均值更重要。

蒙特卡洛模拟解决的就是这个问题：它不给你一个“魔法数字”，而是跑几千次甚至上万次不同的未来，每次都用不同的收益序列。最后告诉你：“你的组合有88%的概率撑过30年，12%的概率提前归零”。这个答案没那么爽，但更有用。我因为发现网上免费的模拟器要么假设太粗糙，要么把核心逻辑藏进付费版，所以自己写了一个。

核心模拟循环

只需要三样东西：一个收益模型、一条提款规则、一个成功判定标准。下面是最基础的版本：

import numpy as np

def simulate_portfolio(
    initial_balance=1_000_000,    # 初始本金，1百万美元
    annual_withdrawal=40_000,     # 每年提取金额，4万美元
    years=30,                     # 退休时长，30年
    mean_return=0.07,             # 预期年化收益率 7%
    volatility=0.15,              # 年化波动率 15%
    inflation=0.03,               # 通胀 3%
    n_simulations=10_000          # 模拟次数 1万次
):
    np.random.seed(42)            # 固定随机种子，方便复现
    # 初始化余额矩阵：每行是一次模拟，每列是一年结束时的余额
    balances = np.zeros((n_simulations, years + 1))
    balances[:, 0] = initial_balance

    for sim in range(n_simulations):
        # 生成随机收益率序列（正态分布）
        returns = np.random.normal(mean_return, volatility, years)
        # 生成随机通胀序列（正态分布，标准差1%）
        infl = np.random.normal(inflation, 0.01, years)
        for year in range(years):
            # 当年提取金额 = 基础提款 × 累计通胀
            wd = annual_withdrawal * (1 + infl[year])
            # 余额变化：先扣提款，再乘以（1+收益率）
            balances[sim, year + 1] = (balances[sim, year] - wd) * (1 + returns[year])
            # 如果余额归零，提前结束该次模拟
            if balances[sim, year + 1] <= 0:
                balances[sim, year + 1:] = 0
                break
    return balances

# 跑一下
results = simulate_portfolio()
# 计算成功率：最后一年余额仍大于0的比例
success_rate = (results[:, -1] > 0).mean()
print(f"Success rate: {success_rate:.1%}")

运行结果：在100万本金、每年4万提款、7%收益、15%波动的情况下，成功率大约85%-90%。关键发现：平均值告诉你组合能涨到200-300万，但失败案例在第23年就归零了——这两个情况同时存在。

警告
所谓的“4%法则”来自 Trinity Study——那也是一个基于美国历史数据的模拟，假设30年退休、60/40股债配比。如果你的退休期超过30年，或者提款率高于3.5%，成功率会大幅下降。请用你自己的真实数据跑模拟。不要迷信4%。

加入相关性：股票和债券不是各玩各的

上面假设股票和债券收益独立——实际上它们有相关性。用np.random.multivariate_normal可以模拟两者的联动：

def simulate_correlated(
    initial=1_000_000,             # 初始本金
    withdrawal=40_000,             # 年提款额
    years=30,                      # 退休年限
    stock_ret=0.09, stock_vol=0.17,   # 股票预期收益9%，波动17%
    bond_ret=0.04, bond_vol=0.06,     # 债券预期收益4%，波动6%
    corr=0.1,                      # 股债相关系数 0.1（弱正相关）
    stock_pct=0.6,                 # 股票仓位60%
    n_sims=10_000                  # 模拟次数
):
    means = np.array([stock_ret, bond_ret])
    # 协方差矩阵
    cov = np.array([
        [stock_vol**2, stock_vol * bond_vol * corr],
        [stock_vol * bond_vol * corr, bond_vol**2]
    ])
    np.random.seed(42)
    balances = np.zeros((n_sims, years + 1))
    balances[:, 0] = initial

    for sim in range(n_sims):
        # 生成每年的股票和债券收益率（服从联合正态分布）
        rets = np.random.multivariate_normal(means, cov, years)
        # 按仓位合成组合收益率
        port = stock_pct * rets[:, 0] + (1 - stock_pct) * rets[:, 1]
        for y in range(years):
            balances[sim, y+1] = (balances[sim, y] - withdrawal) * (1 + port[y])
            if balances[sim, y+1] <= 0:
                balances[sim, y+1:] = 0
                break
    return balances

当股债相关系数为0.1时，债券能让组合的波动率比纯股票降低20%-30%——这就提高了成功率，而且失败时的亏空程度也更小。这就是为什么60/40组合在数学上合理，即使债券预期收益更低。

把结果画出来：只看一个成功率什么也看不出来

一个数字（比如87%）掩盖了所有关键信息。下面三张图能告诉你更多：

import matplotlib.pyplot as plt

def plot_simulation(results, initial, years):
    terminal = results[:, -1]                  # 每次模拟的最终余额
    terminal_ok = terminal[terminal > 0]       # 只取成功的情况
    fail = (results[:, -1] == 0).mean()        # 失败率

    fig, axes = plt.subplots(1, 3, figsize=(16, 5))

    # 左图：随便画100条路径，看趋势
    for i in range(min(100, len(results))):
        axes[0].plot(results[i], alpha=0.1, color='steelblue', lw=0.5)
    axes[0].axhline(initial, color='gray', ls='--', alpha=0.5)
    axes[0].set_title(f'100条路径 (失败率 {fail:.0%})')
    axes[0].set_ylabel('组合价值 ($)')

    # 中图：百分位数路径（5th, 25th, 50th, 75th, 95th）
    for p in [5, 25, 50, 75, 95]:
        axes[1].plot(np.percentile(results, p, axis=0), label=f'{p}th', lw=1.2)
    axes[1].axhline(initial, color='gray', ls='--', alpha=0.5)
    axes[1].set_title('百分位数路径')
    axes[1].legend(fontsize=8)

    # 右图：最终余额分布直方图（只画成功样本）
    axes[2].hist(terminal_ok, bins=50, color='steelblue', edgecolor='white')
    axes[2].axvline(initial, color='gray', ls='--', label=f'初始本金 ${initial:,.0f}')
    axes[2].set_title(f'最终余额分布\n中位数: ${np.median(terminal_ok):,.0f}')
    axes[2].legend(fontsize=8)

    plt.tight_layout()
    plt.savefig('monte_carlo_results.png', dpi=150)
    print(f"5th percentile: ${np.percentile(terminal, 5):,.0f}  "
          f"50th: ${np.median(terminal):,.0f}  "
          f"95th: ${np.percentile(terminal, 95):,.0f}  "
          f"失败率: {fail:.1%}")

分布是右偏的——成功的好结局远远大于失败的坏结局。只看中位数会让你误以为很安全，但你真正需要对抗的是最坏的情况。

提款策略：别死磕固定金额

每年提4万（随通胀调整）是最简单的，但不是最优的。下面三种值得实现：

• 按比例提款：每年提当前余额的4%。你永远不会归零，但每年的收入会波动。
• Guyton-Klinger护栏规则：初始提5%，当组合下跌时不调整通胀，当上涨时增加提款。它能避免在熊市里提太多，效果通常比固定提款好。
• 保底+上限：用确定性收入覆盖基本生活，多余的钱才靠股票提款。

实现方法很简单：把提款逻辑变成一个可插拔的函数：

def fixed_wd(year, balance, params):
    """固定金额提款（随通胀增长）"""
    return params['initial'] * (1 + params['inflation']) ** year

def pct_wd(year, balance, params):
    """按当前余额比例提款"""
    return balance * params['rate']

def guardrails_wd(year, balance, prev, params):
    """Guyton-Klinger护栏：根据前一年提款额和当前余额调整"""
    rate = prev / balance           # 当前提款率
    if rate > params['upper']:      # 提款率过高，减少10%
        return prev * 0.9
    if rate < params['lower']:      # 提款率过低，增加10%
        return prev * 1.1
    return prev                     # 正常范围内不变

把提款函数换掉，重新跑一遍。一个灵活的按比例规则往往能让30年成功率提高10-15个百分点——效果相当于把仓位从60/40换成80/20。

提示
蒙特卡洛结果对输入参数非常敏感。平均收益率差1%、波动率差2%，成功率就能差10个百分点以上。自己写模拟器的价值不在于找到“正确”的数字——没人知道未来的收益——而在于测试你的计划对每个假设有多脆弱。如果你的退休计划只在7%收益下才能成功，那它本质上就不可靠。