In this blogpost we will dive into the world of finance, and specifically option pricing. We are going to illustrate how we can use Monte Carlo simulations to price financial options, by looking at two examples. In both of these examples no explicit solution formula exists, making the Monte Carlo approach useful in pricing these options.

If you are not familiar with Monte Carlo simulations or financial options, or could use a reminder, I suggest that you take a look at the post Introduction to Monte Carlo simulations and option pricing where we introduced all these concepts. In that post we:

Introduced derivative assets and financial options, especially the European call option
Looked at our first model for stocks, namely the standard Black & Scholes model
Introduced the concept of Monte Carlo simulations
Looked at our first pricing example

The examples in this post are definitely more challenging mathematically, so make sure that you have a grasp on the basics from the introductory post.

Let's get started.

Arithmetic basket call in the standard Black & Scholes model

For our first example, we will use the standard Black & Scholes model, i.e. the same one used in the introduction post about Monte Carlo simulations and option pricing. Here the stock dynamics follow a Geometric Brownian Motion:

dS(t) = rS(t)dt + \sigma S(t)dW(t)

with the solution:

S(T) = S(0)e^{(r-\frac{\sigma^2}{2})T + \sigma W(T)} = e^{\log(S(0)) + (r-\frac{\sigma^2}{2})T + \sigma W(T)}

where $W(T)\sim N(0,T)$ , so $S(T)$ has a log-normal distribution. The contract we will look at now, is very similar to a European call option, but instead of having just one underlying stock, we take the arithmetic average of several stocks ( $n$ of them). This contract is called an arithmetic basket call option, and has the following payoff at time $T$ :

\Phi(S_1(T),...,S_n(T)) = \max\left((\frac{1}{n}\sum_{i=1}^n S_i(T)) - K, 0\right)

Even if we know the distribution of every stock $S_i$ , we do not know the distribution of the arithmetic average of $n$ of them (the log-normality is not preserved under summation). We can therefore not find an explicit formula like the Black & Scholes formula for European call options in this case for the option price:

\Pi(0) = e^{-rT}\mathbb{E}\left[\max\left((\frac{1}{n}\sum_{i=1}^n S_i(T)) - K, 0\right)\right]

Our solution will be to use Monte Carlo simulations. For simulation purposes, we will use that the log stock prices $\log(S_1(T)),...,\log(S_n(T))$ follow a multivariate normal distribution. For simplicity, we assume that all stocks have the same initial value $S(0)=s$ , same volatility $\sigma$ and that the correlation $\rho$ between all the stocks are the same. The log stock prices will then follow a multivariate normal distribution with the following expected value vector $\mathbf{\mu}$ and covariance matrix $\Sigma$ :

\boldsymbol{\mu} = \begin{bmatrix} \log(s) + (r-\frac{\sigma^2}{2})T \\ \vdots \\ \log(s) + (r-\frac{\sigma^2}{2})T \end{bmatrix} \quad \text{and} \quad \boldsymbol{\Sigma} = \begin{bmatrix} \sigma^2 & \sigma^2\rho & \dots & \sigma^2\rho \\ \sigma^2\rho & \sigma^2 & \dots & \sigma^2\rho \\ \vdots & \vdots & \ddots & \vdots \\ \sigma^2\rho & \sigma^2\rho & \dots & \sigma^2 \end{bmatrix}T

Note that $\boldsymbol{\mu}$ is an $n\times 1$ -vector and $\boldsymbol{\Sigma}$ is an $n\times n$ -matrix. This is useful because we can now simulate from this multivariate normal distribution and then transform the values using the exponential function to get simulations of the stocks. We thus have the following Monte Carlo sampler for estimating the price of the arithmetic basket call option:

\begin{aligned} &\textbf{for } i=1 \textbf{ to } N \\ &\quad\quad \text{draw } [\log(S_1(T)),...,\log(S_n(T))]_i \sim MVN(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \\ &\quad\quad \text{transform to } [S_1(T),...,S_n(T)]_i \text{ with the exponential function} \\ &\quad\quad \text{compute arithmeric average } (\frac{1}{n}\sum_{j=1}^n S_j(T))_i \\ &\textbf{end} \\ &\Pi(0)_N^{MC} = \frac{1}{N}\sum_{i=1}^N e^{-rT}\cdot\max\left((\frac{1}{n}\sum_{j=1}^n S_j(T))_i - K, 0\right) \end{aligned}

So, assuming that we have $n=10$ underlying stocks, that $s=100$ , the strike price $K=90$ , the risk free interest rate $r=1\%$ , the stock volatility $\sigma=0.2$ , the correlation $rho=0.4$ and that the time to maturity $T=1$ year, we can implement the following Monte Carlo sampler in Python:

import numpy as np from scipy.stats import norm n = 10 # number of underlying stocks in the arithmetic basket call s = 100 # initial stock prices K = 90 # strike price sigma = 0.2 # volatility of the stocks rho = 0.4 # correlations between the stocks r = 0.01 # risk free interest rate (continuously compounded) T = 1 # time to maturity (years) mu = (np.log(s) + (r - 0.5*sigma**2)*T)*np.ones(n) # expected value vector of the log stock prices # covariance matrix for the log stock prices cov_mat = (rho*sigma**2)*np.ones((n,n)) cov_mat = cov_mat - np.diag(np.diag(cov_mat)) + np.diag((sigma**2)*np.ones(n)) cov_mat = cov_mat*T N = 100000 # size of the Monte Carlo sample # generate Monte Carlo sample mc = np.random.multivariate_normal(mu, cov_mat, N) # log stock prices mc = np.exp(mc) # stock prices mc = np.mean(mc, 1) # arithmeric mean of the stock prices mc = np.exp(-r*T)*np.maximum(mc - K, 0) # option prices price = np.mean(mc) # estimated price # confidence interval price_std = np.std(mc) / np.sqrt(N) # standard deviation of the price estimator lo = price - norm.ppf(0.975)*price_std # lower bound of confidence interval for the price hi = price + norm.ppf(0.975)*price_std # upper bound of confidence interval for the price

Running a simulation (with $N=100000$ ) gives us the following price estimate and 95% confidence interval for the Monte Carlo estimated price:

\begin{aligned} \Pi&(0)_N^{MC} = 12.308 \\ I_{\Pi(0)} &= (12.236, 12.381) \end{aligned}

As mentioned in the first post, we will get different values if we were to repeat the simulations. So if you try and implement this by yourself you would not get exactly the same values, but they should however be very similar.

European call option in the Heston model

When comparing actual options prices on various financial markets with the prices given by the standard Black & Scholes model, it becomes clear that they generally do not match very well. The standard Black & Scholes model is just not complex enough to explain the stock's price process. There are several adjustments one could make to get a better model. One major drawback of the standard Black & Scholes model is that the stock volatility is constant and deterministic, while the actual stock volatility is a stochastic process just like the stock's price process. We will therefore now look at a model with stochastic stock volatility, namely the Heston model.

In the Heston model, we have dynamics given as stochastic differential equations, for both the stock and the volatility. These dynamics looks as follows:

\begin{aligned} &dS(t) = rS(t)dt + \sqrt{V(t)}S(t)(\sqrt{1-\rho^2}dW^{(1)}(t) + \rho dW^{(2)}(t)) \\ &dV(t) = \kappa(\theta - V(t))dt + \sigma\sqrt{V(t)}dW^{(2)}(t) \end{aligned}

where $S(t)$ is the stock's price process as before, and $V(t)$ is the stock's volatility process. $W^{(1)}(t)$ and $W^{(2)}(t)$ are here two independent standard Brownian Motions. Also, $r$ is the risk free interest rate, and $\rho$ , $\kappa$ , $\theta$ and $\sigma$ are model parameters.

We return to the European call option that we saw in the introductory post. As you might remember, this means that we want to find the option price:

\Pi(0) = e^{-rT}\mathbb{E}[\max(S(T)-K,0)]

In the standard Black & Scholes model we could find an explicit formula for this price. In the Heston model however, this is not the case. The reason is that we cannot find a closed form solution to the system of stochastic differential equations for $S(t)$ and $V(t)$ above. We will therefore need numerical methods and Monte Carlo simulations to generate a sample of stock prices $S(T)$ that we can use to estimate the option price.

So we cannot solve these stochastic differential equations analytically, and will instead rely on numerical methods. As always when dealing with continuous functions in numerical analysis, we are going to use some kind of discretization scheme to solve these equations. This means that instead of dealing with this as a continuous problem over the interval $[0,T]$ , we are reducing this continuous interval to discrete time points $\{t_0,t_1,...,t_{n-1},t_n\}$ where we have that $t_0=0$ and $t_n=T$ . We will then, using some discretization scheme, calculating the values of $S(t)$ and $V(t)$ at these time points only.

We start with the differential equation for $S(t)$ . We will here use something called the Euler Scheme. For a stochastic differential equation on the form:

dS(t) = g(S(t))dt + h(S(t))dW(t)

where $g$ and $h$ are functions, the Euler scheme is:

S_i = S_{i-1} + g(S_{i-1})\Delta + h(S_{i-1})\Delta W

for $i=1,...,n$ , where $\Delta=t_i-t_{i-1}$ is the size of each time step and $\Delta W=W_{t_i}-W_{t_{i-1}} \sim N(0,\Delta)$ is the increment of the standard Brownian Motion. Comparing with the equation for $S(t)$ , the Euler scheme in this case becomes:

S_i = S_{i-1} + rS_{i-1}\Delta + \sqrt{V_{i-1}}S_{i-1}(\sqrt{1-\rho^2}\Delta W^{(1)} + \rho\Delta W^{(2)})

for $i=1,...,n$ . Note that the stochastic part here are the $\Delta W$ s, so the Monte Carlo aspect will be to simulate these normal random variables. Doing this scheme $N$ times, we will have a Monte Carlo sample of size $N$ of the stock prices $S(T)=S_n$ . We must of course simulate the process for $V(t)$ as well. To ensure that the $V_i$ s take on only positive values, we add another term to the discretization scheme. We will not go into details here but this is called the Milstein scheme (this extra term can be seen in the code in the end). If you are interested you can check out the Wikipedia page.

So, we have the following Monte Carlo sampler for estimating the price of the European call option:

\begin{aligned} &\textbf{for } i=1 \textbf{ to } N \\ &\quad\quad \textbf{for } j=1 \textbf{ to } n \\ &\quad\quad\quad\quad \text{draw } \Delta W_i^{(1)}\sim N(0,\Delta) \text{ and } \Delta W_i^{(2)}\sim N(0,\Delta) \\ &\quad\quad\quad\quad \text{compute } V_j \text{ and } S_j \text{ with respective discretization scheme} \\ &\quad\quad \textbf{end} \\ &\quad\quad \text{let } S(T)_i = S_n \\ &\textbf{end} \\ &\Pi(0)_N^{MC} = \frac{1}{N}\sum_{i=1}^N e^{-rT}\cdot\max(S(T)_i - K, 0) \end{aligned}

We assume that the initial stock price is $S(0)=100$ , the initial volatility is $V(0)=0.16$ , the strike price $K=90$ , the risk free interest rate $r=1\%$ and the time to maturity is $T=1$ year. Also, we assume that our model parameters have the values $\kappa=10$ , $\theta=0.16$ , $\sigma=0.1$ and $\rho=-0.8$ . Then, using $n=252$ time steps in the discretization (typically 252 trading days on a year, so daily time steps as $T=1$ year), we can implement the following Monte Carlo sampler in Python:

Running a simulation (with $N=100000$ ) gives us the following price estimate and 95% confidence interval for the Monte Carlo estimated price:

\begin{aligned} \Pi&(0)_N^{MC} = 21.071 \\ I_{\Pi(0)} &= (20.866, 21.276) \end{aligned}

Conclusion

In this blogpost, we made use of the concepts introduced in Introduction to Monte Carlo simulations and option pricing to price two quite tricky financial options using Monte Carlo simulations: the arithmeric basket call option in the standard Black & Scholes model, and the European call option in the Heston model. To summarize the highlights, we have:

Introduced the arithmetic basket call option
Looked at an example of how we can simulate using multivariate distributions
Introduced the Heston model for stocks
Introduced discretization and the Euler scheme for solving stochastic differential equations
Seen two examples of Monte Carlo samplers for estimating expected values

Published on February 11th 2020

Last updated on March 1st 2023, 10:53