Monte Carlo Simulation

The core idea

Nobody knows how an investment will develop over the next years. A Monte Carlo simulation plays through thousands of possible futures to show what you can expect and how wide the range of outcomes is. First you make assumptions, for example that the investment earns 8 % a year on average and swings by about 20 % a year. Then the computer builds one path, year by year. For each year it draws a random return around the 8 %: one year might bring plus 27 %, the next minus 9 %, the next plus 4 %. Multiplying these yearly returns one after the other gives where this path ends after 30 years. The computer repeats this 10,000 times with new random numbers, so you get 10,000 different futures under the same assumptions. From all paths together you read off what is typical, how bad a bad case is and how likely a loss is. The answer can only be as good as the assumptions behind it.

The idea in one example

Put $10,000 into an investment for 30 years. How much might you end with, and how likely is a loss? A Monte Carlo simulation answers this with the method from the core idea. It got its name in the 1940s from the casino of Monte Carlo, a place that lives on chance.

Figure 1 shows the result for 10,000 paths. Most of them end well below the average, because a few very lucky paths pull it up.

The numbers in Figure 1, such as the typical outcome of $60,814 and the 5 % chance of a loss, are estimates. They come from one run of 10,000 random paths. A second run draws new random returns for every year, so every path and every number comes out a little different. The exact typical outcome of this model is $60,496, so this run is about $300 off. The rest of this module shows how large that error is and how to repeat a run exactly. It also checks the simulation against exact answers and shows a trick that makes the results more precise with the same number of paths.

Many possible futures

Each run of the model is called a path. The model is the rule behind every path: here, a random return for each year around the assumed average.

Read off what matters

Sort the 10,000 end values. The middle one is the median, the typical outcome. The value with 10 % of the paths below it is the 10th percentile, a bad case, and the value with 90 % below it is the 90th percentile, a good case. The share of paths that end below the start is the chance of a loss. Each of these is a count or an average over the paths.

The average and the typical future

Gains and losses multiply (module 1), so a few very lucky paths end extremely high, while most paths gather lower down. The average is pulled up by those few. In the simulation of Figure 1, 70.97 % of the paths end below the average of $111,469, the same effect as the volatility drag of module 10.

More paths, more precision

With 100 paths the median could have come out anywhere from about $42,000 to $87,000, depending on the luck of the draws, as Figure 7 shows. With 10,000 paths it is far more precise. The error shrinks slowly, though: a hundred times as many paths make the numbers ten times as precise, the same rule as for the spread of a random walk (module 6).

A seed makes a run repeatable

A computer makes its random numbers with a fixed rule that starts from one number, the seed. The same seed gives exactly the same numbers and so the same paths, which lets anyone repeat and check a result. The simulation of Figure 1 uses the seed 17, fixed before the first run.

Only as good as the model

The 8 % and the 20 % are assumptions, and so is the bell shape of the yearly returns. More paths make the answer more precise for this model, and they cannot tell whether the model fits the real world. Module 17 shows how real returns break the bell shape.

The double arrow ⟹ reads “and therefore”, the sign √ stands for the square root and the sign ≈ means “about equal”.

One simulated path is a guess with a large error. The average of many paths is much more precise, and its error shrinks with the square root of the number of paths.

import numpy as np                               # numpy draws the random numbers

start, years, paths = 10_000, 30, 10_000         # $10,000, 30 years, 10,000 paths
mu, sigma = 0.08, 0.20                           # average return and volatility a year
rng = np.random.default_rng(17)                  # fixed seed: the same paths every run
g = mu - sigma ** 2 / 2                          # the growth rate of the log returns
log_ret = rng.normal(g, sigma, (paths, years))   # a random log return per path and year
final = start * np.exp(log_ret.sum(axis=1))      # what $10,000 grows into on each path
avg = final.mean()                               # the Monte Carlo estimate, an average
s = final.std(ddof=1)                            # the spread of one path result
se = s / paths ** 0.5                            # the standard error, s over root of M
print(f"spread of one path: ${s:,.0f}")          # print the spread s
print(f"average ${avg:,.0f} +/- ${se:,.0f}")     # the estimate and its standard error

Output
spread of one path: $167,084
average $111,469 +/- $1,671

