Monte Carlo Simulation

Quants and system developers who calculate derivative prices use this method to value products that cannot be solved with a formula. Risk managers use it to find the future exposure (a value found from the distribution of the value of a trade assessed at a future date) in counterparty credit risk calculations, and to calculate risk measures such as value at risk (VaR)¹².

What the Monte Carlo method calculates

In financial practice, the Monte Carlo method means repeating numerical experiments with random numbers many times and finding a target value from the many results¹. For the pricing covered on this page, it approximates the price of a derivative by the average of many results generated with random numbers. Take the call option from the example used throughout the option valuation chapter (underlying asset price 100, strike price 100, 1 year to expiry, interest rate 1%, volatility 20%). Valued with 10,000 paths, the estimate is 8.33. The value given by the Black-Scholes model formula (the analytic price) is 8.43, a difference of 0.10³.

Under risk-neutral valuation, the price is the discounted expected value of the payoff (the amount received) at expiry. The Monte Carlo method replaces this expected value with the average payoff over paths generated with random numbers. The paths are therefore generated with the expected return of the underlying set to the risk-free rate⁴. How this setting differs from the scenarios used to forecast future values is covered in Risk-Neutral Valuation.

In financial practice, the Monte Carlo method is used both to find the mean of a distribution, as with derivative prices, and to find a point in the tail of a distribution, as with value at risk (VaR). Both follow the same steps: generate changes in the underlying price with random numbers, calculate values or profits and losses from them, and aggregate the results¹. The rest of the page takes up pricing, which finds the mean.

How paths are generated and the price is found

The calculation has three steps: generate random numbers, generate paths, and discount and average the payoffs. Uniform random numbers between 0 and 1 are converted into standard normal random numbers, which follow the normal distribution with mean 0 and variance 1. These are used to generate prices of the underlying, from which the value of the derivative is calculated. Repeating this as many times as needed and taking the average gives the price estimate⁵.

Geometric Brownian motion is a process in which the log of the price follows Brownian motion with drift. When the price of the underlying follows geometric Brownian motion, the price $\Delta t$ years after time $t$ can be generated from the price $S_t$ with one normal random number $Z$, using the formula below⁶. $r$ is the interest rate (continuously compounded) and $\sigma$ is the volatility.

$$S_{t+\Delta t} = S_t \exp\!\left( \left( r - \tfrac{1}{2}\sigma^2 \right) \Delta t + \sigma \sqrt{\Delta t}\, Z \right)$$

This formula directly uses the fact that the increment in the log price is normally distributed, so the distribution stays exact even when $\Delta t$ is large⁶. For a European option, whose payoff depends only on the price at expiry, you can set $\Delta t$ to the time to expiry $T$ and reach expiry in a single step, which is what the simulator on this page does. For products that look at prices along the way, you advance one step at a time to each observation date.

If the prices at expiry on $N$ paths are $S_T^{(1)}, \dots, S_T^{(N)}$, the price estimate $\hat{C}$ for a call with strike $K$ is the average of the discounted payoffs $Y_i$ on the paths.

$$\hat{C} = \frac{1}{N} \sum_{i=1}^{N} Y_i, \qquad Y_i = e^{-rT} \max\!\left(S_T^{(i)} - K,\ 0\right)$$

For one path, check the calculation with the common example ($S_0 = 100$, $K = 100$, $r = 1\%$, $\sigma = 20\%$, $T = 1$). If the normal random number is $Z = 1$, the price at expiry is $100 \exp(0.01 - 0.02 + 0.2) = 100 \exp(0.19)$, about 120.925. The payoff is 120.925 − 100 = 20.925, and the discounted value $Y$ is $20.925 \times e^{-0.01}$, about 20.72. If $Z = -1$, the price at expiry is $100 \exp(-0.21)$, about 81.06, which is below the strike price, so the payoff is 0 and so is $Y$. The average of these two paths alone is about 10.36, far from the analytic price of 8.43. The value on a single path varies widely, from 0 to large values, so you take the average over more paths.

Because the method only calculates the payoff on each path, the same steps can value a path-dependent option, whose payoff changes with the price movements before expiry⁷.

Worked example: number of paths and the estimate

Valuing the call from the common example with different numbers of paths gives the estimates and standard errors in the table below. The random seed is 1, and every row uses the same sequence of random numbers from the start³.

Quadrupling the paths from 1,000 to 4,000 roughly halves the standard error, from 0.43 to 0.21. Nine times as many paths (9,000) cut it to 0.14, roughly a third. Meanwhile, the difference between the estimate and the analytic price stays around one standard error in size. More paths bring the estimate closer to the analytic price, but the difference moves up and down with the random numbers, so you cannot judge accuracy from the difference in a single run. Accuracy is expressed by the standard error.

