Future-Proof Your Retirement: A Deep Dive into Monte Carlo Simulation ·

Table of Contents

Python for Financial Freedom - This article is part of a series.

Part 1: Python for Financial Freedom: Code Your Way to Wealth

Part 2: This Article

Planning for the future is filled with uncertainty, especially when it comes to your financial independence. How do you know if your nest egg will last? Will your withdrawal strategy hold up against market volatility? Answering these questions with a single number is not only difficult, it’s misleading. The real world operates on probabilities, not certainties.

This is where Monte Carlo simulation comes in. It’s a powerful analytical technique used by financial professionals to model the probability of different outcomes. Instead of predicting a single future, a Monte Carlo simulation runs thousands, or even hundreds of thousands, of possible futures to see what a realistic range of outcomes looks like.

In this post, we’ll do a deep dive into a Python script that uses a Monte Carlo simulation to project the long-term health of a diversified investment portfolio. We’ll break down the code, visualize the results, and, as a financial analyst, I’ll provide a critique on how the model could be made even more robust.

The Anatomy of my Portfolio
#

First, let’s look at the portfolio we’ll be simulating. It’s a diversified mix of four ETFs with different characteristics. I use this mix of assets myself which has done quite well in recent years.

VTI (Vanguard Total Stock Market ETF): 25%
SCHG (Schwab U.S. Large-Cap Growth ETF): 15%
SCHD (Schwab U.S. Dividend Equity ETF): 40%
CLIP (Global X Treasury Capped 1-3 Year Bond ETF): 20%

This portfolio blends total market exposure, growth, dividends, and short-term bonds. Here is the allocation visualized:

The Monte Carlo Simulation Script
#

The Python script below is the engine of our analysis. It models 100,000 possible futures for our portfolio over a 50-year horizon, assuming a starting value of $1,000,000 and an annual withdrawal of 3% of the current balance.

import numpy as np
import pandas as pd
import json

# --- 1. SIMULATION PARAMETERS ---
N_SIMULATIONS = 100_000      # Number of simulated paths
N_YEARS = 50                 # Simulation duration in years
INITIAL_VALUE = 1_000_000.0  # Starting portfolio value in dollars

# Withdrawal parameters (percentage of current balance each year)
WITHDRAWAL_RATE = 0.03  # 3% of current portfolio balance

# Inflation (still used for real return comparison)
INFLATION_RATE = 0.025

# --- 2. PORTFOLIO & PROXY DEFINITIONS ---
ASSETS = {
    'VTI':  {'mean': 0.10, 'sd': 0.17, 'weight': 0.25},
    'SCHG': {'mean': 0.105, 'sd': 0.18, 'weight': 0.15},
    'SCHD': {'mean': 0.115, 'sd': 0.17, 'weight': 0.40},
    'CLIP': {'mean': 0.040, 'sd': 0.03, 'weight': 0.20}
}

SP500_PROXY = {'mean': 0.10, 'sd': 0.16}

CORR_MATRIX = np.array([
    [1.00, 0.95, 0.95, 0.05], # VTI
    [0.95, 1.00, 0.85, 0.05], # SCHG
    [0.95, 0.85, 1.00, 0.05], # SCHD
    [0.05, 0.05, 0.05, 1.00]  # CLIP
])

# --- 3. PREPARE FOR SIMULATION ---
asset_names = list(ASSETS.keys())
mean_returns = np.array([ASSETS[a]['mean'] for a in asset_names])
std_devs = np.array([ASSETS[a]['sd'] for a in asset_names])
weights = np.array([ASSETS[a]['weight'] for a in asset_names])

cov_matrix = np.outer(std_devs, std_devs) * CORR_MATRIX

portfolio_exp_nominal_return = np.sum(mean_returns * weights)
portfolio_exp_real_return = (1 + portfolio_exp_nominal_return) / (1 + INFLATION_RATE) - 1

# --- 4. RUN THE SIMULATION ---
def run_simulation_percentage(n_sims, n_years, initial_value, withdrawal_rate, annual_returns):
    """
    Simulation where withdrawal is always a % of current balance.
    Returns portfolio value history, final values, and withdrawals per year.
    """
    portfolio_history = np.zeros((n_years + 1, n_sims), dtype=np.float64)
    portfolio_history[0, :] = initial_value
    withdrawals_history = np.zeros((n_years, n_sims), dtype=np.float64)

    for year in range(n_years):
        current_values = portfolio_history[year, :] * (1 + annual_returns[year])
        withdrawals = current_values * withdrawal_rate
        portfolio_history[year + 1, :] = current_values - withdrawals
        portfolio_history[year + 1, portfolio_history[year + 1, :] < 0] = 0
        withdrawals_history[year, :] = withdrawals
    
    final_values = portfolio_history[-1, :]
    return portfolio_history, final_values, withdrawals_history