Going deeper: the math step by step

For readers who want the full math: every step derived, with the full code behind the page at the end.

1. One path of the model

Call the start W₀ = $10,000 and the value after t years W_t, with T = 30 years in all. Each year the log of the value moves by a log return Y_t (module 1), drawn at random from a bell (module 7) with the mean g and the standard deviation σ (sigma), the volatility. The model assumes that the years are independent of each other, and that g and σ stay the same and are known. Log returns add up over time (module 1), and exp(y) is e to the power y. The sign Σ means “add all of these”, here for t from 1 to T.

The model of Figure 1 assumes an average return μ (mu) of 8 % a year, compounded continuously as in module 10, and σ = 20 %. Here μ is defined so that e^μ is the expected growth factor of a year, which is what “compounded continuously” means. With this μ the rule of module 10 holds exactly, so the log returns grow at g = μ − σ²/2 = 0.08 − 0.02 = 0.06 on average. As simple returns, a year then has the average growth factor e^0.08 = 1.0833, an average simple return of 8.33 %, and the median simple return e^0.06 − 1 = 6.18 %. A sum of independent bells is again a bell, and its means and variances add up (module 6). So ln(W_T/W₀) follows a bell with the mean gT = 1.8 and the standard deviation σ√T = 0.2 × √30 = 1.0954, where the sign √ stands for the square root.

2. The Monte Carlo estimate

Run the model M times with fresh random numbers, and read one number X_m off path m, such as its final value. Call the true value we want θ (theta), here the expected final value E[X]. The Monte Carlo estimate θ_M, with the subscript M for the number of paths, is the plain average of the path results.

The expected value of a sum is the sum of the expected values (module 3), so the estimate has no bias: it is right on average. The law of large numbers of module 3 pulls it closer to θ as M grows. Example: the 10,000 final values of Figure 1 average $111,469.

3. The exact answers of this model

This model can also be solved with pen and paper, which lets us check the simulation. The percentile x_u at the share u is the value with the share u of the outcomes below it, so the 10th percentile has u = 0.1 and the median u = 0.5. Φ(z), the capital Greek letter phi, is the share of the standard bell below z (module 7), and z_u is the point with Φ(z_u) = u, for example z_0.9 = 1.2816. The exponential keeps the order of values, so every percentile of the log turns into a percentile of the value. Module 10 showed that a bell Y with the mean a and the standard deviation c has E[e^Y] = e^a + c2/2, and here a + c²/2 = gT + σ²T/2 = μT.

Example: the median is 10,000 e^1.8 = $60,496 (simulated $60,814), the 10th percentile 10,000 e^{1.8 − 1.2816 × 1.0954} = $14,861 (simulated $14,636) and the 90th percentile $246,276 (simulated $246,456). The chance of a loss is Φ(−1.8/1.0954) = Φ(−1.6432) = 5.02 % (simulated 5.00 %, 500 paths), and the average is 10,000 e^2.4 = $110,232 (simulated $111,469). The share of outcomes below the average is Φ(σ√T/2) = Φ(0.5477) = 70.81 % (simulated 70.97 %). Every simulated number lies near its exact value, and step 4 says how near it should be.

4. The standard error

How far is θ_M from θ? Call σ_X the standard deviation of one path result. The paths are independent, so the variance of their sum is M times σ_X² (module 6), and dividing by M divides the variance by M². The standard error SE is the standard deviation of the estimate, the typical distance between the estimate and the true value. The double arrow ⟹ reads “and therefore”, and the sign ≈ means “about equal”.

The true σ_X is usually unknown, so a run measures its own standard deviation s, computed with M − 1 as in module 8, and reports s/√M. Example: the paths of Figure 1 have s = $167,084, so SE = $167,084/√10,000 = $1,671. The variance formula of module 10 gives the exact σ_X = E[W_T] × √(e^σ2T − 1) = $167,904 and the exact SE $1,679. Four times as many paths halve the error: 40,000 paths give $840. A hundred times as many make it ten times smaller: 100 paths give $16,790 and 1,000,000 paths $168.

5. A 95 % interval

