Monte Carlo Pricing of Options on Coupon-Bearing Bonds

# Which model for pricing an option on a coupon bearing bond?


# Which model for pricing an option on a coupon bearing bond?












I'm working on pricing a call option on a coupon bearing bond with a long duration (think something like US Treasury notes) with a payoff given by $$\max(\frac{B(T)}{B(0)}-100\%, 0)$$ rather than a usual call $\max(B(T)-B(0), 0)$ with $B(0)$ and $B(T)$ being bond prices at times $t=0$ and $t=T$ (option maturity) respectively. The market is in the beginning of rates cutting cycle so the expectation is that such options should be worth a while.

Is there a proper way to numerically price an option without applying a Jamshidian decomposition technique? Ideally I would like to just throw a bond in a Monte Carlo pricer to get a range of simulated bond prices at $t=T$ and calculate expected payoff, i.e.

- simulate short rate paths using some short rate model (Vasicek, Hull-White)

- integrate simulated short rates to discounting factors

- cut off all coupon payments happening before option maturity $T$ and discount remaining bond cashflows to maturity date $T$ using discounting factors simulated from a short rate model

- calculate payoff as an average across all simulated forward bond prices

Is it a valid approach or am I missing something important? I implemented an aforementioned approach using a Hull-White 1-Factor model calibrated on historical rates using a maximum likelihood estimation of parameters however model prices seems off market values. The underlying bond barely changes its price relative to spot after a few years making an option almost worthless which contradicts intuition of having an upside due to expected rates cut. I'm trying to justify it by Hull-White having a normal distribution of short rates therefore an average of model generated forward curves isn't much different from current market forward curve and thus the forward bond price is more or less the same.

I'm wondering is it a matter of calibration or are there any fundamental flows? Which model is a standard way to go for pricing (not necessary vanilla) options on coupon bearing bonds?

UPDATE

I decided to include some code snippets in order to better illustrate what I'm talking about. Let us consider bond valuation and skip option part for the sake of simplicity and code clarity. I will also avoid calibration routine and just state model inputs as is.

Set up of Hull-White short rate paths simulation based on QuantLib

```
import QuantLib as ql
import datetime
import numpy as np
import pandas as pd

'''
    QuantLib Hull-White interest rates simulation
'''

a = 0.84327792
sigma = 0.0213146
time_horizon = 18.0 # in years
timesteps = int(time_horizon * 365.)
dcc = ql.Actual365Fixed()
valDate = ql.Date(11, 9, 2025)

yields_g_curve = np.array([15.10, 14.39, 13.97, 13.74, 13.53, 13.60, 13.73, 13.80, 13.81, 13.76, 13.72, 13.71]) / 100. # 11 Sept 2025
tenors_g_curve = np.array([0.25, 0.50, 0.75, 1.00, 2.00, 3.00, 5.00, 7.00, 10.00, 15.00, 20.00, 30.00])
dfs = np.array([1.00] + [np.exp( - yields_g_curve[i] * tenors_g_curve[i]) for i in range(len(yields_g_curve))])
curveDate = ql.Date(11, 9, 2025)

dates = [curveDate, curveDate + ql.Period('3M'), curveDate + ql.Period('6M'), curveDate + ql.Period('9M'), curveDate + ql.Period('1Y'), 
         curveDate + ql.Period('2Y'), curveDate + ql.Period('3Y'), curveDate + ql.Period('5Y'), curveDate + ql.Period('7Y'), 
         curveDate + ql.Period('10Y'), curveDate + ql.Period('15Y'), curveDate + ql.Period('20Y'), curveDate + ql.Period('30Y')]

ql.Settings.instance().evaluationDate = valDate

curr_curve = ql.DiscountCurve(dates, dfs, dcc, ql.NullCalendar())
curr_curve_handle = ql.YieldTermStructureHandle(curr_curve)

hw_process = ql.HullWhiteProcess(curr_curve_handle, a, sigma)
rng = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(timesteps, ql.UniformRandomGenerator()))
seq = ql.GaussianPathGenerator(hw_process, time_horizon, timesteps, rng, False)

def generate_paths(num_paths, timesteps):
    arr = np.zeros((num_paths, timesteps + 1))
    for i in range(num_paths):
        sample_path = seq.next()
        path = sample_path.value()
        time = [path.time(j) for j in range(len(path))]
        value = [path[j] for j in range(len(path))]
        arr[i, :] = np.array(value)
    return np.array(time), arr

def truncateShortRatePaths(shortRatePaths, startDate, truncDate):
    idx = (truncDate - startDate).days
    output = [path[idx:] for path in shortRatePaths]
    return output

num_paths = 1000
time, paths = generate_paths(num_paths, timesteps)

truncated_paths = truncateShortRatePaths(paths, datetime.date(2025, 9, 11), datetime.date(2027, 9, 11))
```

