Monte Carlo Simulation for Investment & Portfolio Planning

What Monte Carlo simulation actually does

The term "Monte Carlo" refers to a broad class of computational methods that use repeated random sampling to obtain numerical results — a name coined by physicists Stanislaw Ulam and Nicholas Metropolis in 1949, borrowed from the famous casino, since the method depends on chance the way roulette does. In the context of portfolio planning, the methodology is conceptually direct: define a return distribution and covariance structure for the asset classes in the portfolio, draw random correlated return sequences for each period across a chosen time horizon, and compute the resulting portfolio value. Repeat this process thousands of times — typically 10,000 paths is standard — and collect the distribution of terminal outcomes.

The return distributions can be constructed in two principal ways. The first is parametric: assume that returns follow a known distribution (most commonly normal, or Gaussian) characterised by a mean and covariance matrix estimated from historical data or set as forward-looking assumptions. The second is historical bootstrapping: draw actual historical return sequences at random, with replacement, preserving the correlation structure across asset classes in each sampled period. Each approach has trade-offs, covered in full in the methods comparison section below.

The output of 10,000 simulated paths is visualised as a fan chart: a band of possible portfolio trajectories that narrows near inception (all paths start from the same point) and widens with time as the paths diverge. The fan is shaded to show percentile bands — the 10th through 90th percentile corridor, the median line, and typically the 25th and 75th percentile inner band. This visual representation communicates something that a single line fundamentally cannot: the outcome is genuinely uncertain, the uncertainty compounds over time, and planning should account for the full distribution rather than a single expected value.

The covariance structure underlying these simulations benefits enormously from rigorous estimation via Ledoit-Wolf shrinkage, which produces better-conditioned correlation and volatility inputs than the raw sample covariance matrix — and which feeds directly into the Cholesky decomposition step below, the mechanism that actually turns a covariance matrix into correlated random paths.

Why Cholesky decomposition preserves correlation structure

A naive simulation that draws each asset's returns independently — ignoring the fact that equities and bonds, or domestic and international stocks, move together to varying degrees — will systematically misstate portfolio risk. If two assets are 70% correlated but simulated as if uncorrelated, the simulation understates how often they fall together and overstates the diversification benefit the portfolio actually has. Cholesky decomposition is the standard linear-algebra tool that solves this: it factors the covariance matrix into a form that lets independent random draws be transformed into correctly correlated ones.

In practice: for each simulated year, the algorithm draws a vector of independent standard normal random numbers Z — one per asset class, all uncorrelated with each other — then multiplies by the Cholesky factor L and adds the mean return vector μ. The result, R = μ + LZ, is a vector of asset returns that reproduces the exact variances and covariances specified in Σ. Portfolio weights are then applied to R to get that period's portfolio return. Repeating this across the full horizon for one path, and repeating the whole process 10,000 times, produces the full simulated distribution.

The method is named after French military engineer and geodesist André-Louis Cholesky, who developed the technique for solving systems of linear equations arising in surveying calculations in the early 1900s; it was published posthumously in 1924 after his death in the First World War. Its adoption in finance followed decades later — Phelim Boyle's 1977 paper "Options: A Monte Carlo Approach" was among the first to apply Monte Carlo methods with explicit covariance structure to financial valuation, laying groundwork that portfolio simulation tools still use today. Cholesky decomposition is preferred over alternatives (such as eigenvalue decomposition) because it is computationally cheaper and numerically stable for the positive semi-definite covariance matrices typical of asset-return data.

One practical failure mode worth flagging: Cholesky decomposition requires the covariance matrix to be positive semi-definite. A sample covariance matrix estimated from a short historical window, or from a large number of highly correlated assets, can fail this condition due to estimation noise — the decomposition then fails outright or produces unstable results. This is precisely the problem Ledoit-Wolf shrinkage is designed to fix: shrinking the noisy sample covariance matrix toward a well-conditioned structured target guarantees a valid, stable input to the Cholesky step.

Sequence-of-returns risk: why order matters as much as average

