Simulation & Monte Carlo

3.5.1 Why simulate

A formula answers one question exactly. A simulation answers almost any question approximately — and, done properly, tells you how approximately.

In 1946, recovering from brain surgery, the mathematician Stanislaw Ulam played solitaire and wondered what fraction of deals could be won. The combinatorics were hopeless; dealing a hundred hands and counting was easy. With John von Neumann he turned the idea into a way to compute neutron diffusion at Los Alamos, and Nicholas Metropolis named it after the casino where Ulam's uncle liked to gamble (Metropolis 1987; LANL 2023). Phelim Boyle brought it to option pricing in 1977 (Boyle 1977).

Monte Carlo simulation estimates a quantity by drawing many random scenarios from a model, computing the quantity in each, and averaging. It is the law of large numbers of Session 3.1 put to work: the average over n independent scenarios converges to the expected value, and its error shrinks like 1/√n. Quants reach for it for three jobs.

The first job is the one formulas are worst at. A stop-loss, a margin call or the maximum drawdown depends on the route the price took, not only on where it ended (§3.5.4). The second is realism: add fat tails, volatility clustering or a strategy's trading rules and the closed form disappears, while a simulation only needs a model it can sample from — or a history it can resample (§3.5.5). The third matters most in this course: formulas have bugs too, so every lab computes its key answers twice, from theory and by simulation, and asserts that they agree — the volatility drag of Session 2.3, the margin-call odds of Session 2.4, the law of large numbers of Session 3.1.

The random numbers are not random. A pseudo-random generator is a deterministic function of its seed. NumPy's default generator, PCG64, has a period of 2¹²⁸, and Generator.spawn derives child streams whose chance of ever colliding is, by NumPy's own pessimistic estimate, about n² · 2⁻¹²⁸ for n streams (NumPy docs). Session 1.5 explained why you seed and why each experiment gets its own stream, so that editing one section never shifts the numbers in another. This session's lab starts the way every lab does, with seven streams for seven experiments:

lab = qn.Lab("3.5")  # seeds the random generators: every run gives the same numbers
(rng_se, rng_gbm, rng_dd, rng_bar, rng_boot, rng_vr, rng_mt) = lab.rng.spawn(7)

3.5.2 The Monte Carlo workflow & its standard error

Every simulation in this course, from a coin flip to a margin call, follows the same five steps. The fourth — attaching a standard error — is the one people skip.

Take a question with a known answer so you can watch the method work. This session's base asset is the one from Session 2.3: drift μ = 8% and volatility σ = 20% a year, no dividends, and a risk-free rate of zero throughout. What is the chance it loses money over a year? Under geometric Brownian motion (§3.5.3) the one-year log return is normal with mean ν = μ − σ²/2 = 6% and standard deviation σ, so the exact answer is Φ(−ν/σ) = 38.21%, where Φ is the standard normal cumulative distribution function. Pretend you don't know that. Step 1 is the model; step 2 draws n one-year log returns; step 3 records a 1 for each losing year and a 0 otherwise.

Step 4 is where the statistics come in. The estimate is the average of the n simulated values $Y_i$, and its Monte Carlo standard error is the standard error of Session 3.1 applied to your own simulated draws: their sample standard deviation $s_Y$ divided by √n. A probability is the expected value of an indicator (1 if the event happens, 0 if not), so its standard error has a closed form, and it is never more than 0.5/√n:

The interval narrows by √10 ≈ 3.2 with every tenfold increase in paths — 9.9 times narrower from 10,000 to 1,000,000. That is the √n rule of Session 3.1, and the economics of Monte Carlo: each extra decimal place costs a hundred times the computing. Notice also that 1 of the 5 estimates misses the truth by more than 2 SE. Nothing is wrong. "Estimate ± 2 SE" promises that about 95% of such intervals contain the truth, not all of them, and the lab checks the promise directly: it reruns the n = 1,000 experiment 4,000 times, and 96.2% of the intervals contain the true value, against 95.6% exactly for this discrete count (95.45% for a normal).

3.5.3 Random walks & GBM: simulate in log space

The workhorse model of a price is a random walk in its logarithm. Simulate it in log space and it is exact at any step size; simulate it the obvious way and it is not.

