Monte Carlo simulation, explained from first principles

The engine, in one sentence

That's the whole subject. Everything else is this sentence plus creativity about what you sample and how. It helps to hear it in three dialects, because different people find different versions obvious - and noticing they are the same thing is itself the first real lesson.

• The gambler. Want a probability? Run the experiment many times and count the fraction of times it happens. This estimates a probability as a frequency.
• The accountant. Want an expected value? Generate many scenarios, compute the quantity in each, average them. Simulate ten thousand futures for a stock price, take the payoff in each, average. This estimates an expectation as a sample mean.
• The mathematician. Want an integral with no closed form? That integral is by definition an expected value, so draw samples and average. This estimates an integral as a sample mean.

Pause on the third one, because it hides a small explosion: integration - a deterministic operation from calculus - can be performed by throwing dice. And for a large class of problems it is the best known way to do it.

The three are one. A probability is just the expected value of an indicator - the fraction of wins is the average of a pile of ones and zeros - and an expectation is an integral. So there is only one engine, seen at three levels of formality:

Every Monte Carlo method you will ever meet is this formula, plus choices about what X is, what f is, and how to draw the samples efficiently.

Where it came from: a sick physicist plays cards

In 1946 the mathematician Stanislaw Ulam was recovering from an illness and playing solitaire. Being a mathematician he asked: what is the probability a randomly dealt game is winnable? He tried combinatorics, discovered the number of game states is astronomical, and then had a lazier and much better idea - just play a hundred games and count the wins. No enumeration. Repetition and counting.

Ulam was at Los Alamos, and John von Neumann saw the same trick could crack a problem they were genuinely stuck on: neutron diffusion. Each individual collision was governed by known probabilities, but the aggregate question - does the chain reaction sustain itself? - had no solvable equations for real geometries. Their move was to simulate individual neutron life stories, using random numbers to decide each scatter and each absorption, then tally thousands of them. The statistics of the simulated population approximate the physics of the real one.

The project needed a code name. Ulam's uncle liked to borrow money from relatives to gamble at the casino in Monte Carlo. It's a good name: the method is, in a precise sense, industrialised gambling.

What Monte Carlo is not

Because the term is broad it collides with its neighbours, and this is where most of the confusion in that Reddit thread actually lives. Fences:

• Monte Carlo simulation - you know the probability model, you can draw independent samples from it directly, and you use them to estimate something. You are simulating from a model you chose. That's this lesson.
• The bootstrap - you have no model, you have data. You resample your own dataset with replacement to estimate a statistic's uncertainty. Same repeat-and-average spirit, but the randomness comes from your data rather than a model.
• MCMC - you have a distribution you cannot sample from directly, typically a Bayesian posterior known only up to a constant. So you build a random walk whose long-run behaviour mimics it. The samples are not independent, which changes both the maths and the diagnostics.
• Monte Carlo algorithms - in the computer-science sense (as opposed to Las Vegas algorithms): randomised algorithms returning a probably-correct answer in bounded time. Related in spirit, different community.

Your first simulation: estimating pi by throwing darts

Everyone's first Monte Carlo estimates pi, and it's also the single most common Monte Carlo interview question - so let's do it as an anatomy lesson rather than a party trick. Every serious simulation you build later has this same skeleton.

Take the unit square. Inscribe a quarter circle of radius 1 in it. The square has area 1; the quarter circle has area pi/4, which is about 0.7854. Now throw a dart that lands uniformly at random in the square, so that any region's chance of catching it is proportional to its area. Then:

Read that right to left and you'll see the trick: pi, a deterministic constant, has been expressed as a probability - and probabilities are things we know how to estimate by counting. Throw N darts, take the fraction that land inside, multiply by 4.

This manoeuvre - rewrite the thing I want as an expected value of something I can simulate - is the fundamental move of the entire field. Here it looks like a curiosity. In part 2 the same move is “rewrite the value of a business as the average NPV over simulated futures”, and in part 3 “rewrite an option price as a discounted average payoff”. Same move every time.

Throw the darts

Blue darts landed inside the quarter circle, red outside. The estimate is four times the blue fraction. Drag N and watch the estimate walk toward pi - slowly, and never in a straight line.

