Monte Carlo Methods in Option Pricing

On this page, I will give the reader a brief idea about Monte Carlo technique and some of its applications in finance.

Generating sample paths
 Assume that the price of an asset, $S_t$, follows a geometric Brownian motion with a drift, $\mu$, and a standard deviation, $\sigma$, as follows

$$\mathrm{d} S_t=\mu S_t\mathrm{d}t+ \sigma S_t\mathrm{d}B_t,$$

where $B_t$ is a standard Brownian motion. We can then use Ito's lemma to get
$$S_{t_{i+1}}=S_{t_i} e^{(\mu - \frac{\sigma^2}{2})~\mathrm{d}t + \sigma \sqrt{\mathrm{d} t} Z_i},$$
where $Z_i$ is a random variable that is standard normally distributed, i.e. $Z_i \sim N(0,1)$. To generate a path over time interval [0,T] (i.e. $S=\{S_{t_0}, S_{t_1},\dots, S_{t_n}\}$) with an increment of $n=\frac{T}{\Delta t}$, we need to generate a sequence of $n$ random numbers generated from a standard normal distribution (i.e. $\{Z_0,Z_1,\dots ,Z_{n-1}\}$). We then calculate $S_{t_i}, ~\forall i\in \{1,2,\dots,n\}$. Similarly, we can generate $m$ paths and denote the $k$-th path by $S^{[k]}=\{S^{[k]}_{t_0}, S^{[k]}_{t_1},\dots, S^{[k]}_{t_n}\}$. Figure 1 shows 10 paths generated by parameters $S_0=36, ~\mu =0.06, ~\sigma=0.1$ and $T=1$ for $n=50$ time steps.

We can then calculate the discounted cash flow at time $t=0$ of the option (take a European call option as an example) for the $k$-path as $$c^{[k]}_{t_0}= c^{[k]}_{t_n} \times e^{-rT} =\mathrm{E} \left[ \mathrm{max} (S^{[k]}_{t_n}-K,0)\right] \times e^{-rT}.$$ Then estimate $c$ by averaging $c^{[k]}_{t_0}$ over all $k$ paths. Note that the price of a standard European option is determined by the price of the termination, $S_{t_n}$, and is independent of the path of $S_t$. Therefore, we are able to sample the value of $S_{t_n}$ directly from $S_{t_0}$ by setting $\Delta t=T$. Asian options are path-dependent and their pay-offs depend on the value of the underlying asset at intermediate time steps.



Discretization methods
We obtained an exact solution of the Geometric Brownian motion, which can be simulated using MC technique. The exact solutions cannot always be obtained, and in such cases, a numerical approximation to the solution can be used. When using a numerical approximation, a numerical error is introduced by the discretization, in addition to the Monte Carlo error of the variance of the estimator. We try to overlap the additional error by using small time steps (i.e. a large $n$) and a good numerical scheme. Euler Scheme is the simplest numerical scheme that is covered on this page. Milstein scheme is also covered and is considered an alternative to Euler scheme.

Assume that $S_t$ satisfies the following SDE $$\mathrm{d} S_t=\mu S_t\mathrm{d}t+ \sigma S_t\mathrm{d}B_t.$$ We can then apply Euler Scheme to the above equation to get $$S_{t_{i+1}}-S_{t_i}= \mu S_{t_i} \Delta t + \sigma S_{t_i} \sqrt{\Delta t} Z_i ~~~~~ , i=\{0,1,\dots,n-1\}.$$ where $Z_i$ is a random number generated from a standard normal distribution. We can also discretize the option value $c$ by approximating $e^{-r\Delta t}$ by $\approx 1 - r\Delta t$, to get $$c_{t_i}=c_{t_{i+1}} (1-r_{t_{i}} \Delta t) ~~~~~ , i=\{n-1, n-2, \dots,0\},$$
We restrict the discretization to a grid with a constant time step $\Delta t= \frac{T}{n}$ such that $t_i= i \times \Delta t= i \times \frac{T}{n}$. We can then calculate the price of a European call option at time 0 by simulating $m$ paths of the price of the underlying security $S$ and do the following for every path $k \in [1,m]$:

1. Simulate $S^{[k]}_{t_n}$. (We will need to simulate $S^{[k]}_{t_i}$ for all $i \in [0,n-1]$, when you price a path-dependent option or when the price of the underlying security is discretized).
2. Calculate the payoff of the European call option $c^{[k]}_T=\max(S^{[k]}_T-E,0)$.
3. Calculate the discounted value of the payoff, i.e. $c^{[k]}_0$.

Then take the average of the discounted cash flows $\overline{c}^{[k]}_0$ over all paths. Note that risk-free interest rate $r$ can follow a stochastic process and violate an assumption in the Black-Scholes model of a constant $r$.

In addition to the statistical error of the estimator, we have the numerical error that is introduced by discretizing the estimator. The numerical error of a call option price (as an example) is represented by the bias of the estimator $\hat{c}$ as: $$\mathrm{E}[\hat{c}]-c=\overline{c}-c=\frac{1}{m} \sum_{k=1}^m c^{[k]}_{0}-c,$$ where $\overline{c}$ is the average of $c^{[k]}_0$ over $m$ paths. This error can be reduced by discretizing over a smaller time interval, $\Delta t$, or by using a more efficient method for discretization. The mean square error $MSE(\hat{c})$ is used as an error measure that balances between the bias and the variance of the estimator. If $\hat{c}$ is an estimator of $c$, then the mean square error equals to the sum of the bias squared and the variance of $\hat{c}$ as follows $$MSE(\hat{c})=\mathrm{E}\left[ (\hat{c}-c)^2 \right]=\mathrm{E} \left[(\hat{c}-\mathrm{E}[\hat{c}])^2 \right]+ \left(\mathrm{E}[\hat{c}]-c \right)^2\\= \mathrm{Var}(\hat{c})+\mathrm{Bias}^2(\hat{c}).$$
Exercise
Heston model violates one of the assumptions of the Black-Scholes model. Heston model does not assume that the volatility is constant, but rather that it follows a stochastic process that is correlated with the share price process. Discretize the following processes introduced by Heston $$\mathrm{d}S_t= \mu S_t\mathrm{d}t + \sigma_t S_t\mathrm{d}B_t^{(1)}\\ \label{Var_Var}\mathrm{d}V_t= k (\theta -V_t)\mathrm{d}t + \alpha \sigma_t\mathrm{d}B_t^{(2)},\\ $$ where $B_t^1$ and $B_t^2$ are standard Brownian motions that are correlated with correlation coefficient $\rho$ that is $\mathrm{d}B_t^1~\mathrm{d}B_t^2 = \rho\mathrm{d}t$. The volatility of the volatility $\alpha$, the mean revision rate $k$, the drift $\mu$ and the long run variance $\theta$ are positive constants. $V_t=\sigma_t^2$ is the variance at time $t$ and $S_t$ is the price of the underlying asset at time $t$.

Milstein Scheme
Milstein has added an additional term to increase the order of the accuracy of the stochastic term in Euler scheme. The derivation of Milstein scheme is based on Ito's lemma. Milstein scheme can be applied to a stochastic process as defined previously, to get the following $$X_{t_{i+1}}-X_{t_i}= b(t,\mathrm{X}_t) \Delta t + v(t,\mathrm{X}_t) \sqrt{\Delta t} Z_i + \frac{1}{2} \frac{\partial v}{\partial\mathrm{X}_t}(t,\mathrm{X}_t) v(t,\mathrm{X}_t) \Delta t (Z_i^2-1) ~~~~~ , i=\{0,1,\dots,n-1\},$$ where each $Z_i$ is a random number that has been taken from the standard normal distribution $N(0,1)$. Milstein scheme can also be applied to the Geometric Brownian motion to get the following $$S_{t_{i+1}}-S_{t_i}= \mu S_{t_i} \Delta t + \sigma S_{t_i} \sqrt{\Delta t} Z_i + \frac{1}{2} \sigma^2 S_{t_i}\Delta t (Z_i^2-1) ~~~~~ , i=\{0,1,\dots,n-1\}.$$ Note that a new term has been introduced in Milstein scheme in comparison with Euler scheme. Milstein scheme has expanded the diffusion term from $O(\sqrt{\Delta t})$ to $O(\Delta t)$. This term is a stochastic term with an expected value of zero and it is uncorrelated with Euler terms.

