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.
