Beating the square-root law

The scorecard, before anything else

Every technique here reports one number: the variance-reduction factor, measured as the naive variance divided by the technique's variance at equal computational budget. It reads directly as “how many times more naive samples would buy the same accuracy”. A factor of 4 halves your error bars for free, forever. A factor of 683 replaces a compute cluster with a laptop. And a factor below 1 means the technique backfired - which is possible, which happens, and which only the scorecard will catch.

Insisting on that number is the whole discipline of this part. Variance reduction is the one area of simulation where folklore runs thickest - people apply a technique because they read that it helps, never measure whether it helped here, and occasionally ship something worse than what they started with.

Antithetic variates: cancellation by construction

Part 1 left a promissory note: independence makes errors cancel eventually, but engineered negative correlation makes them cancel faster. Here is the redemption, and it costs one line. Simulate in pairs: for every draw Z, also evaluate the payoff at minus Z, and average the pair.

If the payoff is monotone - a call rises with Z - then the two halves of each pair move oppositely, their correlation goes negative, and the pair is worth more than two independent draws. A bonus falls out for free: the pair's average shock is exactly zero, so drift-type errors vanish identically.

Measured on the European call at equal budget: correlation between the halves of a pair, minus 0.31; variance-reduction factor, 1.5. Modest - and that is the honest lesson. Antithetic costs one line and reliably helps on monotone payoffs, but it is a 1.5-to-2 times tool, not a revolution. On a non-monotone payoff - a straddle, anything symmetric in Z - that correlation can turn positive and the factor drops below 1. Cheap, mildly useful, occasionally treacherous. Run the scorecard, always.

Control variates: the professional's favourite weapon

Now the revolution. Suppose you're estimating something and you can find a sidekick variable, computed on the same scenarios, with two properties: it is strongly correlated with your target, and its exact mean is already known in closed form. Then the sidekick's simulation error is visible to you - you can see how far its sample mean landed from the truth you know - and since the two err together, that visible error tells you which way your invisible one went.

The correction has mean zero, so the estimator stays unbiased for any beta - and the optimal beta is a regression slope, because a control variate simply is the regression of your target on its sidekick. Look hard at that shrinkage factor, because it is the whole economics of the technique: a correlation of 0.9 buys 5 times, 0.99 buys 50 times, 0.999 buys 500 times. The entire craft is finding a sidekick that is almost your problem, but solvable.

And Part 3 planted one. The arithmetic Asian option - no formula, our patient - has a near-twin: the geometric Asian, whose payoff uses the geometric average, is lognormal-friendly, and does have a closed form. Arithmetic and geometric averages of the same path are nearly the same number, so the correlation should be enormous. The full liturgy runs in three steps.

• Certify the control. Never bet on a formula you haven't tested. The geometric Asian's closed form says 4.1017; a plain simulation of the same payoff says 4.0898 give or take 0.0243 - inside the interval, formula certified. This step looks like bureaucracy and is the thing that catches a mis-transcribed variance term.
• Measure the friendship. The correlation between the arithmetic and geometric payoffs comes out at 0.99927, with an optimal beta of 1.053. Plug that correlation into the shrinkage formula and it predicts a factor near 700 before you have run a single line of the real estimator.
• Collect. Naive arithmetic-Asian price: 4.34602 give or take 0.02558. Control-variate estimate: 4.35852 give or take 0.00098. Measured variance-reduction factor: 683 times.

Run the scorecard

The same European call, the same computational budget, three estimators. The bars are each method's error bar - the thing you are actually buying. The sidekick here is deliberately mediocre: the terminal price itself, whose exact mean we know, correlates about 0.89 with the payoff. Watch what a merely-good friend buys, then remember the Asian option's near-twin scored 0.99927.

The formula predicts the factor from the correlation alone, before you run anything: one over one minus 0.89 squared is about 4.7 - and that is what the scorecard measures. This is why finding a better sidekick, not buying more compute, is where a professional spends the afternoon.

Importance sampling: hunting rare events

You need the probability that a standard normal exceeds 4.5 - a stand-in for every probability-of-ruin, default, or catastrophic-failure question there is. The truth, computable here precisely because we are practising on a solvable case, is about 3.4 in a million. Naive Monte Carlo with a million draws sees a handful of events, and the error bar comes out as large as the answer itself.

And the scaling is hopeless by design. The relative error of a frequency estimate goes like one over the square root of the expected number of hits, so getting ten thousand hits on a one-in-a-million event means ten billion draws. The disease isn't slow convergence; it is that almost every sample lands in the boring region and contributes nothing but certainty about what you already knew.

The cure is audacious: stop sampling from the true distribution. Sample instead from a proposal that puts its mass where the interesting events live - then correct for the lie, sample by sample, with a likelihood ratio weight.

One line of algebra - multiply and divide by the proposal - and it says: the average of f under the true law equals the average of f times w under the proposal law. Samples from the exciting proposal, each carrying an honest weight that says “I was two hundred times more likely under q than under p, so I count for a two-hundredth of a real sample”, give an exactly unbiased estimator. All the cleverness moves into choosing the proposal.

For normal tails the canonical proposal is a mean shift: sample from a normal centred on the interesting region instead of on zero. With the shift set to 4.5 and the same million samples that failed before, the estimate lands at 3.399 in a million give or take 15 in a billion - a measured variance-reduction factor around 68,000 times, and a naive equivalent of roughly seventy billion samples. This is not an incremental improvement. It is the difference between “cannot be done on this hardware” and “runs in two seconds”, and it is the standard machinery behind rare-event numbers you have seen quoted with a straight face: one-in-a-billion network failure rates, deep value-at-risk quantiles, particle-physics detection significances.

