Monte Carlo Simulation: Pricing Paths, Variance and Reproducible Randomness
Price options by simulating thousands of price paths, report every estimate with a standard error, cut that error with antithetic and control variates, and keep every run reproducible.
When there is no formula to call
The last two lessons priced a European call two ways: with a binomial tree and with the Black-Scholes formula. Both work because the payoff depends only on the price at expiry, and the price at expiry has a known distribution. Change the contract slightly and both tools stop fitting. An Asian option pays on the average price over the year, sampled once a month. The average of twelve lognormal prices has no simple distribution, so there is no closed form to look up, and a recombining tree does not track the running average without growing a whole extra dimension.
Monte Carlo simulation handles this by brute force with a statistical guarantee attached. You simulate many possible price paths, compute the payoff on each, and average. The law of large numbers says the average converges to the expected payoff; the central limit theorem says how far off it is likely to be after N paths. That second part is what separates a Monte Carlo price from a random number: every estimate comes with a standard error, and an estimate reported without one is incomplete.
This lesson builds a Monte Carlo pricer from the random numbers up, checks it against Black-Scholes where the answer is known, and only then uses it on the Asian option. It also covers the two most useful variance-reduction techniques and reproducible randomness across threads. All parameters are illustrative; nothing here is investment advice.
A price is an expectation, and an average estimates it
Recall the setting of Black-Scholes. Under the risk-neutral measure, the price at time T of a stock that starts at S0, with interest rate r and volatility sigma, is
S_T = S0 * exp((r - sigma^2 / 2) * T + sigma * sqrt(T) * Z), with Z a standard normal variate.
The - sigma^2 / 2 term (the Ito correction) makes the expected value of S_T equal to S0 * exp(r * T): the stock grows on average at the risk-free rate, which is what "risk-neutral" means. The price of a contract that pays f(S_T) at T is the discounted expectation exp(-r * T) * E[f(S_T)]. For a call with strike K, f(S) = max(S - K, 0).
This formula is exact: one Z gives one terminal price with exactly the right distribution, with no time-stepping error. The recipe for a European call is:
- Draw
Z_1, ..., Z_N, independent standard normals. - For each, compute the terminal price and the discounted payoff
Y_i = exp(-r * T) * max(S_T,i - K, 0). - The estimate is the sample mean
Ybar = (Y_1 + ... + Y_N) / N.
A hand example shows the shape of it. With S0 = 100, K = 100, r = 0.05, sigma = 0.2, T = 1, the exponent is 0.03 + 0.2 * Z. Four draws Z = 1.0, -0.5, 0.2, -1.5 give terminal prices of about 125.9, 93.2, 107.3 and 76.3. The payoffs are 25.9, 0, 7.3 and 0, their mean is 8.3, and discounting gives about 7.9. The true value is 10.4506. Four paths are far too few, and the question is how many are enough and how to tell.
Normal variates you can reproduce
Everything starts with standard normals. The obvious tool is std::normal_distribution. For a system that must reproduce a run, it has a problem: the C++ standard specifies the random number engines exactly, but not the distributions. std::mt19937 must produce the same sequence on every implementation (the standard even states that its 10000th output from a default seed is 4123659995), while std::normal_distribution is only required to produce normally distributed values, by whatever algorithm the library chooses. Nothing guarantees that two standard libraries turn the same engine output into the same normals: same seed, different compiler, possibly a different price.
So we generate normals ourselves from the course's shared Lcg (the multiplier and increment are the constants Knuth gives for MMIX) with the Box-Muller transform. If U1 and U2 are independent and uniform on (0, 1), then with R = sqrt(-2 * ln(U1)) and theta = 2 * pi * U2, the numbers R * cos(theta) and R * sin(theta) are two independent standard normals. The open interval matters: ln(0) is minus infinity, which is why uniform() adds 0.5 before dividing and never returns 0.
#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>
#include <random>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
// Box-Muller: two independent uniforms give two independent standard normals.
struct Normal {
Lcg rng;
double spare = 0;
bool has_spare = false;
double next() {
if (has_spare) { has_spare = false; return spare; }
double r = std::sqrt(-2.0 * std::log(rng.uniform()));
double theta = 2.0 * std::numbers::pi * rng.uniform();
spare = r * std::sin(theta);
has_spare = true;
return r * std::cos(theta);
}
};
int main() {
Normal z{Lcg{42}};
std::print("first four:");
for (int i = 0; i < 4; ++i) std::print(" {:.4f}", z.next());
std::println("");
const int n = 1'000'000;
double sum = 0, sum_sq = 0;
int beyond = 0; // |Z| > 1.96 should happen about 5% of the time
for (int i = 0; i < n; ++i) {
double x = z.next();
sum += x;
sum_sq += x * x;
beyond += std::abs(x) > 1.96;
}
double mean = sum / n;
std::println("mean {:.4f} (expect 0)", mean);
std::println("variance {:.4f} (expect 1)", sum_sq / n - mean * mean);
std::println("|Z|>1.96 {:.4f} (expect 0.0500)", double(beyond) / n);
// The engines are specified exactly; the standard even states this value.
std::mt19937 mt;
for (int i = 1; i < 10000; ++i) mt();
std::println("10000th mt19937 output: {}", mt());
}Press Run to compile and run this example.
first four: 0.1633 1.0506 -0.9081 -0.9719
mean -0.0007 (expect 0)
variance 1.0015 (expect 1)
|Z|>1.96 0.0501 (expect 0.0500)
10000th mt19937 output: 4123659995The sample mean, variance and the fraction beyond 1.96 all land within sampling noise of the targets. The variance, for example, has a standard error of about sqrt(2 / N) = 0.0014 here, so 1.0015 is one standard error away.
Two limits of this generator are worth knowing. uniform() has 31 bits of resolution, so the largest normal it can produce is about 6.7 in absolute value; that does not matter for prices but would for estimating very rare events. And Box-Muller calls log, sqrt, cos and sin, which can differ in the last bit between math libraries. The programs in this lesson print at most four decimals so that such differences cannot show.
Pricing a European call and measuring the error
The sample mean Ybar is itself a random variable: run with another seed and it changes. If the Y_i are independent with variance sigma_Y^2, then
Var(Ybar) = sigma_Y^2 / N, so the standard error of the estimate is SE = sigma_Y / sqrt(N).
We do not know sigma_Y, so we estimate it with the sample standard deviation s (dividing by N - 1) and report SE = s / sqrt(N). By the central limit theorem Ybar is approximately normal for large N, so Ybar +- 1.96 * SE is an approximate 95% confidence interval: in repeated runs, about 95% of such intervals contain the true price.
The program below prices the call with N from 1,000 to 256,000 and compares with the Black-Scholes value. The Welford struct is the online mean and variance from Returns, Volatility and Risk Metrics. The last column, z = (estimate - exact) / SE, says how many standard errors the estimate is from the truth.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
struct Normal { // Box-Muller on the Lcg
Lcg rng;
double spare = 0;
bool has_spare = false;
double next() {
if (has_spare) { has_spare = false; return spare; }
double r = std::sqrt(-2.0 * std::log(rng.uniform()));
double theta = 2.0 * std::numbers::pi * rng.uniform();
spare = r * std::sin(theta);
has_spare = true;
return r * std::cos(theta);
}
};
struct Welford { // running mean and variance, as in lesson 7
long long n = 0;
double mean = 0.0, m2 = 0.0;
void add(double x) {
++n;
double delta = x - mean;
mean += delta / static_cast<double>(n);
m2 += delta * (x - mean);
}
double variance() const { return n > 1 ? m2 / static_cast<double>(n - 1) : 0.0; }
// Standard error of the mean: the spread of `mean` itself, not of one sample.
double std_error() const { return std::sqrt(variance() / static_cast<double>(n)); }
};
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double bs_call(double s, double k, double r, double sigma, double t) {
double sd = sigma * std::sqrt(t);
double d1 = (std::log(s / k) + (r + 0.5 * sigma * sigma) * t) / sd;
double d2 = d1 - sd;
return s * norm_cdf(d1) - k * std::exp(-r * t) * norm_cdf(d2);
}
int main() {
const double s0 = 100, k = 100, r = 0.05, sigma = 0.2, t = 1;
const double drift = (r - 0.5 * sigma * sigma) * t; // risk-neutral, with Ito term
const double vol = sigma * std::sqrt(t);
const double discount = std::exp(-r * t);
const double exact = bs_call(s0, k, r, sigma, t);
std::println("Black-Scholes {:.4f}", exact);
std::println("{:>7} {:>9} {:>8} {:>21} {:>6}", "N", "estimate", "std err",
"95% interval", "z");
for (int n = 1000; n <= 256'000; n *= 4) {
Normal z{Lcg{2024}};
Welford payoff;
for (int i = 0; i < n; ++i) {
double st = s0 * std::exp(drift + vol * z.next()); // exact, one step
payoff.add(discount * std::max(st - k, 0.0));
}
double est = payoff.mean, se = payoff.std_error();
std::println("{:>7} {:>9.4f} {:>8.4f} [{:>7.4f}, {:>7.4f}] {:>6.2f}", n, est,
se, est - 1.96 * se, est + 1.96 * se, (est - exact) / se);
}
}Press Run to compile and run this example.
Black-Scholes 10.4506
N estimate std err 95% interval z
1000 10.6089 0.4551 [ 9.7170, 11.5009] 0.35
4000 10.3555 0.2333 [ 9.8983, 10.8127] -0.41
16000 10.3867 0.1152 [10.1609, 10.6125] -0.55
64000 10.5601 0.0586 [10.4452, 10.6749] 1.87
256000 10.5034 0.0292 [10.4462, 10.5606] 1.81Read the std err column downward: each time N grows by a factor of 4, the standard error halves. That is the 1 / sqrt(N) law, and it is the central fact of Monte Carlo economics. One more correct digit costs 100 times the paths. Every z is below 3 in absolute value, so every estimate is consistent with Black-Scholes. Notice the run with 64,000 paths: it is 1.87 standard errors high, which happens to about 6% of honest runs. A test that demanded |z| < 1.5 would fail here on a correct program.
The cost is O(N) time and O(1) memory: one normal and one exponential per path, and Welford keeps three numbers.
A statistical check can fail by chance
"Within 3 standard errors of Black-Scholes" is a good validation test for a Monte Carlo pricer. It catches a missing discount factor, a wrong drift or a sign error, all of which push z to 10 or more. But it is a statistical statement: for a correct pricer, |z| > 3 still happens with probability about 0.27%, and the 95% interval misses about one run in twenty.
The next program prices the same call 400 times with 10,000 paths, each run with its own seed, and counts how often the interval contains the true value.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
struct Normal { // Box-Muller on the Lcg
Lcg rng;
double spare = 0;
bool has_spare = false;
double next() {
if (has_spare) { has_spare = false; return spare; }
double r = std::sqrt(-2.0 * std::log(rng.uniform()));
double theta = 2.0 * std::numbers::pi * rng.uniform();
spare = r * std::sin(theta);
has_spare = true;
return r * std::cos(theta);
}
};
struct Welford { // running mean and variance, as in lesson 7
long long n = 0;
double mean = 0.0, m2 = 0.0;
void add(double x) {
++n;
double delta = x - mean;
mean += delta / static_cast<double>(n);
m2 += delta * (x - mean);
}
double variance() const { return n > 1 ? m2 / static_cast<double>(n - 1) : 0.0; }
// Standard error of the mean: the spread of `mean` itself, not of one sample.
double std_error() const { return std::sqrt(variance() / static_cast<double>(n)); }
};
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double bs_call(double s, double k, double r, double sigma, double t) {
double sd = sigma * std::sqrt(t);
double d1 = (std::log(s / k) + (r + 0.5 * sigma * sigma) * t) / sd;
double d2 = d1 - sd;
return s * norm_cdf(d1) - k * std::exp(-r * t) * norm_cdf(d2);
}
// SplitMix64 (Steele, Lea and Flood): turns 0, 1, 2, ... into well-mixed seeds.
std::uint64_t splitmix64(std::uint64_t x) {
x += 0x9E3779B97F4A7C15ULL;
x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ULL;
x = (x ^ (x >> 27)) * 0x94D049BB133111EBULL;
return x ^ (x >> 31);
}
int main() {
const double s0 = 100, k = 100, r = 0.05, sigma = 0.2, t = 1;
const double drift = (r - 0.5 * sigma * sigma) * t, vol = sigma * std::sqrt(t);
const double discount = std::exp(-r * t);
const double exact = bs_call(s0, k, r, sigma, t);
const int runs = 400, n = 10'000;
int inside95 = 0, beyond3 = 0;
double worst = 0;
for (int run = 0; run < runs; ++run) {
Normal z{Lcg{splitmix64(run)}};
Welford payoff;
for (int i = 0; i < n; ++i) {
double st = s0 * std::exp(drift + vol * z.next());
payoff.add(discount * std::max(st - k, 0.0));
}
double zscore = (payoff.mean - exact) / payoff.std_error();
inside95 += std::abs(zscore) <= 1.96;
beyond3 += std::abs(zscore) > 3;
worst = std::max(worst, std::abs(zscore));
}
std::println("runs of {} paths: {}", n, runs);
std::println("95% interval holds BS: {} ({:.1f}%)", inside95,
100.0 * inside95 / runs);
std::println("more than 3 std errors off: {} (expected about {:.1f})", beyond3,
runs * 0.0027);
std::println("largest |z|: {:.2f}", worst);
}Press Run to compile and run this example.
runs of 10000 paths: 400
95% interval holds BS: 385 (96.2%)
more than 3 std errors off: 2 (expected about 1.1)
largest |z|: 3.27About 95% coverage, as promised. Two of the 400 correct runs are more than 3 standard errors off, against about one expected; the largest is 3.27. If one of those seeds had been the one in your unit test, the test would have failed on correct code.
A seeded validation test either always passes or always fails, so choose the seed once. If it fails at |z| = 3.3, rerun with more paths before changing anything. A real bug makes z grow like sqrt(N), because its bias stays while the standard error shrinks; bad luck does not.
splitmix64 is the output function of the SplitMix generator (Steele, Lea and Flood, 2014): a bit-mixer that maps run numbers 0, 1, 2, ... to seeds that look unrelated.
Antithetic variates and control variates
Variance reduction lowers sigma_Y itself, so the same N gives a smaller error. Glasserman's "Monte Carlo Methods in Financial Engineering" treats both of the following techniques at length in chapter 4.
Antithetic variates. If Z is standard normal, so is -Z. For each draw, evaluate the payoff at both, Y = f(Z) and Y' = f(-Z), and use their average as one sample. The variance of the pair average is (Var(Y) + Cov(Y, Y')) / 2. When the payoff is monotone in Z, as a call is, Y and Y' move in opposite directions, the covariance is negative, and the pair average varies less than two independent payoffs would. The standard error must come from the N / 2 pair averages, not from N payoffs treated as independent (they are not).
Control variates. Find a quantity X, computed on the same paths, that is correlated with the payoff and whose expectation you know exactly. Then for any constant b,
Y(b) = Y - b * (X - E[X]) has the same mean as Y, and variance Var(Y) - 2 * b * Cov(Y, X) + b^2 * Var(X).
Setting the derivative in b to zero gives the best coefficient b* = Cov(Y, X) / Var(X), and at b* the variance is Var(Y) * (1 - rho^2), where rho is the correlation of Y and X. A correlation of 0.9 cuts the variance by a factor of about 5; 0.99 by about 50. For the call, a natural control is the discounted terminal price X = exp(-r * T) * S_T, whose risk-neutral expectation is exactly S0. In practice b* is estimated from the same sample. That introduces a bias of order 1 / N, which is negligible next to the 1 / sqrt(N) standard error for any realistic N (Glasserman, section 4.1).
The program compares the three estimators with the same number of payoff evaluations. The gain column is (SE_plain / SE)^2: how many times more plain paths you would need to match that standard error.
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <numbers>
#include <print>
#include <vector>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
double uniform() { return (next() + 0.5) / 2147483648.0; } // in (0, 1)
};
struct Normal { // Box-Muller on the Lcg
Lcg rng;
double spare = 0;
bool has_spare = false;
double next() {
if (has_spare) { has_spare = false; return spare; }
double r = std::sqrt(-2.0 * std::log(rng.uniform()));
double theta = 2.0 * std::numbers::pi * rng.uniform();
spare = r * std::sin(theta);
has_spare = true;
return r * std::cos(theta);
}
};
struct Welford { // running mean and variance, as in lesson 7
long long n = 0;
double mean = 0.0, m2 = 0.0;
void add(double x) {
++n;
double delta = x - mean;
mean += delta / static_cast<double>(n);
m2 += delta * (x - mean);
}
double variance() const { return n > 1 ? m2 / static_cast<double>(n - 1) : 0.0; }
// Standard error of the mean: the spread of `mean` itself, not of one sample.
double std_error() const { return std::sqrt(variance() / static_cast<double>(n)); }
};
double norm_cdf(double x) { return 0.5 * std::erfc(-x / std::numbers::sqrt2); }
double bs_call(double s, double k, double r, double sigma, double t) {
double sd = sigma * std::sqrt(t);
double d1 = (std::log(s / k) + (r + 0.5 * sigma * sigma) * t) / sd;
double d2 = d1 - sd;
return s * norm_cdf(d1) - k * std::exp(-r * t) * norm_cdf(d2);
}
const double s0 = 100, k = 100, r = 0.05, sigma = 0.2, t = 1;
const double drift = (r - 0.5 * sigma * sigma) * t, vol = sigma * std::sqrt(t);
const double discount = std::exp(-r * t);
double terminal(double z) { return s0 * std::exp(drift + vol * z); }
double payoff(double z) { return discount * std::max(terminal(z) - k, 0.0); }
void report(const char* name, double est, double se, double se_plain, double exact) {
std::println("{:<12} {:>9.4f} {:>8.4f} {:>7.2f} {:>6.2f}", name, est, se,
se_plain * se_plain / (se * se), (est - exact) / se);
}
int main() {
const int n = 10'000; // every method evaluates the payoff n times
const double exact = bs_call(s0, k, r, sigma, t);
std::println("{:<12} {:>9} {:>8} {:>7} {:>6}", "method", "estimate", "std err",
"gain", "z");
Normal z1{Lcg{7}};
Welford plain;
for (int i = 0; i < n; ++i) plain.add(payoff(z1.next()));
const double se_plain = plain.std_error();
report("plain", plain.mean, se_plain, se_plain, exact);
// Antithetic: each pair (z, -z) is ONE sample, its average; n/2 samples.
Normal z2{Lcg{7}};
Welford pairs;
for (int i = 0; i < n / 2; ++i) {
double z = z2.next();
pairs.add(0.5 * (payoff(z) + payoff(-z)));
}
report("antithetic", pairs.mean, pairs.std_error(), se_plain, exact);
// Control variate: X = discounted S_T has known mean s0 under the
// risk-neutral measure. Estimate b = Cov(Y, X) / Var(X) from the sample.
Normal z3{Lcg{7}};
std::vector<double> ys(n), xs(n);
for (int i = 0; i < n; ++i) {
double z = z3.next();
ys[i] = payoff(z);
xs[i] = discount * terminal(z);
}
Welford ym, xm;
for (int i = 0; i < n; ++i) { ym.add(ys[i]); xm.add(xs[i]); }
double sxy = 0, sxx = 0;
for (int i = 0; i < n; ++i) {
sxy += (ys[i] - ym.mean) * (xs[i] - xm.mean);
sxx += (xs[i] - xm.mean) * (xs[i] - xm.mean);
}
const double b = sxy / sxx;
Welford controlled;
for (int i = 0; i < n; ++i) controlled.add(ys[i] - b * (xs[i] - s0));
report("control", controlled.mean, controlled.std_error(), se_plain, exact);
std::println("b = {:.4f}, correlation(Y, X) = {:.4f}", b,
sxy / std::sqrt(sxx * ym.variance() * (n - 1)));
std::println("Black-Scholes {:.4f}", exact);
}Press Run to compile and run this example.
method estimate std err gain z
plain 10.3584 0.1478 1.00 -0.62
antithetic 10.5991 0.1046 1.99 1.42
control 10.4516 0.0562 6.91 0.02
b = 0.6753, correlation(Y, X) = 0.9248
Black-Scholes 10.4506Antithetic sampling doubles the efficiency here; the control variate gives a factor of about 7. The measured gain matches the theory: 1 / (1 - 0.9248^2) = 6.9. The two methods also combine, by using antithetic pairs and a control on the pair averages.
Both have limits. Antithetic variates can increase variance for a payoff that is not monotone in Z, such as a straddle, which pays |S_T - K|. A control variate is only as good as its correlation: deep out of the money, the call payoff is mostly zero, its correlation with S_T is weak, and the gain shrinks.