Monte Carlo price of European option on ZCB under Vasicek short rate

I'm trying to replicate the analytical result from the closed form Vasicek formula for European options on zero-coupon bonds using Monte-Carlo simulation.

The interest rate paths I've simulated seem to be correct. The bond price paths that I calculate from those also seem correct (today's average of the simulated bond prices matches the analytical bond price formula under Vasicek).

My approach to calculate e.g. a call option's price is: identify the paths where the price at option maturity exceeds the strike price. For those paths, I calculate the present value of the payoff at option maturity as PV(bond price) - PV(strike price). I had already calculated the whole price paths for the bond prices, so for those I just take the value of the price path today. For the PV of the strike, I apply the same interest rate paths that generated the corresponding bond prices to discount the strike's value to today, and thus treat them as a ZCB with strike payoff at option maturity. Finally, I divide the sum of the payoffs today of those paths by the total number of paths generated.

Below you can find the code for the bond option, the simulated bond price paths under Vasicek as well as for the simulated IR paths. Many thanks in advance for you help.

def eurOptZcbVasicekMonteCarlo(simulatedRates,optionMaturity,strike,flag="Call"):

nPaths = len(simulatedRates.columns)

pricePaths = bondPrices(simulatedRates,1,True)
pricePaths = pricePaths[pricePaths.index<=optionMaturity]
strikePaths = bondPrices(simulatedRates,strike,True)
strikePaths = strikePaths[strikePaths.index<=optionMaturity]

if flag.lower()=="call":
sumPathsPV = np.sum(pricePaths.iloc[0,:]-strikePaths.iloc[0,:])

elif flag.lower()=="put":
sumPathsPV = np.sum(strikePaths.iloc[0,:]-pricePaths.iloc[0,:])

else:
Exception("Choice for 'option' can be either 'call' or 'put'. Please specify.")

optionPrice = sumPathsPV/nPaths

return optionPrice

def bondPrices(simulatedRates, faceValue, returndf=False, paths = 100):
TTM = int(simulatedRates.index[-1])
nTimes = simulatedRates.shape[0]
simulatedPrices = np.zeros_like(simulatedRates)
step = TTM/(nTimes-1)
simulatedPrices[-1]=faceValue
for time in np.arange(2,nTimes+1,1):
simulatedPrices[-time]=simulatedPrices[-time+1]*np.exp(-simulatedRates.iloc[-time+1]*step)
simulatedPrices=pd.DataFrame(simulatedPrices,index=simulatedRates.index)

if returndf == True:
return simulatedPrices
else:
simulatedPrices.iloc[:,:paths].plot(figsize=(15,12),legend=False,color='xkcd:sky blue', alpha = 0.5)
simulatedPrices.mean(axis=1).plot()
simulatedPrices.quantile(q=0.025,axis=1).plot()
simulatedPrices.quantile(q=0.975,axis=1).plot()

def OU_processes(years, timestep, num_sims, startRate, kappa, theta, sigma):
"""
timestep has to be defined in years or a fraction of years
e.g. 0.1 => 1/10th of a year; 2 => 2 years
"""
# I first generate a list of the times
times = np.arange(0,years+timestep,timestep)
#I then generate an object of standard random normal variables with a shape that depends on the number of timepoints and the number of simulations
epsilon = np.random.normal(0, 1, (num_sims, len(times)-1))
# Next, I compute the standard deviation of the stochastic component
elt = 0.5 / kappa * (1.0 - np.exp (-2.0 * kappa * timestep))
V = elt * sigma ** 2
sqrt_V = np.sqrt(V)
# I create a storage space for the values of the simulated Ornstein-Uhlenbeck processes
ou = np.zeros((num_sims, len(times)))
# I assign the starting rate to the current rate
ou[:,0] = startRate
# I generate the OU-processes by first assigning the stochastic component and then adding the deterministic component
ou[:, 1:] = np.kron(sqrt_V, np.ones((num_sims, 1))) * epsilon
for i in range(1, ou.shape[1]):
ou[:, i] += theta * (1 - np.exp(-kappa * timestep))
ou[:, i] += np.exp (-kappa * timestep) * ou[:, i-1]
ou = pd.DataFrame(np.transpose(ou))
ou.index = times
return ou



A call option on zero coupon bonds $$P(t, T)$$ imply that its price, at time $$t$$, is given by

$$V(t) = E_t^Q \left[ D(t, T) \cdot \max \left(P(\tau, T) - K, 0 \right) \right],$$

with $$t \leq \tau \leq T$$. $$P(t, T)$$ is the discount bond or zero coupon bond and it is given by:

$$P(t, T) = E_t^Q \left[ \exp \left( - \int_t^T r(s) ds \right) \right].$$

Fortunately, for this particular model we have an analytical expression. Otherwise, you would have to solve a Riccati System of Ordinary Differential Equations. The analytical solution can be found in many textbooks. If you try to solve this value by a Monte Carlo simulation then, when you compute the option price $$V(t)$$, you will get a Monte Carlo on top of another Monte Carlo, which might be computationally prohibitive.

On the other hand, $$D(t, T)$$ is the discount factor, which is given by the following expression:

\begin{aligned} \beta(t) &= \int_t^T r(s) ds, \\ D(t, T) &= \frac{\beta(t)}{\beta(T)}. \end{aligned}

with $$\beta(t)$$ the money market account or bank account.

• Thanks for your feedback. Could you provide more details on why the simulated short rate paths cannot serve both to calculate the bond price and discount factor at the same time? Why would this require a second Monte Carlo? And how would you go about implementing it? Could you maybe provide a reference? Thanks again! – Bart S Oct 12 '20 at 6:13
• Hi @BartS! Well, one thing is to realize that \begin{aligned} V(t) &= E_t^Q \left[ D(t, T) \cdot \max \left(P(\tau, T) - K, 0 \right) \right] \\ V(t) &= E_t^Q \left[ D(t, T) \cdot \max \left(E_\tau^Q \left[ \exp \left( - \int_\tau^T r(s) ds \right) \right] - K, 0 \right) \right] \end{aligned} Then, you have to compute an expectation that has a filtration until $\tau$, with $t < \tau$. You simulate until $\tau$ (the outer expectation) and, for each of those paths, start a new simulation at that point. Hope it helps! – rvignolo Oct 12 '20 at 14:35