Tilt the dice - and watch the smoke detector

Two hundred thousand draws either way. Naive counts how many happened to land past the threshold; importance sampling shifts the proposal there and pays with weights. The exact answer is known, because we are practising where we can check. Then push the tilt too far and watch the technique fail the way it fails in production: silently, and optimistically.

Aim the proposal at the threshold and importance sampling is a miracle. Aim it far past, and the estimator rides on a hundred-odd samples while claiming a million's worth of authority - the interval stays narrow because it is computed from the samples you got, and the samples you got don't know what they're missing. Report the effective sample size next to every importance-sampling estimate the way you report a standard error next to every mean.

Quasi-Monte Carlo: attacking the exponent itself

Every technique so far accepted the square-root law and fought over the constant. Quasi-Monte Carlo picks the other fight. Look at a few hundred random uniform points on a square: clumps here, bare patches there. That is not a flaw in your generator - equal probability everywhere guarantees clumps, the same way fair coins produce runs.

Stratification fought clumping on a fixed grid. Quasi-Monte Carlo asks the audacious follow-up: why not abandon randomness altogether and place the points deterministically, so that every region contains almost exactly its fair share at every prefix length? Sequences achieving this are called low-discrepancy, and the workhorse is the Sobol sequence. Feed Sobol points into any estimator you would have fed uniforms - through the inverse transform, exactly Part 2's pipeline - and the error falls roughly like one over N instead of one over the square root of N. At a million points, that is the difference between three digits and six.

The measurements. On the European call, sweeping the budget and fitting log-log slopes: pseudo-random comes out at minus 0.49, reassuringly close to the theoretical minus 0.5, while scrambled Sobol comes out at minus 0.83 - not a better constant but a visibly steeper line, pulling away without limit. On the 52-dimensional Asian option - the stress test - Sobol's error is 4.1 times smaller, worth about 17 times fewer samples.

Simulation studies as a science

Part 2's power project simulated tidy trials: independent patients, one per row. Real biostatistics is hierarchical - patients nest inside clinics, students inside schools - and often the unit of randomisation is the cluster, because you cannot give half a clinic a new workflow. Patients who share a clinic share a clinic effect, and even a small intraclass correlation changes everything.

Stage the crime under known truth: six clinics per arm, fifty patients each, an intraclass correlation of 0.05, and no true effect whatsoever. The plain patient-level t-test on 300 against 300 patients - what a hurried analyst reaches for - produces a false-positive rate of 28.8% against a nominal 5%. Not inflated: catastrophic. One in three null studies comes out “significant”. Collapse to cluster means and t-test those six against six, and the false-positive rate is 5.0%, exact.

The mechanism is Part 1's law of large numbers colliding with Part 3's correlation: positively correlated observations are worth fewer than their count, so a standard error computed as if they were independent comes out too small - here by a factor that inflates every t-statistic by about 1.9 times. Same disease as the parallel-seed trap. Same disease as MCMC autocorrelation below. Three costumes, one crime: counting correlated things as if they were independent.

That pattern - simulate the null, take the quantile of the maximum - is worth filing under its general name, because it is how multiple-comparison corrections, scan statistics and permutation-based error control all work. One pattern, half the questions in modern applied statistics. You have already run it on this site from the other end, in the multiple-comparisons case, where the maximum runs across contexts rather than across looks.

The bridge to MCMC

Everything in this series assumed one quiet superpower: that you can draw independent samples from your distribution, because it is a named distribution your library stocks, or an inverse CDF you can build, or a model you simulate forward. Bayesian inference demolishes that assumption on contact. A posterior is the likelihood times the prior, up to a constant - you can evaluate the right-hand side at any parameter value in microseconds, but the normalising constant is an integral over the whole parameter space: unknown, and in high dimensions unknowable. So here is the wall: a distribution you can touch pointwise but cannot sample.

The 1953 idea - Metropolis and colleagues, at Los Alamos again, the same building as Part 1's opening story five years on - is to stop trying to leap into the distribution and instead wander through it, taking steps whose local rules need only density ratios, in which the unknown constant cancels. Propose a small random step. Compute the ratio of the unnormalised densities. If it goes uphill, accept. If it goes downhill, accept with probability equal to that ratio; otherwise stay - and record the current point again.

The ritual, one last time. Certify against a case with a known answer: a coin showing 7 heads in 10 tosses under a uniform prior has an exactly known posterior with mean two-thirds. The chain's histogram lies indistinguishably on the true density, and its posterior mean comes out at 0.6657 against a true 0.6667 - inside the error bar. Machinery certified; now it may be pointed at posteriors nobody can write down.

But look at the price. The chain delivers the right distribution with the wrong independence: consecutive samples sit near each other by construction, so they are positively autocorrelated - and Part 1 told you exactly what that costs, because impaired cancellation is impaired cancellation. The exchange rate has a formula: the chain behaves like a smaller number of independent draws, the effective sample size. Measured here: a 58,000-draw chain is worth about 13,500 independent ones - 23% efficiency - and the posterior mean's standard error must be quoted using that number, not the raw count. Quoting MCMC error bars with raw N is pseudo-replication: the same crime, for the third time, as the parallel-seed trap and the clustered t-test.

The production checklist

Both guides compressed into the pre-flight list for any simulation that will touch a decision.

The whole series

That's the guide. Part 1 builds the engine and its error bar, Part 2 the recipe and two real projects, Part 3 paths, options and correlation, and all four sit on the series hub.