Monte Carlo Portfolio Simulation in Python: From Random Walks to Retirement Projections

Every retirement calculator asks the same question: “What rate of return do you expect?” It is a trap. A constant 7% compounded over 30 years produces a single number that looks precise and is almost certainly wrong. Markets do not return 7% every year — they return +24%, then -18%, then +31%. The sequence of returns matters more than the average.

Monte Carlo simulation addresses this by running thousands of possible futures with different return sequences. Instead of “your portfolio will be worth $1.2 million,” you get “your portfolio survives in 88% of simulations and fails in 12%.” The second answer is less satisfying but more useful. I built my own simulator after finding that free calculators either used oversimplified assumptions or hid their methodology behind a paid tier.

The Core Simulation Loop

Three components: a return model, a withdrawal rule, and a success criterion:

import numpy as np

def simulate_portfolio(
    initial_balance=1_000_000, annual_withdrawal=40_000, years=30,
    mean_return=0.07, volatility=0.15, inflation=0.03, n_simulations=10_000
):
    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)
        infl = np.random.normal(inflation, 0.01, years)
        for year in range(years):
            wd = annual_withdrawal * (1 + infl[year])
            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()
success_rate = (results[:, -1] > 0).mean()
print(f"Success rate: {success_rate:.1%}")

With $1M, $40K withdrawals, 7% returns, 15% volatility: success rate is roughly 85-90%. The critical insight: the mean outcome (portfolio grows to $2-3M) and failure cases (zero in year 23) coexist in the same distribution.

Modeling Correlation Between Assets

The basic model assumes stocks and bonds move independently — they do not. A correlated model with np.random.multivariate_normal:

def simulate_correlated(
    initial=1_000_000, withdrawal=40_000, years=30,
    stock_ret=0.09, stock_vol=0.17, bond_ret=0.04, bond_vol=0.06,
    corr=0.1, stock_pct=0.6, 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

With a 0.1 stock-bond correlation, bonds reduce portfolio volatility by 20-30% versus 100% equity — improving success rates and reducing failure severity. This is the mathematical justification for the 60/40 portfolio despite bonds having lower expected returns.

Visualizing the Distribution

A single success rate hides everything worth knowing. Here is a three-panel view:

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))

    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 Paths (Fail: {fail:.0%})')
    axes[0].set_ylabel('Portfolio Value ($)')

    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('Percentile Paths'); 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 (${initial:,.0f})')
    axes[2].set_title(f'Terminal Balance\nMedian: ${np.median(terminal_ok):,.0f}')
    axes[2].legend(fontsize=8)

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

The distribution is right-skewed — upside cases extend far beyond downside. Looking only at the median balance is misleading: the worst cases are what you are planning against.

Withdrawal Strategies Beyond Fixed Dollar

The fixed withdrawal ($40K inflation-adjusted) is simple but suboptimal. Three alternatives worth implementing:

• Percentage of portfolio: Withdraw 4% of current balance. You never hit zero, but your income is volatile.
• Guyton-Klinger guardrails: Withdraw 5% initially, skip inflation adjustments when portfolio drops, increase when it rises. Outperforms fixed withdrawals by avoiding drawdowns in bear markets.
• Floor-and-upside: Cover essentials with guaranteed income; variable equity withdrawals for discretionary spending.

Implement these by replacing annual_withdrawal with a pluggable callable:

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):
    rate = prev / balance
    if rate > params['upper']: return prev * 0.9
    if rate < params['lower']: return prev * 1.1
    return prev

Swap the function and rerun. A flexible percentage rule often adds 10-15 percentage points to 30-year success rates — comparable to shifting from 60/40 to 80/20 allocation.