# Hull/White One-Factor-Model¶

In [3]:
%pylab inline
import pandas as pd
import scipy.stats

Populating the interactive namespace from numpy and matplotlib


In [23]:
y = pd.read_excel("M4L2 - yield curve.xlsx", "Tabelle1", index_col=0)
y

Out[23]:
yield
t
0.25 0.0020
0.50 0.0050
0.75 0.0070
1.00 0.0110
1.50 0.0150
2.00 0.0180
3.00 0.0200
4.00 0.0220
5.00 0.0250
10.00 0.0288
20.00 0.0310
30.00 0.0340
In [5]:
Z = []
logZ = []
for t in y.index:
Z.append(exp(-y["yield"][t]*t))
logZ.append(-y["yield"][t]*t)
Z = pd.Series(Z, index=y.index)
logZ = pd.Series(logZ, index=y.index)


## Model parameters¶

In [6]:
# Model parameters

a = 10     # Speed of reversion
c = 0.1    # Volatility
r = 0.03   # spot interest rates

gamma = 1/a


## Yield curve¶

In [7]:
fig, ax = plt.subplots()
ax.plot(y.index, y["yield"], 'o-')
ax.set(xlabel='time (y)', ylabel='yield (%)',
title='Yield curve')
ax.grid()
plt.show()

In [8]:
def make_spline(x, a):
n = len(x) - 1
h = []
alpha = zeros(n)
for i in range(n):
h.append(x[i+1] - x[i])

for i in range(1, n):
alpha[i] = 3 / h[i] * (a[i+1] - a[i]) - 3 / h[i-1] * (a[i] - a[i-1])

l = zeros(n+1)
l[0] = 1
mu = zeros(n+1)
z = zeros(n+1)
for i in range(1, n):
l[i] = 2 * (x[i+1] - x[i-1]) - h[i-1] * mu[i-1]
mu[i] = h[i] / l[i]
z[i] = (alpha[i] - h[i-1] * z[i-1]) / l[i]

l[n] = 1
z[n] = 0
b = zeros(n)
c = zeros(n+1)
d = zeros(n)