Of all the concepts that Monte Carlo simulation makes legible, sequence-of-returns risk is perhaps the most important for investors approaching or in retirement. The arithmetic seems counterintuitive at first: how can two investors with identical 20-year average annual returns end up with dramatically different terminal wealth? The answer is that in the presence of ongoing withdrawals, the order in which returns are received materially affects the outcome in a way that the average return does not capture.

The asymmetry is structural. During the accumulation phase — when an investor is adding contributions each month — poor early returns are partially mitigated by the ability to purchase additional units at depressed prices. A severe drawdown in year two of a 30-year accumulation horizon is painful but partially self-correcting. The distribution phase reverses this entirely. Once an investor begins withdrawing from a portfolio — in retirement, or during any structured drawdown — poor early returns force the sale of units at depressed prices to fund spending. Those units are permanently gone. Those units cannot participate in the subsequent recovery. The portfolio is permanently impaired.

Monte Carlo simulation captures this risk precisely because Monte Carlo simulation generates thousands of distinct return sequences, including sequences where severe drawdowns cluster early in retirement. A deterministic projection using the average return will always show the same optimistic terminal value. The simulation shows the distribution of outcomes, and critically, the simulation shows how many of those thousands of paths ended in portfolio depletion before the planned end date. That number — the probability of ruin — is one of the most consequential pieces of information available to a retiree, and the probability of ruin is invisible in any single-line projection. This is also why tail risk measurement via CVaR is a natural complement to Monte Carlo analysis: where the simulation tells you the frequency of bad outcomes, CVaR tells you the severity.

How to read a Monte Carlo output

The output of a properly constructed Monte Carlo simulation should be read as a probability distribution, not as a prediction. When 10,000 simulations are run for a portfolio and 2,200 of those paths end in ruin, the result is not a precise prediction of a 22% ruin probability. The result says that under the specified return and volatility assumptions, and under the specified spending rate, 22% of plausible historical return sequences would have caused portfolio depletion. This is qualitatively different from a deterministic forecast, and it should be used differently.

The three numbers to focus on are the 10th percentile planning floor, the 50th percentile central case, and the 90th percentile optimistic scenario. The 10th percentile is the most important of the three. The 10th percentile represents the portfolio value at which only 10% of simulated paths performed worse — it is, in practical terms, the floor against which one should stress-test spending plans. A retirement that remains solvent at the 10th percentile is reasonably robust to adverse sequences. A retirement that runs into difficulty at the 30th percentile is fragile.

The "success rate" figure that Monte Carlo tools report — "your plan has a 78% success rate" — means precisely that 78 of 100 simulated paths achieved the stated goal (most commonly, maintaining a positive portfolio balance through a specified retirement age). A success rate of 95% indicates a robust plan. A success rate of 60% requires adjustment — reduced spending, increased savings, extended working years, or some combination. The simulation does not specify which lever to pull; the simulation specifies how hard the lever needs to be pulled.

One practical point: the success rate figure is sensitive to the assumed time horizon and spending rate. A 4% withdrawal rate on a 30-year horizon produces a materially lower success rate than a 3.5% rate on the same horizon, and the relationship is nonlinear in the tail. Small reductions in the withdrawal rate produce disproportionate improvements in the probability of success because small reductions reduce the frequency of worst-case compounding failures. Modest spending flexibility has outsized impact on tail outcomes — one of the most practically valuable results that Monte Carlo analysis produces.

Distributional assumptions and their limits

Every Monte Carlo simulation is only as good as the distribution it draws from, and this is the single most consequential — and most frequently glossed over — design choice in the entire framework. The overwhelming majority of retail and even institutional planning tools assume returns are normally (Gaussian) distributed. This assumption is analytically convenient: it requires estimating only a mean and covariance matrix, it is well understood, and it decomposes cleanly under Cholesky factorisation. It is also, empirically, wrong.