By the central limit theorem of module 7, the average of many paths is close to a bell around θ with the standard deviation SE. A bell holds 95 % of its area within 1.96 standard deviations of its middle (module 7). So the interval from θ_M − 1.96 SE to θ_M + 1.96 SE holds θ in about 95 % of all runs. The sign ± means “plus or minus” and writes both ends at once.

Example: $111,469 ± 1.96 × $1,671 runs from $108,194 to $114,743, and it holds the exact $110,232. The estimate lies (111,469 − 110,232)/1,671 = +0.74 standard errors above it. The bell of the central limit theorem needs enough paths, the more so the more lopsided the results are. Exercise 4 repeats a simulation 1,000 times: for the final value the interval holds the truth in 89.9 % of the runs with 100 paths and in 94.0 % with 1,000, for the symmetric log in 94.2 % and 94.4 %.

6. A chance is an average too

A chance is the average of a flag that is 1 when the event happens and 0 otherwise. Call the true chance p and its estimate from M paths p_M. A flag times itself is the flag again, so the expected square of a flag is p, and the variance rule of module 6 gives it the variance p − p² = p(1 − p). Step 4 then gives the standard error, and divided by p it shows the error as a share of the chance itself.

Example: the simulation of Figure 3 estimates the chance of a loss as 5.00 %, with SE = √(0.05 × 0.95/10,000) = 0.218 %. The 95 % interval runs from 4.57 % to 5.43 % and holds the exact 5.02 %. As a share of the chance the error is 4.4 %. A chance of 1 % from 10,000 paths rests on only about 100 hits, and its error is 9.9 % of its size. For the same relative error, the number of paths must grow about like 1/p, so rare events, such as the very bad days of module 18, need many paths.

7. The error of a percentile

A percentile moves when paths cross it. At the true percentile x_u, the share of paths below it is an estimated chance with the standard error √(u(1 − u)/M) of step 6. Let f be the density of the outcomes, the height of the curve whose area gives the chances (module 7). Moving a value by a small step δ (delta) changes the share below it by about f(x_u) δ, so an error in the share turns into an error in the value of about the share error divided by f(x_u). For the median of the log, a bell with the standard deviation σ√T, the height in the middle is 1/(σ√T × √(2π)), where π = 3.14159 is the circle number, and u = 1/2. That gives the second formula.

Example: σ√T = 1.0954, so the median of 10,000 paths is uncertain by 1.37 %, the median of 1,000 paths by 4.34 % and the median of 100 paths by 13.73 %. The run of Figure 1 agrees, as Figure 7 shows: the medians of 100 groups of 100 of its paths spread by 13.64 %, from $41,680 to $86,829, the medians of 10 groups of 1,000 by 5.31 %, a rough value from only ten groups, and the median of all 10,000 lies 0.52 % above the exact $60,496.

8. Seeds, and checking the error with other runs

A computer’s random numbers come from a fixed rule that starts at the seed, so they are called pseudo-random. The same seed gives the same paths, and a different seed, or the same generator used again, gives an independent run. Independent runs are a direct check of the standard error: the distance of an estimate from the truth, counted in standard errors, behaves like a draw from the standard bell. The bars | | stand for the size of a number without its sign.

Example: the seed 17 run again gives every path exactly as before. The seeds 1 to 5 give medians from $58,702 to $62,463, between −2.19 and +2.33 standard errors from the exact median, averages from $108,005 to $111,644 and chances of a loss from 4.66 % to 5.24 %. Two of the five medians lie a little more than two standard errors away. A single run lands that far out with a chance of about 5 %, so two of five happen by luck alone in about 2 % of such sets of runs, and exercise 2 checks the error over 100 runs. A seed is fixed before the run and never picked after a look at the results, because trying seeds until one looks good turns luck into a false result.

9. Antithetic pairs

The error can also shrink with smarter random numbers. Each path is driven by draws Z from the standard bell, and the mirrored draws −Z are exactly as likely. Use both and count the pair as one value, the average of the results X(Z) and X(−Z). This is the method of antithetic pairs, the simplest of the methods that reduce the variance of an estimate. With the correlation ρ (rho) of the two halves of a pair (module 6), the variance of the sum of two equal variances gives the formula below. M paths make M/2 pairs.