below is the hands on implementation of Monte Carlo bond valuation based on simulated above short rate paths

```
'''
    Monte Carlo bond pricing based on simulated short rate paths
'''

def bondMonteCarlo(shortRatePaths, valDate, faceValue, maturityDate, couponRate, prevCouponDate, couponDates, n_of_sims=1000):
    allDates = pd.date_range(valDate, maturityDate, freq='d').tolist()
    allDates = [dd.date() for dd in allDates]
    cashflows = []
    for i in range(len(couponDates)):
        if i == 0:
            cf = faceValue * couponRate * (couponDates[i] - prevCouponDate).days / 365.
            cashflows.append([couponDates[i], cf])
        elif couponDates[i] == maturityDate:
            cf = faceValue + faceValue * couponRate * (couponDates[i] - couponDates[i-1]).days / 365.
            cashflows.append([couponDates[i], cf])
        else:
            cf = faceValue * couponRate * (couponDates[i] - couponDates[i-1]).days / 365.
            cashflows.append([couponDates[i], cf])

    bond_prices = []
    for i in range(n_of_sims):
        r_t = shortRatePaths[i]
        price = 0
        for j in range(len(couponDates)):
            dt = (couponDates[j] - valDate).days / 365.
            idx = allDates.index(couponDates[j])
            dt = 1 / 365
            integral_rate = np.trapezoid(r_t[:idx], dx=dt)
            discount_factor = np.exp(-integral_rate)
            price += cashflows[j][1] * discount_factor
        bond_prices.append(float(price))

    return np.mean(bond_prices)

coupons = [datetime.date(2025, 12, 4), datetime.date(2026, 6, 4), datetime.date(2026, 12, 4), datetime.date(2027, 6, 4), datetime.date(2027, 12, 4), 
           datetime.date(2028, 6, 4), datetime.date(2028, 12, 4), datetime.date(2029, 6, 4), datetime.date(2029, 12, 4), datetime.date(2030, 6, 4), 
           datetime.date(2030, 12, 4), datetime.date(2031, 6, 4), datetime.date(2031, 12, 4), datetime.date(2032, 6, 4), datetime.date(2032, 12, 4), 
           datetime.date(2033, 6, 4), datetime.date(2033, 12, 4), datetime.date(2034, 6, 4), datetime.date(2034, 12, 4), datetime.date(2035, 6, 4), 
           datetime.date(2035, 12, 4), datetime.date(2036, 6, 4), datetime.date(2036, 12, 4), datetime.date(2037, 6, 4), datetime.date(2037, 12, 4), 
           datetime.date(2038, 6, 4), datetime.date(2038, 12, 4), datetime.date(2039, 6, 4), datetime.date(2039, 12, 4), datetime.date(2040, 5, 16)]
maturity_bond = datetime.date(2040, 5, 16)

bond_price = bondMonteCarlo(shortRatePaths=paths, valDate=datetime.date(2025, 9, 11), faceValue=1000, maturityDate=maturity_bond, 
                              couponRate=0.1225, prevCouponDate=datetime.date(2025, 6, 4), couponDates=coupons, n_of_sims=num_paths)
print("Monte Carlo spot bond price: {0}".format(bond_price))
print()

truncated_coupons = [datetime.date(2027, 12, 4), datetime.date(2028, 6, 4), datetime.date(2028, 12, 4), datetime.date(2029, 6, 4), 
                     datetime.date(2029, 12, 4), datetime.date(2030, 6, 4), datetime.date(2030, 12, 4), datetime.date(2031, 6, 4), 
                     datetime.date(2031, 12, 4), datetime.date(2032, 6, 4), datetime.date(2032, 12, 4), datetime.date(2033, 6, 4), 
                     datetime.date(2033, 12, 4), datetime.date(2034, 6, 4), datetime.date(2034, 12, 4), datetime.date(2035, 6, 4), 
                     datetime.date(2035, 12, 4), datetime.date(2036, 6, 4), datetime.date(2036, 12, 4), datetime.date(2037, 6, 4), 
                     datetime.date(2037, 12, 4), datetime.date(2038, 6, 4), datetime.date(2038, 12, 4), datetime.date(2039, 6, 4), 
                     datetime.date(2039, 12, 4), datetime.date(2040, 5, 16)]
fwd_bond_price = bondMonteCarlo(shortRatePaths=truncated_paths, valDate=datetime.date(2027, 9, 11), faceValue=1000, maturityDate=maturity_bond, 
                              couponRate=0.1225, prevCouponDate=datetime.date(2027, 6, 4), couponDates=truncated_coupons, n_of_sims=num_paths)
print("Monte Carlo forward bond price after 2 years: {0}".format(fwd_bond_price))
```