A random walk is what you get by adding up independent random steps: $X_k = X_{k-1} + \varepsilon_k$. Prices make a poor random walk — a $1 move means something different at $10 and at $1,000 — but log prices make a good one, because log returns add across time (Session 2.3). Geometric Brownian motion (GBM) is the model in which the log price is a random walk with normally distributed steps. Over a step of Δ years:

Here t and Δ are measured in years, μ is the expected (arithmetic) growth rate and σ the volatility. The log price drifts at ν = μ − σ²/2, lower than μ by the volatility drag of Session 2.3. Because the steps are independent normals, their sum over a horizon of τ years is normal with mean ντ and variance σ²τ, whatever the step. So this update — the exact scheme — is correct for any Δ: one 10-year step gives exactly the same distribution of the final price as 2,520 daily steps. Final prices are lognormal and every path stays above zero. Fig 3.5.3 shows three simulated decades of the base asset.

The obvious way to simulate a price is to draw normal simple returns with mean μΔ and standard deviation σ√Δ and compound them — the Euler scheme (the Euler–Maruyama scheme applied to the price itself):

It converges to GBM as Δ shrinks, but at any finite step it is an approximation, and with big steps a bad one: any draw with $Z \lt -(1 + \mu\Delta)/(\sigma\sqrt{\Delta})$ takes the price to zero or below, which no stock can do (fig 3.5.4).

How bad is the approximation? The lab simulates 40,000 ten-year paths of a volatile asset — 10% drift, 40% volatility, like a single growth stock. For a fair comparison it uses common random numbers: one set of daily shocks per path, summed into longer steps, so every scheme sees the same luck. It tracks two tail quantities with exact answers: the chance of losing more than half, and the 5th percentile of what $1 becomes.

The exact scheme lands on identical final prices whatever the step (the lab asserts it), within sampling error of theory; the first column's standard error is 0.21 percentage points. Euler with yearly steps puts the chance of losing half at 27.6%, 17 standard errors above theory, puts the 5th percentile at $0.035 instead of $0.152, and sends 3.1% of paths through zero. Measured on the same paths against the exact scheme, the excess shrinks with the step: +3.36 points with yearly steps, +0.91 quarterly, +0.24 monthly, and +0.05 ± 0.03 daily, lost in the noise. Euler gets the middle of the distribution nearly right at every step; its error lives in the tails, exactly where risk questions live.

The last row is the most common bug of all: simulating in log space but with drift μ instead of μ − σ²/2. Every path is shifted up by the same factor $e^{\sigma^2\tau/2}$ = 2.23, so the median final wealth comes out at $2.70 instead of $1.22, and the average at $6.07 instead of $2.72. The two schemes as the lab writes them:

def euler_terminal(dW: np.ndarray, step_days: int, mu: float, sigma: float) -> tuple[np.ndarray, np.ndarray]:
    """Naive Euler on PRICES with steps of `step_days` days: P ← P·(1 + μΔ + σΔW).

    dW holds daily Brownian increments (paths × days). Returns (terminal price from P0 = 1,
    ruined flag); a path whose price reaches zero or below is absorbed at 0 (it is bust).
    """
    n_paths = dW.shape[0]
    dW_k = dW.reshape(n_paths, -1, step_days).sum(axis=2)       # Brownian increment per step
    factor = 1.0 + mu * step_days * DT + sigma * dW_k
    ruined = (factor <= 0).any(axis=1)
    terminal = np.where(ruined, 0.0, np.prod(np.maximum(factor, 1e-300), axis=1))
    return terminal, ruined


def exact_terminal(dW: np.ndarray, step_days: int, mu: float, sigma: float) -> np.ndarray:
    """Exact GBM in log space with steps of `step_days` days: ln P += (μ − σ²/2)Δ + σΔW."""
    n_paths = dW.shape[0]
    dW_k = dW.reshape(n_paths, -1, step_days).sum(axis=2)
    return np.exp(((mu - 0.5 * sigma ** 2) * step_days * DT + sigma * dW_k).sum(axis=1))

3.5.4 Drawdowns & barriers: when the path matters

The questions a risk manager asks — how deep can it fall, will it hit the stop — are about the whole path. For the drawdown, simulation is not a check on the formula; it is the formula.