For the final value W₀e^a + cZ and its mirror, the correlation is ρ = −e^−c2, since their covariance e^2a(1 − e^c2) divided by the variance e^2a + c2(e^c2 − 1) leaves −e^−c2. With c² = σ²T = 1.2 this is −0.301, and the error shrinks by √(1 − 0.301) = 0.836, from $1,679 to $1,404 for 10,000 paths. Part 8 of the full code measures a correlation of −0.300 and a standard error of $1,384. For the log of the final value, a straight line in Z, ρ = −1 and the pairs give the exact 1.8 every time. For the chance of a loss the gain is tiny, as exercise 5 shows. Pairs help most for results that move steadily with the random numbers.

Best practice

• Results are usually reported with their standard error or an interval, together with the seed and the number of paths, so that others can repeat them.
• A simulation is only as good as its inputs: the assumed average, the volatility and the shape of the returns decide what comes out.

Common mistakes

• Reporting an estimate without its standard error. The average of $111,469 carries an error of $1,671, and from 100 paths it would carry about $16,790.
• Expecting more paths to fix the model. The standard error measures only the luck of the draws. If real returns have fat tails (module 17), a million paths of a bell answer the wrong question very precisely.
• Picking the seed after seeing the results. Five seeds already give medians from $58,702 to $62,463. Choosing the most pleasing one turns luck into a false result.
• Treating two runs with the same seed as independent. Two generators with the seed 2 gave the same median of $63,401.51 twice, so the second run adds no information.
• Reading the average as the typical outcome. 70.97 % of the paths of Figure 1 end below their average of $111,469, and the median is $60,814.
• Trusting a 95 % interval from few paths of a lopsided result. With 100 paths the interval of the final value held the truth in only 89.9 % of 1,000 runs.
• Using too few paths for a rare event. A chance of 1 % from 10,000 paths is uncertain by 9.9 % of its size.
• Counting an antithetic pair as two independent paths. The two halves are correlated, so each pair is one value, and the standard error comes from the pair values.

Recap

• A Monte Carlo simulation runs a model many times with random numbers and reads every answer off the paths.
• The estimate θ_M is the average of the M path results, and it is right on average.
• Its standard error is σ_X/√M, estimated by s/√M: four times the paths halve it, a hundred times make it ten times smaller.
• θ_M ± 1.96 SE holds the truth in about 95 % of runs, once there are enough paths for the result at hand.
• A chance is the average of a 0/1 flag, with the standard error √(p(1 − p)/M), and rare events need many paths.
• A percentile is uncertain by the error of a share divided by the density there, √(π/2) σ√T/√M for the median of the log.
• A fixed seed makes a run repeatable, and independent runs check the standard error.
• Antithetic pairs shrink the error by √(1 + ρ), which helps when the result moves steadily with the random numbers.

Exercises

Three short calculations first, with pen, paper and a calculator, then five programs with rising difficulty. Try each one before you open its solution.

In Python

Each file runs as it is and prints that it is not solved yet. Fill in the lines marked TODO, run the file again and compare your numbers with the solution.

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Exercise 1 (easy).
# Runs on its own with Python 3 and numpy. Fill in the lines marked TODO.
#
# Throw two dice M times and estimate the chance that their sum is at least
# 10. The exact chance is 6/36 = 1/6: (4,6), (5,5), (6,4), (5,6), (6,5), (6,6).
# a) Write share_and_se(flags): the share of True values and its standard
#    error sqrt(p (1 - p) / M), where M is the number of flags.
# b) Run it for M = 100, 10,000 and 1,000,000 throws.
# c) How far is each estimate from 1/6, counted in standard errors?
import math
import numpy as np

rng = np.random.default_rng(1)                  # fixed seed

def share_and_se(flags):
    # TODO: return (share, standard error) as two floats
    return None

if share_and_se(np.array([True, False])) is None:
    print("Exercise 1 is not solved yet. Fill in the lines marked TODO.")
else:
    for m in (100, 10000, 1000000):
        dice = rng.integers(1, 7, size=(m, 2))      # whole numbers 1 to 6
        share, se = share_and_se(dice.sum(axis=1) >= 10)
        print(f"M = {m:>9,}:  estimate {share:.4%}   SE {se:.4%}   "
              f"off by {(share - 1 / 6) / se:+.2f} SE")

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Solution to exercise 1.
# Runs on its own with Python 3 and numpy.
import math
import numpy as np