# Generate random returns for both portfolios
multi_asset_returns = np.random.multivariate_normal(mean_returns, cov_matrix, size=(N_YEARS, N_SIMULATIONS))
annual_portfolio_returns = np.dot(multi_asset_returns, weights)
annual_sp500_returns = np.random.normal(SP500_PROXY['mean'], SP500_PROXY['sd'], size=(N_YEARS, N_SIMULATIONS))

# Run simulations
portfolio_history, final_portfolio_values, withdrawals_history = run_simulation_percentage(
    N_SIMULATIONS, N_YEARS, INITIAL_VALUE, WITHDRAWAL_RATE, annual_portfolio_returns
)

# --- 5. FUNCTION TO GENERATE CHART DATA ---
def generate_chartjs_data(portfolio_history, n_years, inflation_rate):
    """
    Calculates percentile data for portfolio values over time and prints
    a JSON object compatible with Chart.js.
    """
    inflation_path = (1 + inflation_rate) ** np.arange(n_years + 1)
    real_history = portfolio_history / inflation_path[:, np.newaxis]
    
    percentiles = [5, 25, 50, 75, 95]
    percentile_data = {}
    for p in percentiles:
        percentile_data[p] = np.percentile(real_history, p, axis=1)

    chart_data = {
        "labels": list(range(n_years + 1)),
        "datasets": [
            { "label": "95th Percentile", "data": percentile_data[95].tolist(), "fill": False, "borderColor": "rgba(153, 102, 255, 0.6)", "tension": 0.1, "borderWidth": 2, "pointRadius": 0 },
            { "label": "25th-75th Percentile Range", "data": percentile_data[75].tolist(), "fill": '+1', "backgroundColor": "rgba(75, 192, 192, 0.2)", "borderColor": "rgba(75, 192, 192, 0.6)", "tension": 0.1, "borderWidth": 2, "pointRadius": 0 },
            { "label": "Median (50th)", "data": percentile_data[50].tolist(), "fill": False, "borderColor": "rgba(255, 99, 132, 1)", "tension": 0.1, "borderWidth": 3 },
            { "label": "25th Percentile", "data": percentile_data[25].tolist(), "fill": False, "borderColor": "rgba(75, 192, 192, 0.6)", "tension": 0.1, "borderWidth": 2, "pointRadius": 0 },
            { "label": "5th Percentile", "data": percentile_data[5].tolist(), "fill": False, "borderColor": "rgba(255, 159, 64, 0.6)", "tension": 0.1, "borderWidth": 2, "pointRadius": 0 }
        ]
    }
    
    print("\n--- CHART.JS DATA ---")
    print(json.dumps(chart_data))

# generate_chartjs_data(portfolio_history, N_YEARS, INFLATION_RATE)

The script uses a correlated, multi-asset model, accounts for inflation, and models a percentage-based withdrawal strategy.

Visualizing 100,000 Futures
#

Running the script gives us a wealth of data. The most powerful way to understand it is to visualize the range of outcomes. The interactive chart below shows the growth of the portfolio in real, inflation-adjusted dollars over 50 years.

How to Read the Chart:

Median (50th Percentile): This is the most probable outcome. 50% of the simulations ended up better than this, and 50% ended up worse. After 50 years, the median real portfolio value is around $4.01 million.
25th-75th Percentile Range: This is the “interquartile range,” representing the middle 50% of all outcomes. It gives you a sense of what is “likely” to happen. There’s a 50% chance your portfolio’s real value will be between $2.19 million and $7.28 million after 50 years.
5th and 95th Percentiles: These represent the outlier scenarios. There’s a 5% chance of a “bad” outcome (ending with ~$908k) and a 5% chance of a “great” outcome (ending with ~$16.88M).

Critically, with this withdrawal strategy, the probability of the portfolio depleting to zero is virtually 0%.

Analysis: How to Make the Model Better
#

How to make the script better? I see a few ways to improve its realism. In the near future, I will publish this code on github and ask for your support to improve it over time. Leave me a comment below for your feedback. Thanks

Static Assumptions vs. Historical Data: The script uses hardcoded, static values for mean returns and standard deviation. A better approach is to calculate these parameters from historical market data (e.g., using the yfinance library to pull 30+ years of data). An even more advanced technique is bootstrapping, where you sample directly from historical returns instead of assuming a clean normal distribution. This better captures real-world phenomena like market crashes.
Lack of Rebalancing: The simulation never rebalances the portfolio. As some assets outperform others, the portfolio’s allocation will drift from our intended 25/15/40/20 split in my example above. A real investor would rebalance periodically (e.g., annually). Adding a rebalancing step to the simulation loop would more accurately model a managed portfolio’s behavior and risk profile over time.

Conclusion
#

A Monte Carlo simulation is an indispensable tool for long-term financial planning. It forces us to abandon the idea of a single, predictable future and instead embrace a range of possibilities. By visualizing thousands of potential outcomes, we can make more informed decisions about our savings, withdrawal strategies, and risk tolerance.

While my model has room for improvement, it clearly demonstrates that a diversified portfolio combined with a sensible withdrawal strategy has a very high probability of sustaining wealth over the long term. The key is to plan for the unpredictable and build a strategy that is robust enough to handle the ride.