Realised financial returns exhibit fat tails — extreme events occur far more frequently than a normal distribution predicts. The 1987 crash, the 2008 financial crisis, and the March 2020 drawdown all involved single-day or single-month return realisations that a normal distribution assigns a probability close enough to zero to be considered essentially impossible. A simulation built on normality will therefore generate a fan chart whose lower percentiles are systematically too optimistic — the true 10th percentile outcome, drawn from a fatter-tailed distribution, is worse than what a Gaussian-based simulation shows.

The heatmap above illustrates how ruin probability scales jointly with withdrawal rate and time horizon — and it should be read as a floor, not a ceiling, if the underlying simulation assumes normal returns. Three practical responses mitigate the fat-tail problem. First, use a Student's t-distribution in place of the normal, with degrees of freedom calibrated to historical excess kurtosis — this is a straightforward parametric adjustment that meaningfully improves tail realism without abandoning the Cholesky-based correlated-draw framework. Second, run historical bootstrapping alongside the parametric simulation as a cross-check; because it draws real historical sequences, it captures whatever fat-tail behaviour actually occurred in the sample period. Third, overlay deterministic stress scenarios — a repeat of 2008, a repeat of 1973–74 stagflation — that sit outside what either simulation approach would generate organically, since no calibration on a modern sample will produce a structural regime the sample never experienced.

Practical inputs and their impact on simulation quality

The quality of a Monte Carlo analysis is entirely determined by the quality of its inputs. The most consequential input is the return assumption. Forward-looking return estimates are preferable to simply extrapolating historical averages. Current valuations, yield environments, and structural conditions in markets embed meaningful information about prospective returns that the historical average ignores. Using the historical equity premium mechanically in an environment of elevated valuations and compressed credit spreads overstates the expected return and systematically understates the probability of poor outcomes.

Volatility and correlation inputs are the second most important determinant of simulation quality. The covariance matrix is estimated using Ledoit-Wolf shrinkage, a regularisation technique that improves on the sample covariance matrix by shrinking the sample matrix toward a structured constant-correlation target. The sample covariance matrix — computed directly from historical returns — is notoriously noisy when the number of assets is large relative to the number of observations. Ledoit-Wolf shrinkage produces a better-conditioned estimate that reduces the impact of estimation error on simulation outcomes, and as covered above, it is also what guarantees the covariance matrix remains valid for Cholesky decomposition.

Spending rules and withdrawal rates interact with the simulation in important and nonlinear ways. A fixed dollar withdrawal (e.g., €80,000 per year regardless of portfolio value) produces different and generally worse tail outcomes than a variable spending rule (e.g., withdrawing a fixed percentage of current portfolio value, or using a floor-and-ceiling guardrail approach). Variable spending rules effectively build a form of automatic adjustment into the plan — when the portfolio underperforms, spending adjusts downward, reducing the probability of ruin at the cost of spending variability. The choice of spending rule is as important as the return assumption in determining success rates.

Finally, time horizon sensitivity deserves explicit attention. Extending the planning horizon by five years does not increase risk linearly — extending the horizon increases the probability of encountering a sustained adverse sequence and amplifies the compounding effect of early drawdowns. Planning to the median life expectancy is, by construction, planning to fail half the time.

Safe withdrawal rate calibration

The "4% rule" is the most widely cited figure in retirement planning, and it is also the most widely misapplied. William Bengen introduced it in a 1994 Journal of Financial Planning paper, "Determining Withdrawal Rates Using Historical Data," after testing withdrawal rates against every 30-year rolling period in U.S. market history back to 1926 and finding that a 4% initial withdrawal rate, adjusted annually for inflation, survived every historical sequence for at least 30 years for a 50/50 stock-bond portfolio. The 4% figure is not a law of nature — it is the worst-case survival rate observed in one specific historical dataset, for one specific asset allocation, over one specific set of 30-year windows.

Monte Carlo simulation replaces this single historical worst-case with an explicit, calibrated confidence level. Rather than asking "did 4% survive the worst historical sequence," Monte Carlo asks "what withdrawal rate produces a 90% success rate, given my asset allocation and my forward-looking return assumptions." The two questions can produce meaningfully different answers, particularly when forward-looking equity return assumptions are lower than the long-run historical average — a common conclusion when starting valuations are elevated.