rng = np.random.default_rng(1)                  # fixed seed

def share_and_se(flags):
    # a share is the average of 0/1 flags, so its SE is sqrt(p (1 - p) / M)
    flags = np.asarray(flags, dtype=float)
    p = float(flags.mean())
    return p, math.sqrt(p * (1 - p) / flags.size)

for m in (100, 10000, 1000000):
    dice = rng.integers(1, 7, size=(m, 2))      # whole numbers 1 to 6
    share, se = share_and_se(dice.sum(axis=1) >= 10)
    print(f"M = {m:>9,}:  estimate {share:.4%}   SE {se:.4%}   "
          f"off by {(share - 1 / 6) / se:+.2f} SE")

M =       100:  estimate 18.0000%   SE 3.8419%   off by +0.35 SE
M =    10,000:  estimate 16.4400%   SE 0.3706%   off by -0.61 SE
M = 1,000,000:  estimate 16.7142%   SE 0.0373%   off by +1.27 SE

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Exercise 2 (easy).
# Runs on its own with Python 3 and numpy. Fill in the lines marked TODO.
#
# The model: $10,000 over 30 years, yearly log returns from a bell
# with the mean g = 0.08 - 0.2^2 / 2 = 0.06 and the sd 0.2.
# a) Write median_of_run(rng, paths): simulate paths final values with the
#    generator rng and return their median.
# b) Two generators with the same seed: do they give the same median?
# c) One generator used again gives a new run. Make 100 runs of 1,000 paths
#    and compare the spread of the logs of their medians with the theory
#    sqrt(pi / 2) * 0.2 * sqrt(30) / sqrt(1000).
import math
import numpy as np

start, years, g, sigma = 10000, 30, 0.06, 0.2
rng_a = np.random.default_rng(2)                # two generators, same seed
rng_b = np.random.default_rng(2)

def median_of_run(rng, paths):
    # TODO: draw (paths, years) log returns, add them up per path,
    #       turn them into final values and return their median
    return None

if median_of_run(np.random.default_rng(0), 10) is None:
    print("Exercise 2 is not solved yet. Fill in the lines marked TODO.")
else:
    first, twin = median_of_run(rng_a, 1000), median_of_run(rng_b, 1000)
    print(f"seed 2 in two generators: ${first:,.2f} and ${twin:,.2f}")
    print(f"the first generator again: ${median_of_run(rng_a, 1000):,.2f}")
    logs = np.log([median_of_run(rng_a, 1000) for _ in range(100)])
    theory = math.sqrt(math.pi / 2) * sigma * math.sqrt(years) / math.sqrt(1000)
    print(f"100 runs of 1,000 paths: medians from ${math.exp(logs.min()):,.0f} "
          f"to ${math.exp(logs.max()):,.0f}")
    print(f"spread of their logs {logs.std(ddof=1):.2%}   theory {theory:.2%}")

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Solution to exercise 2.
# Runs on its own with Python 3 and numpy.
import math
import numpy as np

start, years, g, sigma = 10000, 30, 0.06, 0.2
rng_a = np.random.default_rng(2)                # two generators, same seed
rng_b = np.random.default_rng(2)

def median_of_run(rng, paths):
    # every call takes the next random numbers of rng: a new, independent run
    log_ret = rng.normal(g, sigma, (paths, years))
    final = start * np.exp(log_ret.sum(axis=1))
    return float(np.median(final))

first, twin = median_of_run(rng_a, 1000), median_of_run(rng_b, 1000)
print(f"seed 2 in two generators: ${first:,.2f} and ${twin:,.2f}")
print(f"the first generator again: ${median_of_run(rng_a, 1000):,.2f}")
logs = np.log([median_of_run(rng_a, 1000) for _ in range(100)])
theory = math.sqrt(math.pi / 2) * sigma * math.sqrt(years) / math.sqrt(1000)
print(f"100 runs of 1,000 paths: medians from ${math.exp(logs.min()):,.0f} "
      f"to ${math.exp(logs.max()):,.0f}")
print(f"spread of their logs {logs.std(ddof=1):.2%}   theory {theory:.2%}")

