#!/usr/bin/env python
# coding: utf-8

# # Hull/White One-Factor-Model

# In[3]:


get_ipython().run_line_magic('pylab', 'inline')
import pandas as pd
import scipy.stats


# ## Read yield curve

# In[23]:


y = pd.read_excel("M4L2 - yield curve.xlsx", "Tabelle1", index_col=0)
y


# 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)