Problem with pricing a call option using the Monte Carlo Vasicek model

Question

I am trying to price a call option on a zero coupon under the Vasicek Model using Monte Carlo method:

$$C_0 = B(0,\theta) \ \mathbb{E}^{\mathbb{Q}_T}[(B(\theta,T)-K)^{+}]$$

The problem is that the code I have written for this is yielding a price too different from the expected answer which is 0.095 based on the closed form formulas.

What am I possibly doing wrong?

The underlying assumptions and equations I am using:

Let $(Ω, (\mathcal{F}_t)_t, \mathbb{Q})$ such that $\mathbb{Q}$ is risk neutral probability.

Vasicek model: $dr_t = a(b - r_t ) dt + \sigma dW_t^{\mathbb{Q}}$

Price of zero coupon bond: $B(T,S) = \tau(T,S) \ exp(-\eta(T,S)\ r_t)$

Since $r_t$ is gaussian under $\mathbb{Q}_t$ (forward measure): $$r(t) \sim \mathcal{N}(r_0e^{-at} + b(1-e^{-at}) + \frac{\sigma(1 - e^{-at})}{a^2} - \frac{\sigma(1 - e^{-2at})}{2a^2} \ , \ \sigma^{2} \frac{1 - e^{-2at}}{2a}) $$

And my code based on the above: (one can simply copy and paste and then run the code below assuming the requisite external Python libraries are already installed.)

import numpy as np
import pylab as plt
from random import gauss
from math import exp, sqrt, log
from scipy.stats import norm
a = 0.2
b = 0.1
sigma = 0.02
r = 0.1
t = 0
theta = 0.5
T = 1
K = 0.9

#Distribution of rt
r_0 = r
m_rt = r * exp(-a*theta) + b * (1 - exp(-a*theta)) + sigma * ((1 - exp(-a*theta)) / a**2 ) - sigma * ((1 - exp(-2*a*theta))/2*(a**2))
v_rt = (sigma**2) * (1 - exp(2*a*theta))/ 2*a


eta = (1 - exp(-a*(T-theta)))/a
tau = np.exp( (b - (sigma**2)/(2*(a**2))) * (eta - T + theta) - (0.25 * (sigma**2) * (eta**2) / a) )


def VasicekZCPrice(t,T,a,b,sigma,r):
    B = ( 1- exp( -a*(T-t) ) )/ a
    A = exp( ( (b-(0.5*sigma**2)) / a**2 )*( B - T + t ) - ( 0.25*sigma**2 / a ) * B**2)
    ZC = A*exp(-B*r)
    return ZC

def monte_carlo_call(theta, T, a, b, sigma, r, K, nbSimul = 10000, steps = 100):
    dt = (T-theta)/steps
    r__T = r * np.exp(-a*dt) +  b* (1-np.exp(-a*dt)) + (sigma*(1-exp(-a*dt))/a**2) - sigma * ((1-np.exp(-2*a*dt))/2*a**2) + sigma*np.sqrt((1-np.exp(-2*a*dt))/2*a) * np.random.normal(size=(steps,nbSimul))
    B__T = tau * np.exp(-eta*r__T)
    payoffs = np.maximum( B__T - K, 0 )
    option_price = VasicekZCPrice(t, T,a,b,sigma,r) * np.mean(payoffs)
    return option_price

monte_carlo_call(theta, T, a, b, sigma, r,K ,  nbSimul = 10000, steps = 100)

which gives me:

0.03686737433679907

An expanded version of the code and further explanations are accessible through this link on Google Colab.

Thanks very much!

Answer 1

To make sure that I understand the problem: you are trying to price a call option expiring at time 0.5, which will exercise into a unit notional zero-coupon bond with a maturity of 1.0 at a strike (price) of 0.9.

I will follow the notation from "Interest Rate Models – Theory and Practice" by Brigo & Mercurio (2006) which is slightly different than your notation.

From page 58 in the mentioned book, we get the dynamics of the Vasicek model as:

$$ dr(t)=k[\theta-r(t)]dt+\sigma dW(t) $$

And we are choosing the parameters:

k = 0.2
theta = 0.1
sigma = 0.02
r_0 = 0.1

time_now = 0
expiry = 0.5
maturity = 1
strike = 0.9

In QuantLib (C++), the function of interest is already implemented. Unfortunately, this function is not exposed in the Python binding, but we can use it as a benchmark:

Vasicek vasicek(0.1,0.2,0.1,0.02);
std::cout << std::setprecision(16);
std::cout << vasicek.discountBondOption(Option::Type::Call,0.9,0.5,1.0) << std::endl;

It returns a value of 0.0487763759011719. Let us see whether we can get this value with Python.

From page 59 in the same book as before, we can calculate the price of a zero-coupon bond using:

$$ P(t,T)=A(t,T)e^{-B(t,T)r(t)} $$ where $$ A(t,T)=\exp\left\{\left(\theta-\frac{\sigma^{2}}{2k^{2}}\right)\left[B(t,T)-T+t\right]-\frac{\sigma^{2}}{4k}B(t,T)^{2}\right\} $$ and $$ B(t,T)=\frac{1}{k}\left[1-e^{-k(T-t)}\right] $$

To calculate these functions in Python, we import the following packages:

import numpy as np
from math import exp, sqrt, log
from scipy.stats import norm

and implement the functions as follows:

def B(t,T):
    return (1/k)*(1-exp(-k*(T-t)))
def A(t,T):
    return exp((theta-pow(sigma,2)/(2*pow(k,2)))*(B(t,T)-T+t)-(pow(sigma,2)/(4*k))*pow(B(t,T),2))
def P(t,T,r_t):
    return A(t,T)*exp(-B(t,T)*r_t)

From page 60, we have the zero-coupon bond call option formula:

$$ ZBO(t,T,S,X)=P(t,S)\Phi(h)-XP(t,T)\Phi(h-\sigma_{p}) $$

where

$$ \sigma_{p}=\sigma\sqrt{\frac{1-e^{-2k(T-t)}}{2k}}B(T,S) $$

$$ h=\frac{1}{\sigma_{p}}\ln\frac{P(t,S)}{P(t,T)X}+\frac{\sigma_{p}}{2} $$

which we implement as:

def ZBO(t,T,S,X,r_t):
    sigma_p=sigma*sqrt((1-exp(-2*k*(T-t)))/(2*k))*B(T,S)
    h=(1/sigma_p)*log(P(t,S,r_t)/(P(t,T,r_t)*X))+sigma_p/2
    return P(t,S,r_t)*norm.cdf(h, 0, 1)-X*P(t,T,r_t)*norm.cdf(h-sigma_p, 0, 1)

The value of the option calculated through ZBO(time_now,expiry,maturity,strike,r_0) comes out as 0.0487763759011719.

This is the same as the value calculated using QuantLib. I am not sure where you got your closed form formula from, but it seems a bit off.

Implementing a Monte Carlo simulation is a bit tricky because we are interested in not only the simulated $r(t)$ from $0$ to $T$, but also the integrated interest rate $\int_0^Tr(t)dt$. This is because we are trying to evaluate:

$$ e^{-\int_0^Tr(t)dt}\left[P(T,S)-X\right]^+ $$

While the value of the zero-coupon bond within the parenthesis only depends on the terminal interest rate $r(T)$, the discount factor at the beginning depends on the instantaneous spot rate, a.k.a. the short rate.

Now, consider the following approximation:

$$ \int_0^Tr(t)dt\approx\sum_{i=0}^{n}r(i\cdot \Delta t)\cdot \Delta t $$

This means we will divide the time into $n$ pieces and use the interest rate at the beginning of each discrete period to approximate the integrated rate.

First and foremost, we need to simulate the rate itself. From page 59, we have the distribution of the short rate under the risk-neutral measure ($\mathbb{Q}$):

$$ \mathbb{E}\left[ r(t) \middle| \mathcal{F}_{s} \right] =r(s)e^{-k(t-s)}+\theta\left(1-e^{-k(t-s)}\right) $$

$$ \text{Var}\left[ r(t) \middle| \mathcal{F}_{s} \right] =\frac{\sigma^{2}}{2k}\left[1-e^{-2k(t-s)}\right] $$

which implementation in Python is as follows:

def E_r(s,t,r_s):
    return r_s*exp(-k*(t-s))+theta*(1-exp(-k*(t-s)))
def Var_r(s,t,r_s):
    return (pow(sigma,2)/(2*k))*(1-exp(-2*k*(t-s)))

In the Monte Carlo simulation, we use equidistant steps of size $\Delta t$. So, one of the arguments we use in calling these functions is 0,dt. The interest rate input to the formula is updated at every step. These lead us to the following implementation of the Monte Carlo simulation:

np.random.seed(2)
n_steps=100
n_sim=10000

def mc_single(t,T,S,X):
    r=r_0
    dt=(T-t)/n_steps
    time=t
    interest_rate_integral=0
    for x in range(n_steps):
        interest_rate_integral=interest_rate_integral+r*dt
        r=np.random.normal(E_r(0.0,dt,r),sqrt(Var_r(0.0,dt,r)), 1)
        time=time+dt
    disc=exp(-interest_rate_integral)
    return disc*max(P(T,S,r)-X,0)

def mc(t,T,S,X):
    sum=0
    for x in range(n_sim):
        sum=sum+mc_single(t,T,S,X)
    return sum/n_sim

I am setting the seed so that it is possible to reproduce the result. The simulated value calculated through mc(time_now,expiry,maturity,strike) comes out as 0.0488357730948401.

It is not exactly the same, but close enough considering the simplicity of this algorithm.