seed 2 in two generators: $63,401.51 and $63,401.51
the first generator again: $58,902.18
100 runs of 1,000 paths: medians from $52,655 to $66,796
spread of their logs 4.13%   theory 4.34%

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Exercise 3 (medium).
# Runs on its own with Python 3 and numpy. Fill in the lines marked TODO.
#
# 10,000 simulated paths of $10,000 (8 % a year, swings of 20 %), read after
# 1, 5, 10, 20 and 30 years.
# a) Write exact_percentile(t, u): after t years ln(W / W0) follows a bell
#    with the mean g t and the sd sigma sqrt(t), so the percentile at the
#    share u is W0 exp(g t + z sigma sqrt(t)), z = NormalDist().inv_cdf(u).
# b) Write exact_loss(t): the chance to end below W0, which is
#    NormalDist().cdf(-g t / (sigma sqrt(t))).
# c) Compare both with the simulated paths, and watch the ratio of the 90th
#    to the 10th percentile grow with the square root of time.
import math
from statistics import NormalDist
import numpy as np

start, years, paths, g, sigma = 10000, 30, 10000, 0.06, 0.2
rng = np.random.default_rng(3)                  # fixed seed
value = start * np.exp(rng.normal(g, sigma, (paths, years)).cumsum(axis=1))

def exact_percentile(t, u):
    # TODO: start * exp(g t + z sigma sqrt(t)) with z of the bell at u
    return None

def exact_loss(t):
    # TODO: the bell's share below -g t / (sigma sqrt(t))
    return None

if exact_percentile(1, 0.5) is None or exact_loss(1) is None:
    print("Exercise 3 is not solved yet. Fill in the lines marked TODO.")
else:
    print(" t   10th sim / exact      median sim / exact     loss sim / exact")
    for t in (1, 5, 10, 20, 30):
        lo, mid = np.percentile(value[:, t - 1], [10, 50])
        print(f"{t:>2}  {lo:>8,.0f} / {exact_percentile(t, 0.1):>8,.0f}"
              f"   {mid:>8,.0f} / {exact_percentile(t, 0.5):>8,.0f}"
              f"   {(value[:, t - 1] < start).mean():6.2%} / {exact_loss(t):6.2%}")
    for t in (1, 4, 9, 16, 30):
        hi, lo = np.percentile(value[:, t - 1], [90, 10])
        print(f"t = {t:>2}: 90th / 10th = {hi / lo:5.2f}   exact "
              f"{exact_percentile(t, 0.9) / exact_percentile(t, 0.1):5.2f}")

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Solution to exercise 3.
# Runs on its own with Python 3 and numpy.
import math
from statistics import NormalDist
import numpy as np

start, years, paths, g, sigma = 10000, 30, 10000, 0.06, 0.2
rng = np.random.default_rng(3)                  # fixed seed
value = start * np.exp(rng.normal(g, sigma, (paths, years)).cumsum(axis=1))

def exact_percentile(t, u):
    # the log grows by g t on average and spreads by sigma sqrt(t)
    z = NormalDist().inv_cdf(u)
    return start * math.exp(g * t + z * sigma * math.sqrt(t))

def exact_loss(t):
    # below the start means ln(W / W0) below 0
    return NormalDist().cdf(-g * t / (sigma * math.sqrt(t)))

print(" t   10th sim / exact      median sim / exact     loss sim / exact")
for t in (1, 5, 10, 20, 30):
    lo, mid = np.percentile(value[:, t - 1], [10, 50])
    print(f"{t:>2}  {lo:>8,.0f} / {exact_percentile(t, 0.1):>8,.0f}"
          f"   {mid:>8,.0f} / {exact_percentile(t, 0.5):>8,.0f}"
          f"   {(value[:, t - 1] < start).mean():6.2%} / {exact_loss(t):6.2%}")

for t in (1, 4, 9, 16, 30):
    hi, lo = np.percentile(value[:, t - 1], [90, 10])
    print(f"t = {t:>2}: 90th / 10th = {hi / lo:5.2f}   exact "
          f"{exact_percentile(t, 0.9) / exact_percentile(t, 0.1):5.2f}")

 t   10th sim / exact      median sim / exact     loss sim / exact
 1     8,214 /    8,218     10,654 /   10,618   37.56% / 38.21%
 5     7,682 /    7,610     13,560 /   13,499   24.83% / 25.12%