Start with the maximum drawdown, the worst fall from a running peak (Session 2.3). For GBM its distribution has no simple closed form — Magdon-Ismail et al. (2004) express it as an infinite series — so the lab simulates 6,000 ten-year daily paths of the base asset, whose Sharpe ratio is μ/σ = 0.4 at a zero risk-free rate, and computes each path's maximum drawdown with qn.returns.max_drawdown_paths.

The typical decade of a perfectly healthy asset contains a drawdown of −41.0% (batch-means SE 0.2 percentage points), and one decade in twenty goes beyond −63.7%. A drawdown worse than −30% happens in 86% of decades, worse than −50% in 25%. So a kill criterion (Session 1.3) must be set against the distribution, not a gut feeling: a rule that shuts a strategy at a −20% drawdown would have shut this asset down in 99.9% of decades.

Touching a barrier

Now a question with a formula: what is the chance the base asset falls 20% below its starting price at some point within a year — a stop-loss, or the margin call of Session 2.4? For a price watched continuously, the reflection principle gives the answer (fig 3.5.6). Take a driftless path that touches the barrier and ends above it. Reflect everything after the first touch in the barrier, and you get a path that ends below it — just as likely, because from the barrier a driftless random walk is as likely to go down as up. So the touched-and-ended-above paths pair off one-for-one with the ended-below paths, all of which touched, and P(touch) = 2 × P(end below).

With a drift the two halves of the pairing are no longer equally likely — the mirror image of an upward-drifting path drifts downward — and the second probability picks up the factor $e^{2\nu h/\sigma^2}$ (Shreve 2004, §3.7 has the derivation). With $h = \ln(B/P_0) \lt 0$ the log distance down to a barrier B:

This is qn.stats.p_touch. For a 20% fall within a year it gives 18.46%: 7.84% of ending below plus 10.62% of touching and recovering, the drift factor being 0.512. The lab measures the two terms separately at 7.86% and 10.51%.

A simulation, like your broker, only looks at the price at its time steps, and a dip that recovers before the next check goes unseen, so the discretely monitored probability is lower. Broadie, Glasserman & Kou (1997) showed that you can correct for this by moving the barrier further away: replace h by h − βσ√Δt, where Δt is the time between checks and β = −ζ(1/2)/√(2π) ≈ 0.5826 is their constant (ζ is the Riemann zeta function; this β is not a market beta). p_touch(..., dt=Δt) applies it. The lab simulates 200,000 one-year paths on a daily grid and looks at them daily, weekly, monthly and quarterly:

Checking at the daily close instead of continuously lowers the probability from 18.5% to 17.2%, far more than the sampling error, and the corrected formula matches daily checks within it. Weekly it is within a tenth of a point; with monthly and quarterly checks it undershoots, being an approximation built for frequent monitoring. And a drawdown is a barrier that moves up with every new peak, so it is always at least as likely as the same fall from the start: in the same simulation a −20% drawdown within the year happened on 37.4% of paths, against 17.2% for a 20% fall below the starting price.

3.5.5 The bootstrap: resampling history

When you do not trust any model, let the data be the model and resample the history you have. It is one of the most useful ideas in empirical finance — and it cannot see anything your history did not contain.

The bootstrap (Efron 1979) estimates how much a statistic would vary across samples you do not have by treating the sample you do have as the population. Draw a same-length history from it with replacement, recompute the statistic, repeat a few thousand times, and the spread of the results estimates its standard error — no formula needed, which is why Session 3.1 suggested it for fat-tailed returns.

Resampling days one at a time — the iid bootstrap — assumes the days are independent. Returns are not: big moves cluster in time (volatility clustering, Session 5.1). The block bootstrap (Künsch 1989) resamples runs of consecutive days instead, so calm and stormy stretches survive inside each block (fig 3.5.7). The lab uses the circular variant, which wraps around the end of the sample, with the block length chosen by the Politis–White rule (Politis & White 2004, corrected by Patton, Politis & White 2009) as implemented in the arch package: 93 trading days here. The stationary bootstrap of Politis & Romano (1994), with random block lengths, is the common alternative.

