Hull/White One-Factor-Model

In [3]:
%pylab inline
import pandas as pd
import scipy.stats
Populating the interactive namespace from numpy and matplotlib

Read yield curve

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