10     8,282 /    8,102     18,161 /   18,221   17.03% / 17.14%
20    10,725 /   10,552     33,295 /   33,201    8.60% /  8.99%
30    15,414 /   14,861     61,034 /   60,496    4.66% /  5.02%
t =  1: 90th / 10th =  1.68   exact  1.67
t =  4: 90th / 10th =  2.81   exact  2.79
t =  9: 90th / 10th =  4.63   exact  4.65
t = 16: 90th / 10th =  7.66   exact  7.77
t = 30: 90th / 10th = 16.42   exact 16.57

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Exercise 4 (medium).
# Runs on its own with Python 3 and numpy. Fill in the lines marked TODO.
#
# How often does a 95 % interval hold the true value? In the model of step 1,
# ln(W / W0) after 30 years follows a bell with the mean 1.8 and the
# sd 0.2 sqrt(30), so one draw per path is enough. The exact average final
# value is 10,000 e^(0.08 * 30).
# a) Write interval(values): the average plus and minus 1.96 standard
#    errors, with the standard error s / sqrt(M) and s computed with M - 1.
# b) Repeat a simulation of M paths 1,000 times, for M = 100 and 1,000, and
#    count how often the interval holds the exact average of the final value
#    and how often the interval of ln(W / W0) holds the exact 1.8.
import math
import numpy as np

rng = np.random.default_rng(4)                  # fixed seed
start, mean_log, sd_log = 10000, 1.8, 0.2 * math.sqrt(30)
exact_final = start * math.exp(0.08 * 30)

def interval(values):
    # TODO: return (low, high), the average -/+ 1.96 standard errors
    return None

if interval(np.array([1.0, 2.0, 3.0])) is None:
    print("Exercise 4 is not solved yet. Fill in the lines marked TODO.")
else:
    for m in (100, 1000):
        hit_final = hit_log = 0
        for _ in range(1000):
            logs = rng.normal(mean_log, sd_log, m)
            low, high = interval(start * np.exp(logs))
            hit_final += low <= exact_final <= high
            low, high = interval(logs)
            hit_log += low <= mean_log <= high
        print(f"M = {m:>5,}: the interval of the final value holds it in "
              f"{hit_final / 1000:.1%}, the interval of the log in "
              f"{hit_log / 1000:.1%}")

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Solution to exercise 4.
# Runs on its own with Python 3 and numpy.
import math
import numpy as np

rng = np.random.default_rng(4)                  # fixed seed
start, mean_log, sd_log = 10000, 1.8, 0.2 * math.sqrt(30)
exact_final = start * math.exp(0.08 * 30)

def interval(values):
    # the central limit theorem makes the average about a bell around the truth
    values = np.asarray(values, dtype=float)
    se = values.std(ddof=1) / math.sqrt(values.size)
    return values.mean() - 1.96 * se, values.mean() + 1.96 * se

for m in (100, 1000):
    hit_final = hit_log = 0
    for _ in range(1000):
        logs = rng.normal(mean_log, sd_log, m)
        low, high = interval(start * np.exp(logs))
        hit_final += low <= exact_final <= high
        low, high = interval(logs)
        hit_log += low <= mean_log <= high
    print(f"M = {m:>5,}: the interval of the final value holds it in "
          f"{hit_final / 1000:.1%}, the interval of the log in "
          f"{hit_log / 1000:.1%}")

M =   100: the interval of the final value holds it in 89.9%, the interval of the log in 94.2%
M = 1,000: the interval of the final value holds it in 94.0%, the interval of the log in 94.4%

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Exercise 5 (hard).
# Runs on its own with Python 3 and numpy. Fill in the lines marked TODO.
#
# Antithetic pairs: every standard bell draw z is used twice, as z and -z.
# In the model of step 1, W / W0 = exp(1.8 + 0.2 sqrt(30) z) after 30 years.
# a) Write plain(func, z): the average of func(z) and its standard error.
# b) Write antithetic(func, z): the average of the pair values
#    (func(z) + func(-z)) / 2 and their standard error, one value per pair.
# c) Compare both with the same number of paths, 10,000 plain draws against
#    5,000 pairs, for the final value, its log and the chance of a loss.
#    Then repeat both 200 times and measure the real spread of each estimate.
import math
import numpy as np