The practical discipline is straightforward: never accept a single withdrawal-rate recommendation without knowing which success rate it targets and which return assumptions it embeds. A planner quoting "4% is safe" without specifying the confidence level and the forward-looking capital market assumptions underneath it is quoting a historical artifact, not a calibrated plan.

What Monte Carlo tells you that a financial plan does not

A traditional financial plan — even a sophisticated one — is a deterministic document. Monte Carlo simulation replaces the single projected path with a distribution, and in doing so, Monte Carlo simulation surfaces four categories of insight that the deterministic plan cannot produce.

First, Monte Carlo simulation quantifies tail scenarios and ruin probability. The probability that a portfolio is depleted before the end of the planning horizon is a calibration of how much buffer the plan contains. A plan showing 95% success has very different characteristics than one showing 75% success, even if both show the same median terminal value. The difference is entirely in the distribution of outcomes below the median.

Second, Monte Carlo simulation reveals the optimal spending rate as a function of desired confidence level. Rather than applying a rule-of-thumb withdrawal rate, the simulation allows explicit calibration: what withdrawal rate corresponds to a 90% success rate, a 95% rate, or a 99% rate? The answer depends on the specific asset allocation, return assumptions, and time horizon — it is not a universal constant. The Asset Lens tool is designed to support this analysis across different portfolio compositions.

"A single-line projection tells you what will probably happen. Monte Carlo tells you what could happen — and how often the plan survives it."

— Anton Ladnyi, CFA

Third, Monte Carlo simulation makes explicit the value of spending flexibility. A client who is willing to reduce spending by 10% in response to a poor sequence of returns dramatically improves the probability of success. The simulation quantifies this: spending flexibility is a form of risk management, and spending flexibility can be traded directly against asset allocation risk. A more conservative portfolio paired with flexible spending often outperforms an aggressive portfolio with rigid spending at the tail of the distribution.

Fourth, Monte Carlo simulation provides the analytical foundation for comparing the impact of different risk management strategies — such as incorporating alternative allocations, adjusting asset class weights, or evaluating structured products — on the full distribution of outcomes rather than just the expected return. This connects naturally to the Black-Litterman framework for expressing return views within a disciplined optimisation, where the optimised portfolio is subsequently stress-tested through simulation.

Monte Carlo vs. historical bootstrapping vs. parametric VaR

Monte Carlo simulation is one of three broad approaches to modelling portfolio risk, and the three are frequently confused or treated as interchangeable when they answer materially different questions.

In practice, the three are complementary rather than substitutes. Monte Carlo simulation is the right tool when the question is multi-period and goal-oriented — "what is the probability this portfolio sustains these withdrawals for 30 years?" CVaR is the right tool when the question is single-period tail severity — "how bad is the average outcome in the worst 5% of next year's returns?" Historical bootstrapping is best used as a cross-check on either, since it grounds the analysis in return sequences that actually occurred rather than ones a distributional assumption generates. A rigorous risk process runs all three and treats disagreement between them as a signal to interrogate assumptions further, not as noise to average away.

When NOT to use Monte Carlo simulation

Monte Carlo simulation is powerful but not universally appropriate. Four situations warrant a different tool, or at minimum, heavy caveats on the simulation's output.

Who uses Monte Carlo simulation — institutional context

Monte Carlo simulation is standard practice across several distinct institutional contexts, each applying the technique to a different decision.

Python implementation — Cholesky Monte Carlo simulation

The following implementation generates 10,000 Cholesky-decomposed, correlated return paths for a multi-asset portfolio, applies annual contributions or withdrawals, rebalances to target weights, and reports the percentile fan and success rate. It requires only NumPy.

import numpy as np

