How wrong can the estimate be?
A Monte Carlo estimate is the mean of n random draws, so it is itself random. Its typical error is the standard error: the standard deviation of the draws divided by sqrt(n), the same formula as in the t-test. Two consequences matter. First, the error does not depend on how hard the problem looks, only on the spread of the quantity you average. Second, it shrinks slowly: to gain one more digit of accuracy you need 100 times more draws.
import math
import numpy as np
# Estimate E[exp(Z)] for a standard normal Z; the exact value is exp(0.5) = 1.6487
for n in (100, 10000, 1000000):
values = np.exp(np.random.default_rng(5).standard_normal(n))
print(n, round(values.mean(), 4), round(values.std(ddof=1) / math.sqrt(n), 5))
# 100 1.2167 0.14583
# 10000 1.6887 0.0217 the error falls roughly like 1 / sqrt(n)
# 1000000 1.6508 0.00216Error bars and the cost of precision
Report a simulated number with its confidence interval, mean -/+ 1.96 * standard error: the range that contains the true value about 95% of the time. If you need the answer within a tolerance, invert the standard error formula: n = (std / target_error)², rounded up. With a spread of about 2.2 and a target of 0.001 you need almost five million draws; for 0.01 only 47,000. Run a small pilot simulation first to learn std, then decide how many draws the real run needs.
values = np.exp(np.random.default_rng(5).standard_normal(10000))
mean, std = values.mean(), values.std(ddof=1)
std_error = std / math.sqrt(len(values))
print(round(mean - 1.96 * std_error, 4), round(mean + 1.96 * std_error, 4)) # 1.6462 1.7312 contains the exact 1.6487
print(math.ceil((std / 0.01) ** 2)) # 47111 draws for a standard error of 0.01
print(math.ceil((std / 0.001) ** 2)) # 4710946 10 times the precision, 100 times the drawsCheaper accuracy: antithetic variates
Besides drawing more, you can draw smarter. For a standard normal Z, the mirror image -Z is just as valid a draw. If f is monotonic, f(Z) and f(-Z) move in opposite directions, so their average is much less noisy than either one. The antithetic estimator averages (f(z) + f(-z)) / 2 over the draws. In the example the standard error falls from 0.0217 to 0.0119 for the same number of random numbers, and for an odd function such as 3z + 1 the noise vanishes completely. This is one of several variance-reduction tricks.
z = np.random.default_rng(5).standard_normal(10000)
plain = np.exp(z)
antithetic = (np.exp(z) + np.exp(-z)) / 2 # each draw is paired with its mirror image
print(round(plain.mean(), 4), round(plain.std(ddof=1) / 100, 5)) # 1.6887 0.0217
print(round(antithetic.mean(), 4), round(antithetic.std(ddof=1) / 100, 5)) # 1.6565 0.01192 half the noise, same 10,000 numbers
print(round(math.exp(0.5), 4)) # 1.6487 the exact value# The running mean of one long simulation settles down as n grows
values = np.exp(np.random.default_rng(5).standard_normal(100000))
running = np.cumsum(values) / np.arange(1, len(values) + 1)
print([round(float(running[k - 1]), 4) for k in (10, 100, 1000, 100000)]) # [1.5522, 1.2167, 1.634, 1.6502]- Standard error =
std / sqrt(n): 100 times more draws give 10 times less error; reportmean -/+ 1.96 * se. - Draws needed for a target error:
ceil((std / target) ** 2); run a small pilot to learnstdfirst. - Antithetic variates: average
(f(z) + f(-z)) / 2to cut the noise at no extra random numbers; always check convergence with another seed.
Next up: present value and discounting, the first step into fixed income: a euro tomorrow is worth less than a euro today.
Lab · what's graded
- Measure the noise of an estimate with
standard_error25 pts - Put error bars on the answer with
mc_confidence_interval25 pts - Budget the simulation with
samples_needed25 pts - Get the same accuracy for less with
antithetic_mean25 pts