Vasicek_and_CIR_calibration.ipynb

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

def zero_coupon(tau, r0, kappa, theta, sigma, model): #tau=(T-t),kappa=a
    #tau = tau(:); 잔존만기
    if model == 'Vasicek':
        B = (1 - np.exp(-kappa * tau)) / kappa
        A = (theta - sigma**2 / (2 * kappa**2)) * (B - tau) - (sigma**2 / (4 * kappa)) * B**2
    elif model == 'CIR':
        g = np.sqrt(kappa**2 + 2 * sigma**2)  
        tmp = 2 * kappa * theta / sigma**2
        tmp1 = kappa * tau / 2
        tmp2 = g * tau / 2
        
        A = tmp * np.log(np.exp(tmp1) / (np.cosh(tmp2) + (kappa / g) * np.sinh(tmp2)))
        #B = 2. / (kappa + g * (1. / np.tanh(g * tau / 2)))
        tanh = np.tanh(g * tau / 2)
        B = 2. * tanh / (kappa * tanh + g)
    else:
        print('zero_coupon: model must be "Vasicek" or "CIR"!')
        return -1

    #n1 = len(r0)
    n2 = len(A)
    #A_ = repmat(A,1,n1)
    A_ = A
    #B_ = repmat(B,1,n1)
    B_ = B

    #r_ = repmat(r0,n2,1)
    r_ = np.repeat(r0, n2) #r0를 n2번 반복시킨 배열을 생성

    p = np.exp(A_ - B_ * r_)
    return p

def swapRates(tau, p, mat):
    
    tmax = mat[-1]
    
    ttemp = np.arange(0.5, tmax + 0.5, 0.5)
    ptemp = np.interp(ttemp, tau, p) 
    #선형보간 x 에서 평가 된 주어진 이산 데이터 포인트 ( xp , fp )가 
    #있는 함수에 대한 1 차원 조각 별 선형 보간을 반환합니다 .

    dis = np.cumsum(ptemp)
    #np.cumsum: 각 원소들의 누적 합을 나타냄
    #dis = dis(:);

    # linear interpolation
    pmat = np.interp(mat, tau, p)

    index = (2 * mat).astype(int) - 1
    S = 100 * 2 * (1 - pmat) / dis[index]
    
    return S

def liborRates(tau, p, mat):
    pmat = np.interp(mat, tau, p)
    L = 100 * (1. / pmat - 1) / mat
    return L

def objFunc1(params, tau, LIBOR, SWAP, model):
    # unpack params
    r0 = params[0]
    kappa = params[1]
    theta = params[2]
    sigma = params[3]

    p = zero_coupon(tau, r0, kappa, theta, sigma, model)
    # now that we have zero-coupon bond prices p(t,T)
    # now it is time to calculate MODEL LIBOR rates and SWAP rates
    S = swapRates(tau, p, SWAP[:,0])
    L = liborRates(tau, p, LIBOR[:,0])

    # the goal is to minimize the distance between model rates and market rates
    rel1 = (S - SWAP[:,1]) / SWAP[:,1]
    rel2 = (L - LIBOR[:,1]) / LIBOR[:,1]

    #rel1 = (S-SWAP(:,2))
    #rel2 = (L-LIBOR(:,2))

    #mae = (sum(abs(rel1))+sum(abs(rel2)))
    mae = np.sum(rel1**2) + np.sum(rel2**2)
    
    return mae

def calibration(fun, param_0, tau, LIBOR, SWAP, model):
     # change tolerance
    opt = {'maxiter':1000, 'maxfev':5e3}
    sol = minimize(objFunc1, params1, args = (tau, LIBOR, SWAP, model), method='Nelder-Mead', options=opt)
    print(sol.message)
    par = np.array(sol.x)
    print('parameters = ' + str(par))
    r_star = par[0]
    kappa_star = par[1]   #알파
    theta_star = par[2]    #베타
    sigma_star = par[3]   #시그마는 왜 나오는거지 횡단면 자료...? curve fitting한것 같기도
    p = zero_coupon(tau, r_star, kappa_star, theta_star, sigma_star, model)
    L = liborRates(tau, p, LIBOR[:, 0])
    S = swapRates(tau, p, SWAP[:, 0])
    return p, L, S

Data

LIBOR = np.array([[1/12, 1.49078, 2.2795],
    [2/12, 1.52997, 2.33075],
    [3/12, 1.60042, 2.43631],
    [6/12, 1.76769, 2.63525],
    [12/12, 2.04263, 2.95425]])

SWAP = np.array([[2, 2.013, 3.0408],
    [3, 2.1025, 3.1054],
    [5, 2.195, 3.1332],
    [7, 2.2585, 3.1562],
    [10, 2.3457, 3.199],
    [15, 2.4447, 3.2437],
    [30, 2.5055, 3.227]])

Parameters

# we are going out 30 years every month (1/12)
tau = np.arange(0, 30 + 1/12, 1/12) #np.arange(시작점, 끝점, step size)

# initial starting point
# r0 = 0.02 kappa0 = 2 theta0 = 0.5 sigma0 = 0.1
params = np.array([0.02, 1.0, 0.05, 0.3]);

params1 = np.array([0.25, 5, 0.2, 0.1]);
params2 = np.array([0.25, 5, 0.2, 0.1]);

Vasicek (2 dates)

model = 'Vasicek'
# first date
whichOne = 1
p11, L11, S11 = calibration(objFunc1, params1, tau, LIBOR[:,[0,whichOne]], SWAP[:,[0,whichOne]], model)

Optimization terminated successfully. parameters = [0.01445667 0.9474955 0.02918339 0.0952191 ]

# second date
whichOne = 2
p12, L12, S12 = calibration(objFunc1, params1, tau, LIBOR[:,[0,whichOne]], SWAP[:,[0,whichOne]], model)

Optimization terminated successfully. parameters = [2.13256819e-02 2.82039824e+00 3.21754393e-02 1.27705822e-05]

CIR