def run_monte_carlo(
    weights,            # e.g. np.array([0.6, 0.4]) for 60/40
    mu,                 # annualised expected returns, e.g. np.array([0.07, 0.03])
    cov,                # annualised covariance matrix (n_assets x n_assets)
    start_value=1_000_000,
    withdrawal=40_000,  # fixed annual withdrawal (4% of start_value)
    years=30,
    n_paths=10_000,
    seed=42
):
    rng = np.random.default_rng(seed)
    n_assets = len(weights)

    # Cholesky factor: Sigma = L @ L.T
    L = np.linalg.cholesky(cov)

    terminal_values = np.zeros(n_paths)
    ruin_year = np.full(n_paths, -1)
    path_history = np.zeros((n_paths, years + 1))
    path_history[:, 0] = start_value

    for p in range(n_paths):
        value = start_value
        for t in range(1, years + 1):
            # Draw independent standard normal shocks, then correlate via L
            z = rng.standard_normal(n_assets)
            asset_returns = mu + L @ z          # R = mu + L*Z
            port_return = weights @ asset_returns

            value = value * (1 + port_return) - withdrawal
            path_history[p, t] = max(value, 0)

            if value <= 0 and ruin_year[p] == -1:
                ruin_year[p] = t
                value = 0  # portfolio depleted, stays at zero

        terminal_values[p] = value

    success_rate = np.mean(ruin_year == -1)
    percentiles = {
        p: np.percentile(path_history[:, -1], p)
        for p in [10, 25, 50, 75, 90]
    }

    return {
        "success_rate": success_rate,
        "ruin_probability": 1 - success_rate,
        "terminal_percentiles": percentiles,
        "path_history": path_history,
    }


if __name__ == "__main__":
    # 60/40 global portfolio — illustrative capital market assumptions
    weights = np.array([0.60, 0.40])
    mu = np.array([0.075, 0.035])          # equities 7.5%, bonds 3.5%
    cov = np.array([
        [0.028, 0.001],                     # equity variance ~16.7% vol
        [0.001, 0.0016],                    # bond variance ~4.0% vol
    ])

    result = run_monte_carlo(weights, mu, cov,
                              start_value=1_000_000,
                              withdrawal=40_000,   # 4% withdrawal rate
                              years=30, n_paths=10_000)

    print(f"Success rate:      {result['success_rate']:.1%}")
    print(f"Ruin probability:  {result['ruin_probability']:.1%}")
    for pct, val in result['terminal_percentiles'].items():
        print(f"  {pct}th percentile: ${val:,.0f}")

Three implementation details matter beyond the code itself. First, np.linalg.cholesky will raise a LinAlgError if the covariance matrix is not positive semi-definite — in production code this should be caught and handled by re-estimating the covariance matrix with Ledoit-Wolf shrinkage rather than silently failing. Second, the random number generator should be seeded for reproducibility during testing but re-seeded (or left unseeded) for genuine Monte Carlo runs used in actual planning decisions. Third, this implementation uses a fixed-dollar withdrawal for clarity; production implementations should support the variable and guardrail spending rules discussed in the practical inputs section, since the choice of spending rule materially changes the output.

Limitations: model risk and what Monte Carlo cannot capture

Monte Carlo simulation is a planning tool, not a forecast. Every simulation result is conditional on its assumptions, and those assumptions introduce model risk that can be as consequential as the market risk being modelled. The distributional assumption and fat-tail problem is covered in full above; two further limitations deserve explicit mention.

The first is parameter estimation error. A 1% error in the assumed real return compounds enormously over a 30-year horizon — success rates should be treated as order-of-magnitude indicators, not precise forecasts to the decimal point. The second is structural breaks and regime changes. Simulations calibrated on post-1990 U.S. data will not produce sequences resembling 1970s stagflation or Japan's post-1990 deflationary balance-sheet recession, because those regimes are not represented in the calibration sample regardless of how many paths are simulated. The appropriate response, as discussed above, is to pair Monte Carlo outputs with deterministic scenario analysis for adverse regimes that sit outside the simulation's assumptions.

Common implementation mistakes

A brief history of Monte Carlo simulation