The lab's history is 10 years of daily returns from a GARCH model with Student-t shocks — synthetic, so the truth is known, and built to have the clustering and fat tails of real returns (Session 5.5). Clustering shows up as positive autocorrelation of squared returns: a big move today makes a big move tomorrow more likely.

Averaged over lags 1 to 10, the history's autocorrelation of squared returns is 0.124; iid resamples average 0.000 and block resamples 0.118. Why it matters: ask each bootstrap how precisely ten years pin down the annual volatility, and compare with the truth — the spread of the volatility estimate across 500 fresh ten-year histories drawn from the model.

The iid bootstrap understates the uncertainty 8-fold. The block bootstrap does better yet still falls short: volatility shocks here fade over months, so even 93-day blocks cut some of the memory, and one history can only show the storms it happened to contain.

What no bootstrap can fix

Two things, both visible in the lab. First, the bootstrap is centred on your sample, not on the truth: this history's volatility is 13.3%, the model's long-run volatility is 15.9%, and every resample inherits the gap. Second, it cannot invent what it has not seen. Keep only the first two years: their worst day is −2.23%, so no bootstrap year, iid or block, can contain a worse one — the lab finds a probability of exactly 0%. Under the model that generated the data, 89.1% of years contain a worse day.

3.5.6 Variance reduction: same precision, fewer paths

The √n rule makes precision expensive. Variance reduction beats it honestly: same model, same answer, a smaller standard error for the same number of random draws.

Variance reduction is any technique that lowers the variance of a Monte Carlo estimator at the same cost. Its payoff is the variance-reduction ratio (VRR): plain variance over new variance at an equal number of random draws — how many times more paths the plain method would need for the same standard error. The lab tests two classics on the expected payoff of a one-year call option on the base asset struck at its starting price, E[max(P_τ − K, 0)] with P₀ = K = $100. That expectation has a closed form, $13.11 — the Black–Scholes formula with the drift in place of the interest rate and no discounting, a black box until Session 10.3 derives it — which here serves only as the answer key.

Antithetic variates (Hammersley & Morton 1956) use every normal draw twice, as Z and −Z, and average the pair (fig 3.5.9): when f(Z) and f(−Z) are negatively correlated, their errors cancel. Control variates subtract a quantity whose mean you know exactly — here the final price $P_\tau$, with $\E[P_\tau] = P_0 e^{\mu\tau}$ — scaled to soak up the noise the two share:

Here ρ is the correlation between f(Z) and f(−Z), and $\rho_{YC}$ the correlation between the payoff Y and the control C. Estimating b from the same paths adds a bias of order 1/n, negligible at these sizes.

The control variate turns 200,000 draws into the precision of about 8.7 times as many, for the price of one subtraction — and Boyle's 1977 paper already used both tricks. The ratios match their formulas, so one correlation from a small pilot run tells you in advance whether a trick will pay.

3.5.7 Many checks, one lab: the multiple-testing trap

A lab that asserts dozens of Monte Carlo results against theory is running dozens of statistical tests. Set the bar at 2 or 3 standard errors and a perfectly correct lab will fail by luck.

Each assert is a test. If the code is right, an estimate's error measured in standard errors is roughly standard normal, so a check at c standard errors fails by pure luck with probability 2Φ(−c), the test's size. A lab with k independent checks fails at least one of them with probability

This is the family-wise error rate of multiple testing, which Session 5.2 treats properly. This course learned it the hard way: in review, one lab with about 27 asserts at 3 SE failed on 8% of alternative seeds (Quant Notebook authoring guide, §12a). The formula predicts 7.0% — close, for checks that are not exactly independent or exactly normal. The lab confirms the formula by simulating 200,000 imaginary labs with perfect code and 27 checks each:

Hence the rule every lab in this course follows: assert Monte Carlo results against theory at 4 standard errors. This lab routes all 22 of its statistical checks through one helper that counts them; on the committed seed the largest miss was 2.61 standard errors. Two refinements keep it honest. A looser threshold catches only bigger bugs, so buy sensitivity with more paths, not a tighter c. And a tolerance wider than 4 SE is allowed only for a known, written-down bias in the theory, like the 0.15-point allowance the lab gives the Broadie–Glasserman–Kou correction at weekly checks (§3.5.4):