Exercise
Show that the stochastic term in Milstein scheme has an expected value of zero and that it is uncorrelated with the Euler terms.

Solution.
We know that $Z_{i}$ follow s a normal distribution with zero mean and a standard deviation of 1. We need to calculate the expectations of $Z_{i}^2$ and $Z_{i}^3$ since we will use them in the solution. $\mathrm{E}(Z_{i}^2)=1$ since $\text{Var}(Z_{i})=\mathrm{E}(Z_{i}^2)-[\mathrm{E}(Z_{i})]^2$. $\mathrm{E}(Z_{i}^3)=\int^{\infty}_{-\infty} z^3 f(z)\mathrm{d}z=0$, since we are integrating an odd function over $(-\infty,\infty)$ where the probability density function $f(z)=\frac{1}{\sqrt{2 \pi}} e^{-\frac{z^2}{2}}$. Therefore the expected value of the diffusion term is $$\mathrm{E} \left ( \frac{1}{2} \sigma^2 S_{t_i}\Delta t (Z_i^2-1) \right )=\frac{1}{2}\sigma^2 \Delta t\mathrm{E}( S_{t_i}) \left (\mathrm{E}(Z_i^2)-1) \right)=0.$$
Set $$A=\mu S_{t_i} \Delta t + \sigma S_{t_i} \sqrt{\Delta t} Z_i,$$ $$B=\frac{1}{2} \sigma^2 S_{t_i} \Delta t (Z_i^2-1),$$
then $$\mathrm{Cov}(A,B)=\mathrm{E}[A~B]-\mathrm{E}[A]~\mathrm{E}[B]$$ $$\label{Cov_AB}=\mathrm{E} \left[ \left(\mu S_{t_i} \Delta t + \sigma S_{t_i} \sqrt{\Delta t} Z_i \right) \frac{1}{2} \sigma^2 S_{t_i}\Delta t (Z_i^2-1) \right]-\mathrm{E}[A]~\mathrm{E}[B].$$ The covariance is therefore $$\mathrm{Cov}(A,B)=\mathrm{E} \left[ \left(\mu S_{t_i} \Delta t + \sigma S_{t_i} \sqrt{\Delta t} Z_i \right) \frac{1}{2} \sigma^2 S_{t_i} \Delta t (Z_i^2-1) \right]\\=\frac{1}{2} \sigma^2 \Delta t \left( \mu \Delta t\mathrm{E}[S_{t_i}^2]\mathrm{E}[Z_i^2-1]+ \sigma \sqrt{\Delta t}\mathrm{E}[S_{t_i}^2]\mathrm{E}[Z_i^3-Z_i] \right).$$ Use $\mathrm{E}(Z_{i})=0$, $\mathrm{E}(Z_{i}^2)=1$ and $\mathrm{E}(Z_{i}^3)=0$, and that $S_{t_i}$ and $Z_i$ are independent random variables to get $$\mathrm{Cov}(A,B)=0.$$
Graphs

Figure 1. Geometric Brownian motion

Figure 2. Geometric Brownian motion, a comparison between exact and numerical approximations

Figure 3. Mean reverting process

Figure 4. Geometric Brownian motion with an increasing risk-free interest rate

---

References

1. Boyle P.P., (1977), Options: A Monte Carlo Approach, Journal of Financial Economics 4 (1977) 323-338.
   
2. Boyle P.P., Broadie M. and Glasserman P., (1997), Monte Carlo methods for security pricing, Journal of Economic Dynamics and Control, 21 (1997) 1267-1321
   
3. Heston H., (1993), A Closed-form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options, The Review of Financial Studies, 6(2) P 327-343.
   
4. Robert C.P and Casella G., (2004), Monte Carlo Statistical Methods, 2nd edition, Springer.
   
5. Glasserman P., 2004, Monte Carlo Methods in Financial Engineering, Springer.
   
6. Higham D.J., (2001), An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev., 43(3), 525-546.
   
7. Milshtein G. N., (1974), Approximate Integration of Stochastic Differential Equations, Theory of Probability and its Application 19:557-562. 267.
   
8. Uhlenbeck G. E. and Ornstein L. S., (1930) , On the Theory of the Brownian Motion, Physical Review 36, 823-841.