Monte Carlo Simulation in Trading: Stress-Testing Strategies and Risk

Monte Carlo simulation in trading is the discipline of replacing a single backtest path with a distribution of plausible paths, so that drawdown, risk of ruin, and metric uncertainty become estimable quantities rather than single lucky draws. The technique is simple to code and easy to misuse: the resampling scheme encodes an assumption about the data- generating process, and the wrong scheme produces confidently wrong tail estimates. This article covers the main methods, the dependence structure they preserve or destroy, the convergence of the tail estimates themselves, and the bias that contaminates every simulation built on a flawed backtest.

What you are actually estimating

A backtest equity curve is one realization of a stochastic process. The final compounded return is roughly order-invariant, but path-dependent statistics — maximum drawdown, time underwater, risk of ruin — depend heavily on the ordering and on the draw of returns. Monte Carlo estimates the sampling distribution of these path statistics. The key insight: the quantity you care about (the 99th-percentile drawdown) is itself an estimate with its own sampling error, and tail quantiles converge slowly, so the number of simulations is not a free parameter.

Method 1: IID bootstrap and what it assumes

The basic bootstrap resamples per-trade or per-period returns with replacement. It treats returns as i.i.d. draws from the empirical distribution — which preserves the marginal distribution (including its fat tails and skew) but destroys all serial dependence. That assumption is the method's defining limitation.

import numpy as np

def bootstrap_drawdowns(returns, n_sims=20000, horizon=None, seed=0):
    rng = np.random.default_rng(seed)
    r = np.asarray(returns)
    h = horizon or len(r)
    max_dd = np.empty(n_sims)
    for i in range(n_sims):
        sample = rng.choice(r, size=h, replace=True)
        equity = np.cumprod(1.0 + sample)
        peak = np.maximum.accumulate(equity)
        max_dd[i] = (equity / peak - 1.0).min()
    return max_dd

def report(max_dd):
    qs = [50, 90, 95, 99]
    return {f"p{q}": np.percentile(max_dd, q) for q in qs}

The distinction between with replacement (project to a new sample, e.g. a longer forward horizon) and without replacement / permutation (reuse the exact trades, isolating pure sequence risk) matters. Permutation answers "how bad could the drawdown have been with the same trades reordered"; replacement answers "what risk should I budget going forward." They are different questions and the permutation final return is fixed while the bootstrap final return varies.

Convergence: how many simulations for a stable tail

Tail quantile estimates have sampling error that you can quantify and should report. The standard error of an estimated q-quantile is approximately sqrt(q(1-q)) / (n f(x_q)), where f is the density at the quantile — which is small in the tail, so tail quantiles are noisy. The robust move is to bootstrap the quantile estimate itself* to get its confidence interval, and to use enough sims that the 99th percentile is stable.

def quantile_ci(max_dd, q=99, n_boot=1000, seed=1):
    rng = np.random.default_rng(seed)
    n = len(max_dd)
    ests = [np.percentile(rng.choice(max_dd, n, replace=True), q) for _ in range(n_boot)]
    return np.percentile(ests, [2.5, 97.5])

A few hundred simulations give a hopelessly noisy p99; 20,000+ is a sane floor for a stable, narrow-CI tail. If the CI on your p99 drawdown is wide, you do not yet have an answer — increase the count before drawing leverage conclusions.

Method 2: block bootstrap for dependence

Real returns have autocorrelation and volatility clustering — calm and stormy periods arrive in runs. The i.i.d. bootstrap shatters these runs and therefore understates the probability of long, grinding drawdowns. The block bootstrap resamples contiguous blocks, preserving short-range dependence. The stationary bootstrap (Politis-Romano) uses random geometric block lengths, which avoids the artificial periodicity of fixed blocks and is the better default.

def stationary_bootstrap(returns, expected_block=20, n_sims=20000, seed=0):
    rng = np.random.default_rng(seed)
    r = np.asarray(returns)
    n = len(r)
    p = 1.0 / expected_block          # geometric block-length parameter
    out = np.empty((n_sims, n))
    for s in range(n_sims):
        idx = np.empty(n, dtype=int)
        idx[0] = rng.integers(n)
        for t in range(1, n):
            if rng.random() < p:
                idx[t] = rng.integers(n)            # start a new block
            else:
                idx[t] = (idx[t - 1] + 1) % n       # continue, wrap around
        out[s] = r[idx]
    return out

Block length is the critical tuning knob: too short reintroduces the i.i.d. problem; too long leaves too few distinct blocks and underestimates variability. Politis-White provide a data-driven rule based on the return autocorrelation; absent that, scale the expected block with the persistence you observe, not a round number.

Method 3: parametric and model-based simulation

You can fit a model and simulate from it — useful for forward projection and for imposing fat tails explicitly. A Student-t with low degrees of freedom captures the leptokurtosis a Gaussian misses; pairing it with a GARCH variance process captures clustering properly, which a static-sigma simulation cannot.