def check_mc(name: str, estimate: float, truth: float, se: float, n_se: float = N_SE,
             allowance: float = 0.0) -> float:
    """Assert that a Monte Carlo estimate lies within n_se standard errors of its theoretical value.

    `allowance` widens the tolerance by a known approximation error of the THEORY (never used to
    rescue noise; every use carries a comment). Returns the z-score and remembers it, so the lab
    can count how many statistical checks it makes.
    """
    z = (estimate - truth) / se
    MC_CHECKS.append((name, float(z)))
    assert abs(estimate - truth) < n_se * se + allowance, (
        f"{name}: {estimate:.6g} vs theory {truth:.6g} is {z:+.2f} SE away")
    return float(z)

3.5.8 Try it

Open the companion lab, predict before you run, and watch the standard errors do what the formulas say.

3.5.9 Glossary · Check yourself · Go deeper

Check yourself

📚 Go deeper

1. Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering, ch. 3 "Generating Sample Paths" and ch. 4 "Variance Reduction Techniques". Springer. ISBN 978-0-387-00451-8
2. Efron, B. & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall. ISBN 978-0-412-04231-7
3. Efron, B. (1979). "Bootstrap Methods: Another Look at the Jackknife." Annals of Statistics 7(1), 1–26. doi:10.1214/aos/1176344552
4. Künsch, H. R. (1989). "The Jackknife and the Bootstrap for General Stationary Observations." Annals of Statistics 17(3), 1217–1241. doi:10.1214/aos/1176347265
5. Politis, D. N. & Romano, J. P. (1994). "The Stationary Bootstrap." Journal of the American Statistical Association 89(428), 1303–1313. doi:10.1080/01621459.1994.10476870
6. Broadie, M., Glasserman, P. & Kou, S. (1997). "A Continuity Correction for Discrete Barrier Options." Mathematical Finance 7(4), 325–349. doi:10.1111/1467-9965.00035
7. Shreve, S. E. (2004). Stochastic Calculus for Finance II: Continuous-Time Models, §3.7 "Reflection Principle". Springer. ISBN 978-0-387-40101-0

Sources for the numbers

1. Metropolis, N. & Ulam, S. (1949). "The Monte Carlo Method." Journal of the American Statistical Association 44(247), 335–341. doi:10.1080/01621459.1949.10483310
2. Metropolis, N. (1987). "The Beginning of the Monte Carlo Method." Los Alamos Science 15 (special issue), 125–130. https://mcnp.lanl.gov/pdf_files/Article_1987_LAS_Metropolis_125--130.pdf
3. Los Alamos National Laboratory (2023). "Hitting the Jackpot: The Birth of the Monte Carlo Method." Actinide Research Quarterly. https://www.lanl.gov/media/publications/actinide-research-quarterly/1123-hitting-the-jackpot-the-birth-of-the-monte-carlo-method
4. Boyle, P. P. (1977). "Options: A Monte Carlo Approach." Journal of Financial Economics 4(3), 323–338. doi:10.1016/0304-405X(77)90005-8
5. Hammersley, J. M. & Morton, K. W. (1956). "A New Monte Carlo Technique: Antithetic Variates." Mathematical Proceedings of the Cambridge Philosophical Society 52(3), 449–475. doi:10.1017/S0305004100031455
6. Magdon-Ismail, M., Atiya, A. F., Pratap, A. & Abu-Mostafa, Y. S. (2004). "On the Maximum Drawdown of a Brownian Motion." Journal of Applied Probability 41(1), 147–161. doi:10.1239/jap/1077134674
7. Politis, D. N. & White, H. (2004). "Automatic Block-Length Selection for the Dependent Bootstrap." Econometric Reviews 23(1), 53–70; with Patton, A., Politis, D. N. & White, H. (2009). "Correction." Econometric Reviews 28(4), 372–375. doi:10.1081/ETC-120028836 · doi:10.1080/07474930802459016
8. NumPy developers (2026). "Parallel random number generation" (SeedSequence spawning; PCG64). NumPy documentation. https://numpy.org/doc/stable/reference/random/parallel.html
9. Quant Notebook (2026). Authoring guide, §12a "Lessons from review" (the 27-assert lab that failed on 8% of alternative seeds). https://github.com/georgevrs/quant-notebook/blob/main/AUTHORING.md