#
# Calibration of CIR85 model to Euribor Rates
#
#
import sys
sys.path.append('10_mcs')
import math
import numpy as np
np.set_printoptions(suppress=True,
formatter={'all': lambda x: '%7.6f' % x})
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.family'] = 'serif'
import scipy.interpolate as sci
from scipy.optimize import fmin
from CIR_zcb_valuation_gen import B
#
# Market Data: Eonia rate (01.10.2014) + Euribor rates
# Source: http://www.emmi-benchmarks.eu
# on 30. September 2014
#
t_list = np.array((1, 7, 14, 30, 60, 90, 180, 270, 360)) / 360.
r_list = np.array((-0.032, -0.013, -0.013, 0.007, 0.043,
0.083, 0.183, 0.251, 0.338)) / 100
factors = (1 + t_list * r_list)
zero_rates = 1 / t_list * np.log(factors)
r0 = r_list[0]
#
# Interpolation of Market Data
#
tck = sci.splrep(t_list, zero_rates, k=3) # cubic splines
tn_list = np.linspace(0.0, 1.0, 24)
ts_list = sci.splev(tn_list, tck, der=0)
de_list = sci.splev(tn_list, tck, der=1)
f = ts_list + de_list * tn_list
# forward rate transformation
def plot_term_structure():
plt.figure(figsize=(8, 5))
plt.plot(t_list, r_list, 'ro', label='rates')
plt.plot(tn_list, ts_list, 'b', label='interpolation', lw=1.5)
# cubic splines
plt.plot(tn_list, de_list, 'g--', label='1st derivative', lw=1.5)
# first derivative
plt.legend(loc=0)
plt.xlabel('time horizon in years')
plt.ylabel('rate')
plt.grid()
#
# Model Forward Rates
#
def CIR_forward_rate(opt):
''' Function for forward rates in CIR85 model.
Parameters
==========
kappa_r: float
mean-reversion factor
theta_r: float
long-run mean
sigma_r: float
volatility factor
Returns
=======
forward_rate: float
forward rate
'''
kappa_r, theta_r, sigma_r = opt
t = tn_list
g = np.sqrt(kappa_r ** 2 + 2 * sigma_r ** 2)
sum1 = ((kappa_r * theta_r * (np.exp(g * t) - 1)) /
(2 * g + (kappa_r + g) * (np.exp(g * t) - 1)))
sum2 = r0 * ((4 * g ** 2 * np.exp(g * t)) /
(2 * g + (kappa_r + g) * (np.exp(g * t) - 1)) ** 2)
forward_rate = sum1 + sum2
return forward_rate
#
# Error Function
#
def CIR_error_function(opt):
''' Error function for CIR85 model calibration. '''
kappa_r, theta_r, sigma_r = opt
if 2 * kappa_r * theta_r < sigma_r ** 2:
return 100
if kappa_r < 0 or theta_r < 0 or sigma_r < 0.001:
return 100
forward_rates = CIR_forward_rate(opt)
MSE = np.sum((f - forward_rates) ** 2) / len(f)
# print opt, MSE
return MSE
#
# Calibration Procedure
#
def CIR_calibration():
opt = fmin(CIR_error_function, [1.0, 0.02, 0.1],
xtol=0.00001, ftol=0.00001,
maxiter=300, maxfun=500)
return opt
#
# Graphical Results Output
#
def plot_calibrated_frc(opt):
''' Plots market and calibrated forward rate curves. '''
forward_rates = CIR_forward_rate(opt)
plt.figure(figsize=(8, 7))
plt.subplot(211)
plt.grid()
plt.ylabel('forward rate $f(0,T)$')
plt.plot(tn_list, f, 'b', label='market')
plt.plot(tn_list, forward_rates, 'ro', label='model')
plt.legend(loc=0)
plt.axis([min(tn_list) - 0.05, max(tn_list) + 0.05,
min(f) - 0.005, max(f) * 1.1])
plt.subplot(212)
plt.grid(True)
wi = 0.02
plt.bar(tn_list - wi / 2, forward_rates - f, width=wi)
plt.xlabel('time horizon in years')
plt.ylabel('difference')
plt.axis([min(tn_list) - 0.05, max(tn_list) + 0.05,
min(forward_rates - f) * 1.1, max(forward_rates - f) * 1.1])
plt.tight_layout()
def plot_zcb_values(p0, T):
''' Plots unit zero-coupon bond values (discount factors). '''
t_list = np.linspace(0.0, T, 20)
r_list = B([r0, p0[0], p0[1], p0[2], t_list, T])
plt.figure(figsize=(8, 5))
plt.plot(t_list, r_list, 'b')
plt.plot(t_list, r_list, 'ro')
plt.xlabel('time horizon in years')
plt.ylabel('unit zero-coupon bond value')
plt.grid()
