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

Question

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":
        mask = pricePaths.iloc[-1,:].values>=strike
        pricePaths = pricePaths.iloc[:,mask]
        strikePaths = strikePaths.iloc[:,mask]
        sumPathsPV = np.sum(pricePaths.iloc[0,:]-strikePaths.iloc[0,:])

    elif flag.lower()=="put":
        mask = pricePaths.iloc[-1,:].values<=strike
        pricePaths = pricePaths.iloc[:,mask]
        strikePaths = strikePaths.iloc[:,mask]
        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

Answer 1

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.