Pricing financial options using Monte Carlo simulations

Fredrik Olsson


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 definetely 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+σS(t)dW(t)dS(t) = rS(t)dt + \sigma S(t)dW(t)

with the solution:

S(T)=S(0)e(rσ22)T+σW(T)=elog(S(0))+(rσ22)T+σW(T)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)N(0,T)W(T)\sim N(0,T), so S(T)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 (nn of them). This contract is called an arithmetic basket call option, and has the following payoff at time TT:

Φ(S1(T),...,Sn(T))=max((1ni=1nSi(T))K,0)\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 SiS_i, we do not know the distribution of the arithmetic average of nn 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:

Π(0)=erTE[max((1ni=1nSi(T))K,0)]\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(S1(T)),...,log(Sn(T))\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)=sS(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:

μ=[log(s)+(rσ22)Tlog(s)+(rσ22)T]andΣ=[σ2σ2ρσ2ρσ2ρσ2σ2ρσ2ρσ2ρσ2]T\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×1n\times 1-vector and Σ\boldsymbol{\Sigma} is an n×nn\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:

for i=1 to Ndraw [log(S1(T)),...,log(Sn(T))]iMVN(μ,Σ)transform to [S1(T),...,Sn(T)]i with the exponential functioncompute arithmeric average (1nj=1nSj(T))iendΠ(0)NMC=1Ni=1NerTmax((1nj=1nSj(T))iK,0)\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=10n=10 underlying stocks, that s=100s=100, the strike price K=90K=90, the risk free interest rate r=1%r=1\%, the stock volatility σ=0.2\sigma=0.2, the correlation rho=0.4rho=0.4 and that the time to maturity T=1T=1 year, we can implement the following Monte Carlo sampler in Python:

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

Π(0)NMC=12.308IΠ(0)=(12.236,12.381)\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:

dS(t)=rS(t)dt+V(t)S(t)(1ρ2dW(1)(t)+ρdW(2)(t))dV(t)=κ(θV(t))dt+σV(t)dW(2)(t)\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)S(t) is the stock's price process as before, and V(t)V(t) is the stock's volatility process. W(1)(t)W^{(1)}(t) and W(2)(t)W^{(2)}(t) are here two independent standard Brownian Motions. Also, rr 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:

Π(0)=erTE[max(S(T)K,0)]\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)S(t) and V(t)V(t) above. We will therefore need numerical methods and Monte Carlo simulations to generate a sample of stock prices S(T)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][0,T], we are reducing this continuous interval to discrete time points {t0,t1,...,tn1,tn}\{t_0,t_1,...,t_{n-1},t_n\} where we have that t0=0t_0=0 and tn=Tt_n=T. We will then, using some discretization scheme, calculating the values of S(t)S(t) and V(t)V(t) at these time points only.

We start with the differential equation for S(t)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)dS(t) = g(S(t))dt + h(S(t))dW(t)

where gg and hh are functions, the Euler scheme is:

Si=Si1+g(Si1)Δ+h(Si1)ΔWS_i = S_{i-1} + g(S_{i-1})\Delta + h(S_{i-1})\Delta W

for i=1,...,ni=1,...,n, where Δ=titi1\Delta=t_i-t_{i-1} is the size of each time step and ΔW=WtiWti1N(0,Δ)\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)S(t), the Euler scheme in this case becomes:

Si=Si1+rSi1Δ+Vi1Si1(1ρ2ΔW(1)+ρΔW(2))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,...,ni=1,...,n. Note that the stochastic part here are the ΔW\Delta Ws, so the Monte Carlo aspect will be to simulate these normal random variables. Doing this scheme NN times, we will have a Monte Carlo sample of size NN of the stock prices S(T)=SnS(T)=S_n. We must of course simulate the process for V(t)V(t) as well. To ensure that the ViV_is 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:

for i=1 to Nfor j=1 to ndraw ΔWi(1)N(0,Δ) and ΔWi(2)N(0,Δ)compute Vj and Sj with respective discretization schemeendlet S(T)i=SnendΠ(0)NMC=1Ni=1NerTmax(S(T)iK,0)\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)=100S(0)=100, the initial volatility is V(0)=0.16V(0)=0.16, the strike price K=90K=90, the risk free interest rate r=1%r=1\% and the time to maturity is T=1T=1 year. Also, we assume that our model parameters have the values κ=10\kappa=10, θ=0.16\theta=0.16, σ=0.1\sigma=0.1 and ρ=0.8\rho=-0.8. Then, using n=252n=252 time steps in the discretization (typically 252 trading days on a year, so daily time steps as T=1T=1 year), we can implement the following Monte Carlo sampler in Python:

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

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


In this plogpost, 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



Sockerkokaregatan 21

222 36 Lund


Copyright © 2021
Backtick Technologies

Backtick Technologies