Standard error and the 95% range

The standard error measures how much the estimate varies because of the random numbers. It is found by dividing the standard deviation $s$ of the discounted payoffs on the paths by the square root of the number of paths.

$$\mathrm{SE} = \frac{s}{\sqrt{N}}, \qquad s^2 = \frac{1}{N-1} \sum_{i=1}^{N} \left( Y_i - \hat{C} \right)^2$$

$s$ barely changes as the number of paths grows, and in the common example it is around 13.4³. The standard error therefore shrinks in proportion to $1/\sqrt{N}$: halving it takes four times as many paths, and gaining one more digit of accuracy takes 100 times as many⁸. To bring the standard error in the common example down to 0.01 takes (13.4 ÷ 0.01) squared, or about 1.8 million paths³. Because $s$ is the spread of the payoffs on the paths, it depends on the product's terms. With volatility at 40%, $s$ is about 30, and the standard error with 10,000 paths is 0.30³.

When the number of paths is large, the range of 1.96 standard errors either side of the estimate contains the true price under the model about 95% of the time. The true price is the risk-neutral expected value, which in the common example is the analytic price. 1.96 is the 97.5th percentile of the standard normal distribution⁹. For the common example with 10,000 paths, 8.33 ± 1.96 × 0.13 gives a range of 8.08 to 8.58, which contains the analytic price of 8.43³. Building the range the same way for each seed from 1 to 100, the analytic price fell inside it in 96 of the 100 runs. In the other 4, the difference between the estimate and the analytic price was still within 3 standard errors³.

Deep diveWhy the standard error is inversely proportional to the square root of the number of paths

When the values on the paths $Y_1, \dots, Y_N$ are independent and identically distributed, the variance of their average $\hat{C}$ equals the variance of a single value $\mathrm{Var}(Y)$ divided by $N$¹⁰.

$$\mathrm{Var}(\hat{C}) = \frac{\mathrm{Var}(Y)}{N}$$

The standard error is the square root of this, $\sqrt{\mathrm{Var}(Y)} / \sqrt{N}$. Since $\mathrm{Var}(Y)$ is unknown, actual calculations replace it with $s^2$, computed from the values on the paths.

Random numbers and seeds

The random numbers used in the Monte Carlo method are pseudo-random numbers, which a computer generates one after another with a formula. A sequence is often regarded as random enough, and called pseudo-random, when its values do not depend on earlier values (statistical independence) and appear without bias (uniformity)¹¹. For example, the linear congruential generator produces the next value from the previous one with the formula $X_n = aX_{n-1} + c \pmod m$¹².

Because the numbers come from a formula, the same starting value produces the same sequence. In this site's simulator, the input that sets this starting value is called the seed, and the same seed gives the same estimate. A different seed changes the sequence, and the estimate moves within a range about the size of the standard error.

Random numbers generated by a formula return to the same sequence after a fixed period. Using more random numbers than the period only reuses the same numbers and does not improve the accuracy of the estimate¹². A 2002 paper from the Bank of Japan's Institute for Monetary and Economic Studies showed that the built-in random number generator of a spreadsheet program had a short period that may not be enough for large calculations. It gave the Mersenne Twister, with a long period of $2^{19937}-1$, as an example¹³.

Variance reduction: smaller errors without more paths

Variance reduction methods use random numbers in ways that reduce the spread of the values on the paths ($s$ in the standard error formula). They leave the $1/\sqrt{N}$ rate unchanged and reduce its proportionality constant¹⁵. Two common variance reduction methods are antithetic variates and control variates.

The antithetic variates method pairs a path generated with a normal random number $Z$ with a path generated with the opposite sign, $-Z$, and uses the average of the two values as one value¹⁵. The more the two values in a pair move in opposite directions, the smaller the variance of their average. When the payoff is monotonically increasing or decreasing in the random number (as with a European call), the covariance of the two values is negative, which guarantees a smaller variance¹⁶. However, antithetic variates rarely reduce the variance by a large amount¹⁷.

The control variates method calculates, on the same paths, a value $W$ that moves together with the value to be estimated $Y$ and whose expected value is known from a formula. It then uses the error in the estimate of $W$ to correct the estimate of $Y$¹⁸. If $c$ in the formula below is set from the covariance of $Y$ and $W$, the stronger the correlation between $Y$ and $W$, the smaller the variance¹⁹.