def garch_t_paths(omega, alpha, beta, nu, horizon, n_sims=20000, seed=0):
    rng = np.random.default_rng(seed)
    lr_var = omega / (1 - alpha - beta)
    paths = np.empty((n_sims, horizon))
    for s in range(n_sims):
        var = lr_var
        for t in range(horizon):
            z = rng.standard_t(nu) * np.sqrt((nu - 2) / nu)   # unit-variance t
            shock = np.sqrt(var) * z
            paths[s, t] = shock
            var = omega + alpha * shock**2 + beta * var       # update conditional var
    return paths

The trade-off: parametric simulation is only as honest as the model. Misspecify the tail or the dependence and you inherit a precise, wrong answer. The bootstrap makes no distributional assumption but cannot extrapolate beyond the observed support — it can never simulate a crash worse than the worst one in your sample, which is a serious limitation for tail risk.

What it tells you, and a worked illustration

Output	What it answers
Drawdown distribution	p95/p99 drawdown to budget capital and leverage
Risk of ruin	probability of touching a fatal loss level
Sharpe confidence interval	whether a 1.5 Sharpe is real or could be 0.6
Time-to-recovery distribution	how long underwater periods plausibly last

The recurring lesson: the single backtested drawdown is usually optimistic relative to the simulated median, because the realized path was one draw and you tend to remember the strategies that drew well. Sizing to the backtested max drawdown rather than the simulated p99 systematically under-capitalizes the strategy.

def risk_of_ruin(return_paths, ruin_level=-0.5):
    # return_paths: per-period returns, shape (n_sims, horizon)
    equity = np.cumprod(1.0 + return_paths, axis=1)
    trough = (equity / np.maximum.accumulate(equity, axis=1) - 1.0).min(axis=1)
    return np.mean(trough <= ruin_level)

Ruin is absorbing — once you blow up, favorable long-run expectancy is irrelevant — so professionals target risk of ruin well under 1% over a realistic horizon. Combine the drawdown distribution with the Kelly criterion to choose a leverage fraction whose p99 drawdown is survivable, then validate against risk of ruin directly. The professional workflow inverts the amateur one: fix the worst acceptable outcome first, then let it cap size.

The permutation test for strategy significance

Monte Carlo also answers a different question: is the strategy's edge distinguishable from luck? A permutation (or label-shuffling) test builds the null distribution of a performance metric under "no skill" by repeatedly destroying the timing relationship between signal and return, then asks where the realized metric falls. For a long/short signal, shuffle the signal relative to forward returns; for a strategy on returns, randomize the entry timing. The empirical p-value is the fraction of permutations that beat the real result.

def permutation_pvalue(signal, fwd_returns, n_perm=5000, seed=0):
    rng = np.random.default_rng(seed)
    real = np.mean(signal * fwd_returns)        # realized mean strategy return
    sig = np.asarray(signal)
    count = 0
    for _ in range(n_perm):
        if np.mean(rng.permutation(sig) * fwd_returns) >= real:
            count += 1
    return (count + 1) / (n_perm + 1)

This is more honest than a t-test when returns are non-normal and autocorrelated, because it makes no distributional assumption — but note it tests this strategy in isolation. If the strategy was selected from many candidates, the single-strategy p-value is optimistic and you still need a multiple-testing adjustment.

Preserving cross-asset structure

For a multi-asset book, resampling each asset's returns independently destroys the cross-correlation that drives portfolio risk — and as discussed, that correlation rises in stress. Resample with the columns kept together: draw the same time indices across all assets (multivariate bootstrap), so each simulated path preserves the contemporaneous co-movement, including the joint tail. With the stationary bootstrap, draw block start points once and apply them to every asset's series in lockstep. This keeps the simulated drawdowns honest about the scenario where everything sells off together, which independent resampling makes spuriously rare.

The bias that contaminates everything

Monte Carlo quantifies risk given the input returns, so every bias in the source backtest is faithfully propagated and often amplified. Three to police:

• Look-ahead bias. If the backtested trades used future information, the return

distribution is too good, and so is every simulation. No resampling fixes a poisoned input.

• Survivorship bias. Backtesting on a universe that dropped delisted names inflates

the mean and thins the left tail; the bootstrap then cannot draw the catastrophes that the survivors avoided. See backtesting biases.

• Selection bias. If the strategy was chosen because it had the best backtest among

many, its returns are an upward-biased sample, and the simulation inherits that optimism. Pair Monte Carlo with multiple-testing correction.

Honest limits

Monte Carlo is a coherence check on the risk implied by your model of the world; it cannot warn you about risks the model never contemplated. The bootstrap cannot extrapolate beyond observed history, so its tails are bounded by your worst observed day. Parametric simulation can extrapolate but only as well as its distribution is specified. Both assume the future return distribution resembles the past, which a regime change voids overnight. Use it to convert point estimates into distributions, to size to a survivable tail rather than a lucky path, and to stress inputs — bump volatility, fatten tails, lengthen blocks — and watch whether your sizing still holds. It is an essential companion to VaR, walk-forward optimization, and disciplined risk management — but it is a magnifying glass on your assumptions, not a crystal ball.