rng = np.random.default_rng(5)                  # fixed seed
start, mean_log, sd_log = 10000, 1.8, 0.2 * math.sqrt(30)
tasks = [("final value", ",.0f", lambda z: start * np.exp(mean_log + sd_log * z)),
         ("log of W / W0", ".4f", lambda z: mean_log + sd_log * z),
         ("chance of a loss", ".4f", lambda z: (mean_log + sd_log * z < 0) * 1.0)]

def plain(func, z):
    # TODO: (average of func(z), its standard error)
    return None

def antithetic(func, z):
    # TODO: pair values (func(z) + func(-z)) / 2, then (average, standard error)
    return None

line = tasks[1][2]
if plain(line, np.ones(4)) is None or antithetic(line, np.ones(4)) is None:
    print("Exercise 5 is not solved yet. Fill in the lines marked TODO.")
else:
    for name, spec, func in tasks:
        a, se_a = plain(func, rng.standard_normal(10000))
        b, se_b = antithetic(func, rng.standard_normal(5000))
        print(f"{name:<17} plain {a:{spec}} (SE {se_a:{spec}})   antithetic "
              f"{b:{spec}} (SE {se_b:{spec}})")
    for name, spec, func in tasks:
        runs_a = [plain(func, rng.standard_normal(10000))[0] for _ in range(200)]
        runs_b = [antithetic(func, rng.standard_normal(5000))[0] for _ in range(200)]
        print(f"{name:<17} spread of 200 runs: plain {np.std(runs_a, ddof=1):{spec}}"
              f"   antithetic {np.std(runs_b, ddof=1):{spec}}")

# QuantFinanceToGo, module 16: Monte Carlo Simulation. Solution to exercise 5.
# Runs on its own with Python 3 and numpy.
import math
import numpy as np

rng = np.random.default_rng(5)                  # fixed seed
start, mean_log, sd_log = 10000, 1.8, 0.2 * math.sqrt(30)
tasks = [("final value", ",.0f", lambda z: start * np.exp(mean_log + sd_log * z)),
         ("log of W / W0", ".4f", lambda z: mean_log + sd_log * z),
         ("chance of a loss", ".4f", lambda z: (mean_log + sd_log * z < 0) * 1.0)]

def plain(func, z):
    # one value per draw, SE = s / sqrt(M)
    x = func(z)
    return float(x.mean()), float(x.std(ddof=1) / math.sqrt(x.size))

def antithetic(func, z):
    # the two halves of a pair are not independent, so the pair is one value
    x = (func(z) + func(-z)) / 2
    return float(x.mean()), float(x.std(ddof=1) / math.sqrt(x.size))

for name, spec, func in tasks:
    a, se_a = plain(func, rng.standard_normal(10000))
    b, se_b = antithetic(func, rng.standard_normal(5000))
    print(f"{name:<17} plain {a:{spec}} (SE {se_a:{spec}})   antithetic "
          f"{b:{spec}} (SE {se_b:{spec}})")

for name, spec, func in tasks:
    runs_a = [plain(func, rng.standard_normal(10000))[0] for _ in range(200)]
    runs_b = [antithetic(func, rng.standard_normal(5000))[0] for _ in range(200)]
    print(f"{name:<17} spread of 200 runs: plain {np.std(runs_a, ddof=1):{spec}}"
          f"   antithetic {np.std(runs_b, ddof=1):{spec}}")

final value       plain 113,128 (SE 1,674)   antithetic 110,737 (SE 1,406)
log of W / W0     plain 1.8020 (SE 0.0109)   antithetic 1.8000 (SE 0.0000)
chance of a loss  plain 0.0498 (SE 0.0022)   antithetic 0.0532 (SE 0.0022)
final value       spread of 200 runs: plain 1,717   antithetic 1,355
log of W / W0     spread of 200 runs: plain 0.0110   antithetic 0.0000
chance of a loss  spread of 200 runs: plain 0.0023   antithetic 0.0024

Interview prep

Questions on this module in the style of quant interviews, written in my own words and ordered from easy to hard. Write down your own answer first, then open the worked answer and compare.