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