%pylab inline
import pandas as pd
import scipy.stats
Populating the interactive namespace from numpy and matplotlib
y = pd.read_excel("M4L2 - yield curve.xlsx", "Tabelle1", index_col=0)
y
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 |
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
a = 10 # Speed of reversion
c = 0.1 # Volatility
r = 0.03 # spot interest rates
gamma = 1/a
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()
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
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()
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()
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()
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()
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()
Hull-White SDE for evolution of interest rates: $$ \mathrm d r = \left(\eta(t) - \gamma r\right)\mathrm d t + c\, \mathrm d X $$
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()
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()
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()
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()
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()
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()
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()
# 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)
20.899908991971405