I’m simulating interest rates via the HullWhite One factor model. To simulate the short rate I’m using the code from the Quantlib Python Cookbook, chapter 15 and beyond (By Goutham Balaraman and Luigi Ballabio). The code can be found here:

The code generates the short rate. Unfortunately the code ends with having only the simulated short rate while my goal is to create a whole new term structure based on the simulated curve.

To create a term structure I need to calibrate theta based on the initial term structure. Besides that I also need to calculate A(t) and B(t) such that I’m able to create the term structure.

Questions: Does somebody have an algorithm to calculate theta, A(t) and B(t)?


You can check out here a blog post on simulating the yield term structure for the HullWhite model.

The basic idea is that once you have the paths for the short rate, you can simply integrate (approximately) the short rate throughout each path to obtain the discount factors.

The average of the simulations should match the initial term structure.

Here is an example based on that excellent blog post, with an improvement on execution time by vectorizing the integration.

import QuantLib as ql
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import copy 

nPaths = 500
years = 30

startDate = ql.Date(3, 12, 2018)
endDate = startDate + ql.Period(years, ql.Years)
tenor = ql.Period(1, ql.Days)
schedule = ql.MakeSchedule(startDate, endDate, tenor)

dates = [dt for dt in schedule]

times = [ql.Actual360().yearFraction(startDate, dt) for dt in dates]

curve = ql.YieldTermStructureHandle(ql.FlatForward(startDate, 0.04875825, ql.Actual365Fixed()))
reversionSpeed = 0.01
rateVolatility = 0.001
process = ql.HullWhiteProcess(curve, reversionSpeed, rateVolatility)

periods = (endDate - startDate) + 1
sequenceGenerator = ql.UniformRandomSequenceGenerator(periods, ql.UniformRandomGenerator())
gaussianSequenceGenerator = ql.GaussianRandomSequenceGenerator(sequenceGenerator)
pathGenerator = ql.GaussianPathGenerator(process, years, periods, gaussianSequenceGenerator, False)

paths = np.zeros(shape = (nPaths, periods))
for i in range(nPaths):
    path = pathGenerator.next().value()
    paths[i, :] = np.array([path[j] for j in range(periods)])
dfs = np.zeros(shape=paths.shape)
dt = years / len(schedule)
integral = 0
for j in range(periods):
    if j == 0:
        dfs[:, 0] = 1
        integral = copy.deepcopy(paths[:, 0])
        integral += paths[:, j]        
        dfs[:, j] += np.exp(-integral * dt)
simulatedCurve = ql.DiscountCurve(dates, dfs.mean(0), ql.Actual365Fixed(), ql.NullCalendar())
simulatedDfs = np.array([simulatedCurve.discount(dt) for dt in dates])
dfs_curve = np.array([curve.discount(dt) for dt in dates])

fig, ax = plt.subplots(1,2, figsize=(10,3))

ax[0].set_title('Discount Factors')
ax[0].plot(times, simulatedDfs, linestyle = 'dashed', label = 'simulated curve')
ax[0].plot(times, dfs_curve, linestyle = 'solid', label = 'original curve')

ax[1].plot(paths.T, linewidth=0.5)

The output would be:

enter image description here

  • $\begingroup$ Can you explain what the 0.04875825 represents and how it is determined? $\endgroup$
    – Andrew
    May 29 '20 at 13:56
  • $\begingroup$ It's a flat forward rate. It's not determined and is just a fixed input in this case. Alternatively, you could bootstrap a whole curve from market instruments $\endgroup$ May 29 '20 at 14:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.