# Quantlib: How do I price a ZC bond using the Hull White model?

I am trying to use QuantLib to model short rate and looks like QL has some material here http://gouthamanbalaraman.com/blog/hull-white-simulation-quantlib-python.html

I have been able to simulate using the actual term structure, although in the attached link for simplicity GB has assumed flat forwards. Which is what I have used below.

I am calculating the price of the zero coupon bond at a time t using the below code

def price(T, t, r0):
tau = T - t
B = (1 - np.exp(-a * tau)) / a
A = np.exp(-r0 * tau - B * r0 - sigma ** 2 / (4 * a ** 3) *
(np.exp(-a * T) - np.exp(-a * t)) * (np.exp(2 * a * t) - 1))
return A * np.exp(-r0 * B)


but I feel my zc bond prices are all over the place. Any clues what an I doing wrong here?

Does QuantLib provide a wrapper to calculate the zc prices using the HW model by any chance?

The complete code is below

import QuantLib as ql
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

sigma = 0.1
a = 0.1
length = 30 # in years
steps_per_year = 12
timestep = length * steps_per_year
forward_rate = 0.05
day_count = ql.Thirty360()
todays_date = ql.Date(15, 1, 2015)

ql.Settings.instance().evaluationDate = todays_date

spot_curve = ql.FlatForward(todays_date, ql.QuoteHandle(ql.SimpleQuote(forward_rate)), day_count)
spot_curve_handle = ql.YieldTermStructureHandle(spot_curve)

hw_process = ql.HullWhiteProcess(spot_curve_handle, a, sigma)
rng = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(timestep, ql.UniformRandomGenerator()))
seq = ql.GaussianPathGenerator(hw_process, length, timestep, rng, False)

def generate_paths(n_scenarios):
arr = np.zeros((n_scenarios, timestep+1))
for i in range(n_scenarios):
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

n_scenarios = 1
time, paths = generate_paths(n_scenarios)
rates = pd.DataFrame(paths).T

#price a zero coupon bond

n_years = length
num_steps = timestep
dt = 1 / steps_per_year
prices = np.empty_like(rates)

def price(T, t, r0):
tau = T - t
B = (1 - np.exp(-a * tau)) / a
A = np.exp(-r0 * tau - B * r0 - sigma ** 2 / (4 * a ** 3) *
(np.exp(-a * T) - np.exp(-a * t)) * (np.exp(2 * a * t) - 1))
return A * np.exp(-r0 * B)

prices[0] = price(n_years,dt*0,rates.values[0])

for steps in range(1,num_steps+1):
prices[steps] = price(n_years, dt * steps, rates.values[steps])

plt.plot(prices)
plt.show()

$$$$


Here is the price in HW[4] for a ZCB at time $$t$$:

\begin{align} P(t,T) &= A(t,T) e^{-B(t,T) r(t)}\\ A(t,T) &= {\frac {P(0,T)} {P(0,t)}} \exp \Bigl( B(t,T)F(0,t) - {\frac {\sigma^2} {4a}} B(t,T)^2(1-e^{-2at})\Bigr)\\ B(t,T) &= {\frac {1-e^{-a(T-t)}} {a}} \end{align}

You seem to be simulating to rate $$r(t)$$ at time $$t$$ and putting that into your function argument r0. I'm not exactly sure what r0 is supposed to be, but note that the equations above include $$F(0,t)$$ - which is just a constant in your simulation (as you've used a FlatForward curve which I think should replace most of your r0 terms, but the final line should read return A * np.exp(-rt * B), where rt is the current value of rate. There is also a sign error on the B * r0 term, which should be + instead of -

When I make these changes (and slightly reduce the rates vol, as 10% for 30Y is a long time and a lot of vol for a short rates model), and average the prices along the paths, I see the following, which nicely shows the bond's pull-to-par as its maturity approaches:

Full modified code:

import QuantLib as ql
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

sigma = 0.01
a = 0.1
length = 30 # in years
steps_per_year = 12
timestep = length * steps_per_year
forward_rate = 0.01
day_count = ql.Thirty360()
todays_date = ql.Date(15, 1, 2015)

ql.Settings.instance().evaluationDate = todays_date

spot_curve = ql.FlatForward(todays_date, ql.QuoteHandle(ql.SimpleQuote(forward_rate)), day_count)
spot_curve_handle = ql.YieldTermStructureHandle(spot_curve)

hw_process = ql.HullWhiteProcess(spot_curve_handle, a, sigma)
rng = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(timestep, ql.UniformRandomGenerator()))
seq = ql.GaussianPathGenerator(hw_process, length, timestep, rng, False)

def generate_paths(n_scenarios):
arr = np.zeros((n_scenarios, timestep+1))
for i in range(n_scenarios):
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

n_scenarios = 1000
time, paths = generate_paths(n_scenarios)
rates = pd.DataFrame(paths).T

#price a zero coupon bond
n_years = length
num_steps = timestep
dt = 1 / steps_per_year
prices = np.empty_like(rates)

def price(T, t, f0, rt):
tau = T - t
B = (1 - np.exp(-a * tau)) / a
A = np.exp(-f0 * tau + B * f0 - sigma ** 2 / (4 * a ** 3) *
(np.exp(-a * T) - np.exp(-a * t)) * (np.exp(2 * a * t) - 1))
return A * np.exp(-rt * B)

prices[0] = price(n_years,dt*0, forward_rate, rates.values[0])

for steps in range(1,num_steps+1):
prices[steps] = price(n_years, dt * steps, forward_rate, rates.values[steps])

fig = plt.figure()
plt.figure(figsize=(16,10))

plt.subplot(2, 1, 1)
plt.title("All ZCB paths")
plt.plot(prices)

plt.subplot(2, 1, 2)

ave = [np.mean(x) for x in prices]
plt.plot(ave)
plt.title("Average ZCB paths")

plt.show()


[4]: Hull, J., White, A., “Pricing Interest-Rate-Derivative Securities”, Review of Financial Studies, Volume 3, Issue 4, pp. 573-592, 1990

• thanks @StackG, that was really helpful. – TRex Aug 11 '20 at 7:44
• thanks again, this was really helpful. The equation makes sense now, although had a question and excuse my ignorance as I haven't dabbled with this before. but instead of defining a constant forward rate I pass the actually current term structure, whats the most efficent way to get F(0,t), I believe its instantaneous fwd rate as seen at time 0 betw t and T? – TRex Aug 11 '20 at 9:30
• $F(0,t)$ is the initial rate from $0$ to $t$, you can get it from the ql.YieldTermStructureHandle object (here, you've assumed it is constant). It doesn't depend on the ZCB's maturity time $T$ – StackG Aug 11 '20 at 10:02
• thanks again man! – TRex Aug 11 '20 at 12:33
• looks like F(0,t) has been a struggle. Assuming that I get the current term structure from here https://quant.stackexchange.com/questions/55992/quantlib-how-do-i-price-a-bond-after-having-built-a-term-structure/55997?noredirect=1#comment80321_55997. So my spot_curve_handle in the above code will change to spot_curve_handle = ql.YieldTermStructureHandle(yieldcurve)`. now, how do I get new F(0,t) for each num_steps? Appreciate your help again. – TRex Aug 12 '20 at 13:22