We are looking at a bond with semi-annual $12,25\%$ coupons maturing on 16 May 2040. One can run the code and check that the spot price is $91,28\%$ and the forward price as of September 2027 is $91,44\%$. There is barely any upside despite an inverted yield curve and market expectations of rates cut.

Let us double check this implementation with QuantLib spot and forward pricing based just on the current forward curve

```
'''
    QuantLib coupon bearing bond pricing
'''

valuationDate = ql.Date(11, 9, 2025)
ql.Settings.instance().evaluationDate = valuationDate
compounding = ql.Compounded
calendar = ql.NullCalendar()
coupon = 0.1225
couponFrequency = ql.Semiannual
issueDate = ql.Date(15, 5, 2024)
maturityDate = ql.Date(16, 5, 2040)
settlementDays = 2
settlementDate = calendar.advance(issueDate, ql.Period(settlementDays, ql.Days))
dayCount = ql.Actual365Fixed()
schedule = ql.Schedule(ql.Date(4, 12, 2024), maturityDate, ql.Period(couponFrequency), calendar, ql.Unadjusted, ql.Unadjusted, ql.DateGeneration.Forward, True)
fixedRateBond = ql.FixedRateBond(settlementDays, 1000, schedule, [coupon], ql.Actual365Fixed())
curve = curr_curve
handle = ql.YieldTermStructureHandle(curve)
bondEngine = ql.DiscountingBondEngine(handle)
ofz = fixedRateBond.setPricingEngine(bondEngine)
print('QuantLib spot bond price:', fixedRateBond.NPV())
print()

ql.Settings.instance().evaluationDate = ql.Date(11, 9, 2027)
fwd_yts = ql.ImpliedTermStructure(handle, ql.Date(11, 9, 2027))
fwd_handle = ql.YieldTermStructureHandle(fwd_yts)
fwdBondEngine = ql.DiscountingBondEngine(fwd_handle)
fwd_ofz = fixedRateBond.setPricingEngine(fwdBondEngine)
print('QuantLib forward bond price:', fixedRateBond.NPV())
```

the results are $91,17\%$ for spot and $91,36\%$ for forward prices which are pretty close to Hull-White Monte Carlo.