$$\hat{C}_{\mathrm{cv}} = \frac{1}{N} \sum_{i=1}^{N} \left( Y_i + c \left( W_i - E[W] \right) \right), \qquad c = -\frac{\mathrm{Cov}(Y, W)}{\mathrm{Var}(W)}$$

For example, when valuing an Asian option, whose payoff is set by the average price over a period, you can use as $W$ the discounted payoff of an ordinary European call calculated on the same paths. The expected value of $W$ is the price given by the Black-Scholes model formula²⁰.

When choosing a method, compare not only the variance but also the computational effort per path. A method that halves the variance but takes more than twice as long per path gives lower accuracy for the same computing time²¹. Some methods also reduce the number of trials per unit of time and can increase the computing time²².

Deep diveWhy antithetic variates reduce the variance

Let $Y$ be the value on the path generated with $Z$ and $\tilde{Y}$ the value on the path generated with $-Z$. The variance of the average of the two is as follows²³.

$$\mathrm{Var}\!\left( \frac{Y + \tilde{Y}}{2} \right) = \frac{1}{4} \left( \mathrm{Var}(Y) + \mathrm{Var}(\tilde{Y}) + 2\,\mathrm{Cov}(Y, \tilde{Y}) \right)$$

If $Y$ and $\tilde{Y}$ are independent, the covariance is 0 and the variance is $\mathrm{Var}(Y)/2$, the same as for the average of two independent paths. If the covariance is negative, the variance is smaller than that. Pairing $Z$ with $-Z$ also makes the mean of the random numbers used exactly 0, and makes the odd-order moments 0 as well²³.

Where it is used and how it compares with other methods

The Monte Carlo method needs more computation than other methods to reach high accuracy. According to a 2002 paper from the Bank of Japan's Institute for Monetary and Economic Studies, practitioners chose methods in this order: analytic formulas such as the Black-Scholes model, then lattice methods such as the binomial model and finite difference methods, then the Monte Carlo method. In the front office in particular, where fast calculation is required, the Monte Carlo method was in many cases regarded as a last resort⁷.

The same paper states that, among options exercisable only at expiry, the Monte Carlo method is used in the following cases⁷: options whose model of underlying price movements is complex enough to make other methods laborious to build, options with path dependence, and options on four or more underlyings. For lattice and finite difference methods, the amount of computation rises sharply as the number of underlyings grows, and about three underlyings is the practical limit. With the Monte Carlo method, adding an underlying only adds random numbers in proportion to the number of paths²⁴.

For a barrier option, an example of a path-dependent product, the price on each path is checked at every date on which the barrier is observed. Take a down-and-out call with an underlying price of 100, a strike price of 100, and a barrier of 90, observed 50 times before expiry. The estimate with 10,000 paths is 7.44 (standard error 0.13). This is higher than the analytic price of 6.88, which assumes continuous monitoring, and close to the analytic price of 7.32 adjusted for the number of observations³. The adjustment shifts the barrier in the continuous-monitoring formula away from the underlying, to suit a contract monitored only on set dates²⁵. Because it moves the barrier from 90 to 88.53, the adjusted analytic price of 7.32 is higher than 6.88.

Besides derivative prices, the method is also used to calculate counterparty credit risk. The expected exposure used in XVA calculations can be found with simplified methods, but the Monte Carlo method is used to find it precisely while capturing complex trades and the effect of netting. Its drawback is the heavy computational load². The whole picture of how the exposure feeds into valuation is covered in Overview of XVA.

Calculations that find a point in the tail of a distribution, such as VaR, can have larger errors than calculations of the mean for the same number of paths. In the normal distribution example in the same 2002 Institute for Monetary and Economic Studies paper, the error of the 99.5th percentile was larger than the error of the mean. Most variance reduction methods are techniques for reducing the error in calculations of the mean, and they cannot always be applied as they are to calculations of tail points²⁶.

Notes

1. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), Section 2 (1) "モンテカルロ法の定義と実務での適用対象" (Definition of the Monte Carlo method and its uses in practice) (page 54) (in Japanese). ↩ ↩² ↩³

2. Takashi Kawakami, "ヘッジ取引におけるデリバティブ信用評価調整の影響についての考察" (A study of the effects of credit valuation adjustments of derivatives on hedging transactions), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (April 2018), the paragraph on calculating expected exposure on page 126 (page 8 of the PDF) (in Japanese). ↩ ↩²

