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

Question

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()

```

Answer 1

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