Pricing financial options using Monte Carlo simulations - Part 2
2020-02-11 Fredrik Olsson

In the first part of this blog post about pricing financial options using Monte Carlo simulations, we did the following:

  • 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
  • Here in the second part of this blog post we will look at two, more complicated examples, where we unlike the example in the first part cannot find an explicit formula, and we thus need Monte Carlo simulations to price these option contracts.

    Arithmetic basket call in the standard Black&Scholes model

    For our first example, we will stay in the standard Black&Scholes model, i.e. 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 will we look at now, is very similar to the European call options, 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((\frac{1}{n}\sum_{i=1}^n S_i(T)) - K, 0)

    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}[\max((\frac{1}{n}\sum_{i=1}^n S_i(T)) - K, 0)]

    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] \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}

    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((\frac{1}{n}\sum_{j=1}^n S_j(T))_i - K, 0) \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 ρ=0.4\rho=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 part of this 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 had in the first part of this 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:

    S(t)=g(S(t))dt+h(S(t))dW(t) S(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)ΔW S_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 ΔWs\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 VisV_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}

    So, in this second part of this blog post we continued with pricing financial options using Monte Carlo simulations, and we looked at two examples that was a bit more complicated compared to the one in part 1. 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 our second model for stocks, namely the Heston model
  • Introduced discretization and the Euler scheme for solving stochastic differential equations
  • Seen more examples of Monte Carlo samplers for estimating expected values