for j in range(n-1, -1, -1):
c[j] = z[j] - mu[j] * c[j+1]
b[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3
d[j] = (c[j+1] - c[j])/(3*h[j])

def spline(t):
if t < x[0]:
t = x[0]
if t >= x[-1]:
t = x[-1]
for j in range(n):
if t >= x[j] and t <= x[j+1]:
s = a[j] + b[j]*(t-x[j]) + c[j]*(t - x[j])**2 + d[j]*(t - x[j])**3
s1 = b[j] + 2 * c[j] * (t - x[j]) + 3 * d[j] * (t - x[j])**2
s2 = 2 * c[j] + 6 * d[j] * (t - x[j])
return (s, s1, s2)
return None
return spline


## Interpolate yield curve with cubic spline¶

In [9]:
yield_spline = make_spline(array(y.index), array(y["yield"]))

T = arange(y.index[0], y.index[-1], 0.1)
S = []
for t in T:
(s, s1, s2) = yield_spline(t)
S.append(s)

fig, ax = plt.subplots()
ax.plot(y.index, y["yield"], 'o', T, S)
ax.set(xlabel='time (y)', ylabel='yield (%)',
title='Yield curve (Spline interpolated)')
ax.grid()
plt.show()


## Zero coupon prices¶

In [10]:
fig, ax = plt.subplots()
ax.plot(y.index, Z, 'o-')

ax.set(xlabel='time (y)', ylabel='Z',
title='zero coupon prices')
ax.grid()
plt.show()


## calculate log(Z)¶

In [11]:
fig, ax = plt.subplots()
ax.plot(y.index, logZ, 'o-')

ax.set(xlabel='time (y)', ylabel='log(Z)',
title='log Z')
ax.grid()
plt.show()


## Interpolate log(Z) with cubic splines¶

In [22]:
Z_spline = make_spline(array(y.index), array(Z))
logZ_spline = make_spline(array(y.index), array(logZ))

T = arange(y.index[0], y.index[-1], 0.1)
S = []
S1 = []
S2 = []
for t in T:
(s, s1, s2) = logZ_spline(t)
S.append(s)
S1.append(s1)
S2.append(s2)

fig, ax = plt.subplots()
ax.plot(y.index, logZ, 'o', T, S)
ax.set(xlabel='time (y)', ylabel='log Z',
title='log Z (Spline interpolated)')
ax.grid()
plt.show()


## Instantaneous forward rates¶

In [13]:
T = arange(y.index[0], y.index[-1], 0.01)
S1 = []
for t in T:
(s, s1, s2) = logZ_spline(t)
S1.append(-s1)

fig, ax = plt.subplots()
ax.plot(T, S1, y.index, y["yield"])
ax.set(xlabel='time (y)', ylabel='f',
title='Inst. forward rates')
ax.grid()
plt.show()


## Set model parameters¶

Hull-White SDE for evolution of interest rates: $$\mathrm d r = \left(\eta(t) - \gamma r\right)\mathrm d t + c\, \mathrm d X$$

## Calculate $\eta(t)$¶

$$\eta^*(t) = - \frac{\partial^2}{\partial t^2}\operatorname{log}(Z_M(t^*; t)) - \gamma\frac{\partial}{\partial t}\operatorname{log}(Z_M(t^*; t)) + \frac{c^2}{2\gamma}\left(1-\operatorname{e}^{-2\gamma(t-t^*)}\right)$$
In [14]:
def eta(t):
(s, s1, s2) = logZ_spline(t)
return -s2 - gamma*s1 + c**2/(2*gamma)*(1-np.exp(-2*gamma*t))

T = arange(y.index[0], y.index[-1], 0.1)
E = []
for t in T:
E.append(eta(t))

fig, ax = plt.subplots()
ax.plot(T, E)
ax.set(xlabel='time (y)', ylabel='eta(t)',
title='Eta')
ax.grid()
plt.show()


## Plot long term mean level¶

In [15]:
T = arange(y.index[0], y.index[-1], 0.1)
r_mean = []
for t in T:
r_mean.append(eta(t)*a)

fig, ax = plt.subplots()
ax.plot(T, r_mean)
ax.set(xlabel='time (y)', ylabel='eta(t)',
title='r_mean')
ax.grid()
plt.show()


## Evolution of interest rates¶

In [16]:
dt = 0.1

times = arange(y.index[0], y.index[-1], dt)
rates = []

rates = []
r_ = r
for t in times:
dr = (eta(t) - gamma * r_) * dt + c * np.sqrt(dt) * np.random.normal()
r_ = r_ + dr
rates.append(r_)

fig, ax = plt.subplots()
ax.plot(times, rates)
ax.set(xlabel='time (y)', ylabel='r',
title='interest rates')
ax.grid()
plt.show()


## Zero curve according to Hull White model¶

In [17]:
def B(t, T):
return 1/gamma*(1-exp(-gamma*(T-t)))

def A(t, T):
Z_T = Z_spline(T)[0]
Z_t = Z_spline(t)[0]
logZ_t1 = logZ_spline(t)[1]
return log(Z_T/Z_t)-B(t, T)*logZ_t1 - c**2/(4*gamma**3)*(np.exp(-gamma*T) - np.exp(-gamma*t))**2*(np.exp(2*gamma*t)-1)

def HullWhite_Z(t, T):
return exp(A(t, T) - r*B(t, T))

times = arange(y.index[0], y.index[-1], 0.1)
HWZ = []

for t in times:
HWZ.append(HullWhite_Z(0, t))

fig, ax = plt.subplots()
ax.plot(y.index, Z, 'o', T, HWZ)
ax.set(xlabel='time (y)', ylabel='Z',
title='Hull White Z')
ax.grid()
plt.show()


## Implementation according to Brigo¶

In [18]:
def eta(t):
(s, s1, s2) = logZ_spline(t)
return s2 + gamma*s1 + c**2/(2*gamma)*(1-np.exp(-2*gamma*t))

T = arange(y.index[0], y.index[-1], 0.1)
E = []
for t in T:
E.append(eta(t))

fig, ax = plt.subplots()
ax.plot(T, E)
ax.set(xlabel='time (y)', ylabel='eta(t)',
title='Eta')
ax.grid()
plt.show()


## Bond pricing¶

In [19]:
a = gamma  # Gamma is named a in Brigo

dt = 0.1

def B_(t, T):
global a
return 1/a*(1-exp(-a*(T-t)))

def A_(t, T):
P_T = Z_spline(T)[0]
P_t = Z_spline(t)[0]
f_t = -logZ_spline(t)[1]
return P_T/P_t*exp(B_(t, T)*f_t - c**2/(4*a)*(1-exp(-2*a*t))*B_(t, T)**2)

def Brigo_Z(t, T):
return A_(t, T)*exp(-r*B_(t, T))

times = arange(y.index[0], y.index[-1], dt)
BZ = []

for t in times:
BZ.append(Brigo_Z(0, t))

fig, ax = plt.subplots()
ax.plot(y.index, Z, 'o', T, BZ)
ax.set(xlabel='time (y)', ylabel='Z',
title='Brigo Z')
ax.grid()
plt.show()


## Yield curve¶

In [20]:
dt = 0.1
times = arange(y.index[0], y.index[-1], dt)
BY = []

for t in times:
z = Brigo_Z(0, t)
by = -log(z)/t
BY.append(by)

fig, ax = plt.subplots()
ax.plot(y.index, y["yield"], 'o', times, BY)
ax.set(xlabel='time (y)', ylabel='Z',
title='Hull White (Brigo) yield curve')
ax.grid()
plt.show()


## Option pricing¶

In [21]:
# European Bond Option for Zero coupon bond
def ZBC(t, T, S, X):
sigma_p = c * sqrt((1-exp(-2*a*(T-t)))/(2*a)) * B_(T, S)
h = 1/sigma_p * log(Brigo_Z(t, S)/(Brigo_Z(t, T) * X)) + sigma_p / 2
return Brigo_Z(t, S) * scipy.stats.norm.cdf(h) - X * Brigo_Z(t, T) * scipy.stats.norm.cdf(h - sigma_p)

def ZBP(t, T, S, X):
sigma_p = c * sqrt((1-exp(-2*a*(T-t)))/(2*a)) * B_(T, S)
h = 1/sigma_p * log(Brigo_Z(t, S)/(Brigo_Z(t, T) * X)) + sigma_p / 2
return X * Brigo_Z(t, T) * scipy.stats.norm.cdf(-h + sigma_p) - Brigo_Z(t, S) * scipy.stats.norm.cdf(-h)

def Cap(t, T, N, X):
# t = settlement date
# T = array of reset dates (> t)
# N = Nominal
# X = Strike
s = 0
for i in range(1, len(T)):
tau = T[i-1] - T[i]
s += (1 + X*tau) * ZBP(t, T[i-1], T[i], 1/(1+X*tau))
return N * s

ZBC(0, 1, 2, 0.4)

Cap(0, [0.25, 0.5, 0.75, 1], 100, 0.25)

Out[21]:
20.899908991971405