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.

warning

The 4% rule emerged from the Trinity Study — a similar simulation on US historical data. It assumes 30-year retirement at 60/40 allocation. For retirements longer than 30 years or withdrawal rates above 3.5%, success rates drop materially. Run your own simulation with your actual numbers. Do not treat 4% as universal.

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.

info

Monte Carlo results are sensitive to inputs. A 1% change in mean return or 2% in volatility shifts success rates by 10+ percentage points. The value of building your own simulator is not getting the "right" number — nobody knows future returns — but testing how sensitive your plan is to each assumption. If your retirement only works at 7% returns, it does not actually work.

---

Originally published at pickuma.com. Subscribe to the RSS or follow @pickuma.bsky.social for new reviews.