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

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!

• Sure, thanks for your remark. I I'll clean it up. However, I think I made some reasoning error: using the distribution of r_T under Q instead of its distribution under Q_T. But then, I don't know if I should write the payoff in function of B_T or r_T. Sep 11, 2022 at 11:56
• Thx for commenting. I recommend the link instead. It is more detailed and easier to understand Sep 11, 2022 at 16:16

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.

• Excellent, thank you for the detailed answer. Actually, I found the same the result in my other Monte carlo functions. So now, I know that my analytical result was wrong since the beginning :( Sep 15, 2022 at 10:53
• When compiling you code, I've got NameError: name 'E_r' is not defined. Could I take look on your function E_r and Var_ r to be sure ? please. Correct me if I am wrong because it's surprising that you have two exact values. It's supposed to be a Monte Carlo and it implies some errors, right ? Sep 15, 2022 at 13:43
• Thanks for the suggestion, Alper. I now show what packages I am using. Apologies, had forgot to insert the Python implementation of E(r) and Var(r). They are now included in the answer. I am not sure I understand your second question. I have not explicitly calculated the standard error of the Monte-Carlo simulation, but the number is of course uncertain. If you comment out the np.random.seed(2) you will get different numbers every time you run the simulation as the random numbers in the simulation will change. Sep 16, 2022 at 18:41
• @Feynman_kac The answer coming from the QuantLib (C++) implementation is the same as that from the first implementation with Python, if that is what you mean, because I understand both are based on the closed form. Sep 16, 2022 at 21:31
• Thanks, Alper. mc() inputs are now corrected. Both the discountBondOption C++ implementation and ZBO Python implementation are based on the closed-form formula, so they should of course be the same, as Alper correctly has stated. Sep 17, 2022 at 14:39