It is perfectly normal for the estimate at 200,000 darts to be momentarily worse than it was at 150,000. Convergence in probability is a drunkard weaving toward home, not a train on rails.

Why it works: the law of large numbers

The dartboard works because of a theorem you already believe: as you average more and more independent draws of the same random quantity, the average converges to that quantity's true expected value. Flip a fair coin ten times and you might see seven heads; flip it ten million times and the fraction will be flatteringly close to a half.

The intuition worth carrying is cancellation. Each draw misses the true value by some amount, sometimes high, sometimes low. Those misses are independent, so they have no reason to agree with each other, and when you add them up they partly cancel. The sum of N draws grows like N, but the accumulated error grows only like the square root of N - so the error per draw, which is what an average is, shrinks. Averaging doesn't remove randomness; it lets randomness eat itself.

How wrong are you? The square-root law

The law of large numbers promises you arrive eventually. It says nothing about when, and “eventually” does not pay for compute. The question every practitioner actually needs is: my simulation gave me a number - if I ran it again, how different would the number be?

Start by treating the estimate as what it is: not a number, but a random variable with its own distribution, its own centre and its own spread. The central limit theorem tells you the shape of that distribution - bell-shaped, centred on the truth - and, more usefully, its width:

That is the standard error: the typical distance between your estimate and the truth. Two consequences run the whole field.

• Accuracy is bought in bulk. The square root means halving your error costs four times the runs, and each extra correct digit costs a hundred times more. This is why brute force runs out of road, and why the professional toolkit is mostly about shrinking the numerator instead of buying the denominator.
• It does not care about dimension. Nothing in that formula mentions how many variables you have. Deterministic grid methods die of combinatorial explosion as dimensions grow; Monte Carlo keeps plodding at the same rate whether the problem has 2 uncertain inputs or 200. We don't use Monte Carlo because it's fast. We use it because it's the only thing left standing when problems get big.

Run it again: the estimate is a random variable

This is the picture beginners never think to ask for. Each bar is one complete run of the dart experiment at your chosen N. The spread of the histogram is what “how wrong am I?” means - and the theoretical square-root prediction is printed beside what the simulation actually did.

Two things to notice. The measured spread tracks the formula's prediction - that is the central limit theorem being right in front of you. And quadrupling the darts per run only halves the width: drag the slider from 1,000 to 4,000 and watch the histogram tighten by half, not by four.

Where random numbers come from

A confession: there is no randomness in your computer. Every number your simulation calls random came out of a deterministic algorithm - a pseudo-random generator - which, given a starting point, produces the identical stream forever. The starting point is the seed.

This sounds like a compromise. It is a superpower. Fixing the seed makes a simulation reproduce exactly: your results become auditable, a colleague can rerun your analysis and get your number to the last digit, and a bug becomes findable because the failure happens the same way twice. A simulation that can't be reproduced can't really be reviewed - which is why every number quoted in this series comes from a seeded run you can repeat.

What's next in this series

• Part 2 - the recipe, and two real projects. Turning uniform random numbers into any distribution you want, the five-step recipe every simulation follows, then a valuation DCF rebuilt as a distribution of value, and statistical power derived by simulation instead of a formula.
• Part 3 - paths and options. What changes when time enters: random walks, geometric Brownian motion, and the option-pricing workflow certified against Black-Scholes. Correlated inputs via Cholesky, and the dependence modelling error that helped detonate 2008.
• Part 4 - beating the square-root law. Antithetic and control variates, importance sampling for rare events, quasi-Monte Carlo, simulation studies as a science, and the bridge to MCMC.

Already run a Monte Carlo on this site without noticing? The look-everywhere simulator in our multiple-comparisons case does exactly this: simulate the null many times, take the quantile of the maximum. That pattern is how multiple-comparison corrections and permutation tests work, and Part 4 takes it professional.

Sources

• r/AskStatistics - the question this lesson answers.
• Metropolis, N. and Ulam, S. (1949). The Monte Carlo Method. Journal of the American Statistical Association - the paper that named and framed the method.
• Eckhardt, R. (1987). Stan Ulam, John von Neumann, and the Monte Carlo Method. Los Alamos Science - the solitaire story, from the people who were there.