3. Values from the same calculation as this page's simulator, using the function computeMonteCarlo in src/lib/simulators/monte-carlo.ts. The standard deviation $s$ is the standard error × $\sqrt{N}$, which is 13.41 with 10,000 paths. The analytic price is the value from the function in src/lib/simulators/black-scholes.ts (8.4333). ↩ ↩² ↩³ ↩⁴ ↩⁵ ↩⁶ ↩⁷ ↩⁸

4. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), Example 2 (pages 4–5). The price is the expected discounted payoff under the risk-neutral probability $Q$, under which the stock price is assumed to follow geometric Brownian motion with drift $r$. ↩

5. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), steps 1 to 5 at the start of Section 3 (page 57) (in Japanese). ↩

6. Martin Haugh, "Generating Random Variables and Stochastic Processes" (Columbia University, 2017), Section 7.2 (pages 18–19). The source's formula is $S_{t+\Delta t} = S_t e^{(\mu - \sigma^2/2)\Delta t + \sigma(B_{t+\Delta t} - B_t)}$, where $B$ is Brownian motion. The text replaces $\mu$ with $r$ and the increment $B_{t+\Delta t} - B_t$ with $\sqrt{\Delta t}\, Z$. ↩ ↩²

7. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), Appendix 3 "派生商品のプライシングを行う際の数値計算手法" (Numerical methods for pricing derivatives) (page 85) (in Japanese). ↩ ↩² ↩³

8. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), the paragraph after equation (A-5) in Appendix 1 (page 78) and note 14 on page 55 (in Japanese). ↩

9. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), Section 1 (page 1). The source writes the confidence interval as $\hat{\theta}_n \pm z_{1-\alpha/2}\,\hat{\sigma}_n/\sqrt{n}$. The text sets $\alpha = 0.05$ and $z_{0.975} = 1.96$. ↩ ↩²

10. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), equations (A-4) and (A-5) in Appendix 1 (page 78) (in Japanese). ↩

11. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), note 5 in Section 1 (page 52) (in Japanese). ↩

12. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), Section 3 (1) "乱数の周期と Excel の組み込み乱数" (The period of random numbers and Excel's built-in random numbers) and note 18 (page 57) (in Japanese). ↩ ↩²

13. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), Section 3 (1) and (2) and "図表 2" (Figure 2) (pages 57–58) (in Japanese). ↩

14. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), note 20 on page 58 (in Japanese). ↩

15. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), the opening of Appendix 2 "分散減少法とモーメント調整法" (Variance reduction and moment matching methods) and "図表 A-1" (Figure A-1) (page 80), and "①対称変量法" (1. Antithetic variates) (pages 81–82) (in Japanese). ↩ ↩²

16. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), Section 3, Theorem 1 in Section 3.1, and Section 3.3 (pages 7–9). ↩

17. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), the paragraph after Example 8 in Section 3 (page 10). ↩

18. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), control variates in "図表 A-1" (Figure A-1) (page 81) (in Japanese). The paper writes $Z = X - \alpha(Y - E[Y])$. The text's $Y$, $W$, and $c$ correspond to the paper's $X$, $Y$, and $-\alpha$. ↩

19. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), Section 2.1 (page 3). The text writes the source's $Z$ as $W$. The source says that when the covariance and variance are unknown, they are first estimated with a small pilot run. ↩

20. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), Example 2 (pages 4–5). ↩

21. Martin Haugh, "Simulation Efficiency and an Introduction to Variance Reduction Methods" (Columbia University, 2017), Section 1.1 (page 2). The source compares the efficiency of methods by the product of the variance and the computational effort per path. ↩

22. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), item ④ in Section 1 (page 53) (in Japanese). ↩

23. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), equations (A-11) to (A-13) in Appendix 2 (page 82) (in Japanese). ↩ ↩²

24. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), item ③ in Appendix 3 (page 86) (in Japanese). ↩

25. Mark Broadie, Paul Glasserman, and Steven Kou, "A Continuity Correction for Discrete Barrier Options" (Mathematical Finance, October 1997), Section 1.1 (pages 325–326). The corrected analytic price is calculated with the barrier set to 88.53 by the same paper's correction (multiplying the barrier by $e^{\pm 0.5826\,\sigma\sqrt{\Delta t}}$). ↩

26. Tatsuya Ishikawa and Yoshihiko Uchida, "モンテカルロ法によるプライシングとリスク量の算出について" (Pricing and risk measurement with the Monte Carlo method), Kinyu Kenkyu, Institute for Monetary and Economic Studies, Bank of Japan (June 2002), Section 2 (2) and (3) (pages 55–56) (in Japanese). The larger error of the 99.5th percentile is the result of the example in "図表 1" (Figure 1), which calculates the expected value and the 99.5th percentile of a normally distributed value. ↩