How to Build a Simple Monte Carlo Retirement Simulator in Python

If you have ever tried to plan for early retirement with a single "assume 7 percent per year" spreadsheet, you already know the frustration. The number the spreadsheet spits out is a single point estimate. Real markets do not deliver single point estimates. They deliver sequences, some of which are wonderful and some of which are terrible.

A Monte Carlo simulator solves this by running thousands of plausible futures and telling you what percentage of them left money on the table at age 90. You do not need a subscription tool for this. You can build the useful minimum in about 40 lines of Python, and the exercise of building it teaches you more about your own plan than any calculator ever will.

What we are actually building

The simulator will take a few inputs: a starting portfolio balance, an annual withdrawal amount in today's dollars, an expected return, a return volatility, an inflation rate, and a time horizon in years. For each simulated future, it will draw a random annual return from a normal distribution, adjust the withdrawal for inflation, subtract, and see if the portfolio is still alive at the end.

Then it runs that entire sequence 10,000 times and reports the success rate.

That is genuinely all it does. The complexity that commercial tools add on top of this base is worth having, but the base is what teaches you why the 4 percent rule works most of the time and fails the specific times it fails.

The dependencies

Only two, both in the standard scientific Python stack:

import numpy as np
import matplotlib.pyplot as plt

If you do not have these, pip install numpy matplotlib gets you there. Documentation for both is at numpy.org and matplotlib.org respectively.

The core simulation loop

Here is the whole thing, deliberately compact so you can see the entire logic on one screen:

def simulate_one_run(start_balance, annual_withdrawal, years,
                     mean_return, return_std, inflation):
    balance = start_balance
    withdrawal = annual_withdrawal
    for _ in range(years):
        yearly_return = np.random.normal(mean_return, return_std)
        balance = balance * (1 + yearly_return) - withdrawal
        withdrawal *= (1 + inflation)
        if balance <= 0:
            return False
    return True

That is the whole per-run function. Two things worth noting. First, the withdrawal happens after the return is applied, which is a slight simplification. In real life you take money out throughout the year. This matters at the margins but not for the overall picture. Second, the withdrawal grows by the inflation rate every year, which is what "inflation-adjusted withdrawal" actually means in the 4 percent rule literature.

Running it 10,000 times

def run_monte_carlo(start_balance, annual_withdrawal, years,
                    mean_return, return_std, inflation, trials=10000):
    successes = 0
    for _ in range(trials):
        if simulate_one_run(start_balance, annual_withdrawal, years,
                            mean_return, return_std, inflation):
            successes += 1
    return successes / trials

Call it with your own numbers:

success_rate = run_monte_carlo(
    start_balance=1_000_000,
    annual_withdrawal=40_000,
    years=40,
    mean_return=0.07,
    return_std=0.15,
    inflation=0.03,
)
print(f"Success rate: {success_rate * 100:.1f}%")

For a million-dollar portfolio with a 4 percent withdrawal over 40 years, expect a number in the low 80s to high 80s. That is the "4 percent for 40 years is not as safe as 4 percent for 30 years" result the FIRE community has been rediscovering for a decade. The Bogleheads community has hundreds of threads working through exactly this comparison.

Adding the 10th percentile output

Success rate is one number. The distribution is much more informative. Rewrite to track final balances:

def run_monte_carlo_full(start_balance, annual_withdrawal, years,
                         mean_return, return_std, inflation, trials=10000):
    final_balances = []
    for _ in range(trials):
        balance = start_balance
        withdrawal = annual_withdrawal
        for _ in range(years):
            yearly_return = np.random.normal(mean_return, return_std)
            balance = balance * (1 + yearly_return) - withdrawal
            withdrawal *= (1 + inflation)
            if balance <= 0:
                balance = 0
                break
        final_balances.append(balance)
    return np.array(final_balances)

Now you can query specific percentiles:

results = run_monte_carlo_full(1_000_000, 40_000, 40, 0.07, 0.15, 0.03)
print(f"Median final balance: ${np.median(results):,.0f}")
print(f"10th percentile: ${np.percentile(results, 10):,.0f}")
print(f"Success rate: {(results > 0).mean() * 100:.1f}%")

The 10th percentile output is the number that actually matters. It is what happens in the worst 10 percent of futures. If your 10th percentile is a positive number, your plan survives even ugly sequences. If it is zero, your plan has a real chance of running out.

Plotting the distribution

plt.hist(results / 1_000_000, bins=50, edgecolor='black')
plt.xlabel('Final balance (millions)')
plt.ylabel('Number of simulated runs')
plt.title('Distribution of Retirement Outcomes')
plt.axvline(0, color='red', linestyle='--', label='Portfolio depleted')
plt.legend()
plt.show()

The shape of the resulting histogram is more useful than any single number. You will see a big lump in the middle around your median outcome, a fat tail to the right where sequences went great, and a spike near zero where things went poorly. The spike near zero is your failure rate visualized.

What this simulator does not model, and why that matters

Three intentional omissions to be honest about.

First, this uses a normal distribution for returns. Real market returns have fatter tails than normal, meaning extreme years happen more often than the model predicts. Fixes exist (log-normal, historical bootstrap, mixture distributions) and are worth adding once the base works. The classic reference on this is any of the volatility analysis at the SciPy documentation for statistical distributions.

Second, it assumes returns are independent year to year. They are not, exactly. There is some autocorrelation, and there are regime shifts. A more sophisticated model uses historical block bootstrap sampling instead of independent draws.

Third, and most importantly, it does not model the "sequence-of-returns risk" that early retirees actually care about most. To model that specifically, add a check: what happens if the first five years are all in the bottom quintile of returns. Rerun the simulation with that constraint and compare. The gap between the two success rates is your personal sequence-risk exposure.

Sanity-checking against a reference tool

Once you have your numbers, run the same inputs through a well-tested calculator to see if you are in the right neighborhood. The free tools at EvvyTools include a personal-finance section that will let you sanity-check the output against a professional Monte Carlo implementation. Numbers should be within a few percentage points if your model is correct.

Where to take this next

Two useful extensions once the base works.

First, add asset allocation. Model stocks and bonds as separate return streams with different mean and volatility parameters, then compute the portfolio return as a weighted sum. This lets you actually compare 100/0 versus 80/20 versus 60/40 for your specific horizon.

Second, add spending guardrails. Instead of a fixed inflation-adjusted withdrawal, reduce the withdrawal by 10 percent when the portfolio falls below a threshold. This is roughly the Guyton-Klinger approach, and it dramatically improves success rates in bad sequences.

For a plain-English walkthrough of what these results mean and how to translate them into actual decisions about your retirement plan, see the guide at EvvyTools. It covers the interpretation side without diving further into code.

The point of building it yourself

A calculator someone else built is a black box. When it says "78 percent success rate," you have to trust the model. When you built the model, you know exactly what "78 percent" means: 7,800 out of 10,000 simulated futures had a positive balance at year 40, given the specific return, volatility, and inflation assumptions you plugged in. If you change those assumptions, you can see immediately how much the result depends on them.

That is a very different relationship to have with your retirement plan than the one most people have. And it takes about an afternoon.