In [156]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.stats as st
import seaborn as sns
sns.set(font=['IPAmincho'])
In [5]:
#Simulate interest rate oath by the Vasicek model
def vasicek(r0, K, theta, sigma, T=1., N=10, seed=777):
np.random.seed(seed)
dt = T/float(N)
rates = [r0]
for i in range(N):
dr = K*(theta-rates[-1]*dt) + sigma*np.random.normal()
rates.append(rates[-1] + dr)
return range(N+1), rates
In [6]:
#Simulate interest rate oath by the Vasicek model
def CIR(r0, K, theta, sigma, T=1., N=10, seed=777):
np.random.seed(seed)
dt = T/float(N)
rates = [r0]
for i in range(N):
dr = K*(theta-rates[-1]*dt) + sigma*math.sqrt(rates[-1])*np.random.normal()
rates.append(rates[-1] + dr)
return range(N+1), rates
In [7]:
x, y = vasicek(0.01875, 0.20, 0.0001, 0.012, 10., 200)
plt.figure(figsize=(10,5))
plt.plot(x, y)
Out[7]:
In [8]:
x, y = CIR(0.01875, 0., 0.0001, 0.012, 10., 200)
plt.figure(figsize=(10,5))
plt.plot(x, y)
Out[8]:
In [9]:
x, y = vasicek(0.01875, 0.20, 0.01, 0.012, 10., 200)
plt.figure(figsize=(10,5))
plt.plot(x, y)
Out[9]:
In [10]:
x, y = vasicek(0.01875, 0.20, 0.01, 0.012, 10., 200, seed=666)
plt.figure(figsize=(10,5))
plt.plot(x, y)
Out[10]:
In [11]:
x, y = vasicek(0.01875, 0.20, 0.01, 0.012, 10., 200, seed=888)
plt.figure(figsize=(10,5))
plt.plot(x, y)
Out[11]:
In [12]:
# Get zero coupon bond price by Vasicek model
def exact_zcb(theta, kappa, sigma, tau, r0=0.):
B = (1 - np.exp(-kappa*tau))/kappa
A = np.exp((theta - (sigma**2)/(2*(kappa**2)))*(B - tau) - (sigma**2)/(4*kappa)*(B**2))
return A*np.exp(-r0*B)
In [13]:
Ts = np.r_[0.0:25.5:0.5]
zcbs = [exact_zcb(0.5, 0.02, 0.03, t, 0.015) for t in Ts]
plt.figure(figsize=(10,5))
plt.title("Zero Coupon Bond (ZCB) Value by Time")
plt.plot(Ts, zcbs, label='ZCB')
plt.ylabel("Value ($)")
plt.xlabel("Time in years")
plt.legend()
plt.grid(True)
plt.show()
In [14]:
def exercise_value(K, R, t):
return K*math.exp(-R*t)
In [15]:
Ts = np.r_[0.0:25.5:0.5]
Ks = [exercise_value(0.95, 0.015, t) for t in Ts]
zcbs = [exact_zcb(0.5, 0.02, 0.03, t, 0.015) for t in Ts]
plt.figure(figsize=(10,5))
plt.title("Zero Coupon Bond (ZCB) Value by Time "
"and Strike(K) Values by Time")
plt.plot(Ts, zcbs, label='ZCB')
plt.plot(Ts, Ks, label='K', linestyle="--", marker=".")
plt.ylabel("Value ($)")
plt.xlabel("Time in years")
plt.legend()
plt.grid(True)
plt.show()
債券の発行者はコールを行使する権利を保有しているので、コール付きゼロクーポン債券の価格は、
$$callable~zero~coupon~bond~price~=~min(ZCB, K（割引）)$$
となる.
この債券価格は現在の金利の水準が与えられた時の近似である.
行使価格は未来に行使すればするほど価値が下がる.
上の図だと3年ぐらいで変わる.
債権の価値は3年より前でZCBだが、3年よりあとで行使価格になる.
次のステップとして、policy iterationの形で、早期行使を考慮することが考えられる.
In [16]:
x = [1, 1, 0]
all((j==1 for j in x))
Out[16]:
In [154]:
class VasicekCZCB:
def __init__(self):
self.norminv = st.distributions.norm.ppf
#使ってない.
self.nor = st.distributions.norm.cdf
#メインの関数。これを呼び出すことで、金利とそれに対応する価値を計算する.
#r0: t=0における短期金利.
#R: 「行使価格」に対する金利. 割り引くため.
#ratio: 行使価格（割合として与えられる）.
#sigma: 短期金利のボラティリティ.
#theta: 短期金利のmean.
#M: the number of steps. 時間の分割.
#prob: 上限と下限を決める確率.
#grid_struct_const: 状態の変化幅...ではない.
def vasicek_czcb_values(self, r0, R, ratio, T, sigma, kappa, theta, M, prob=1e-6, max_policy_iter=10,
grid_struct_const=0.25, rs=None):
#vasicek modelのパラメーターを生成する.
#N: 状態の分割数.
#N*Mのグリッドになる.
r_min, dr, N, dtau = self.vasicek_params(r0, M, sigma, kappa, theta, T, prob, grid_struct_const, rs)
#状態の初期設定.
#計算には使わない.
r = np.r_[0:N]*dr + r_min
#Nの数は異なってもrの最小値と最大値は一致.
#r_0は初期値. 気にするな.
#値を入れる配列をあらかじめ作っておく. これにより計算時間を減らせる.
v_mplus1 = np.ones(N)
#時間を進めていく.
#iterateでは、ある時間が与えられた時に、翌期の
for i in range(1, M+1):
#割り引かれた行使価格. i*dtauは行使までの時間.
K = self.exercise_call_price(R, ratio, i*dtau)
#今行使したら.
#状態個分の行使価格（割引）.
eex = np.ones(N)*K
#陰的差分法の各項を決める.
#ここに書く必要はない. 書かない方がいい. forループの上に出すべき?
#時間を固定してすべての状態について導出.
subdiagonal, diagonal, superdiagonal = self.vasicek_diagonals(sigma, kappa, theta, r_min, dr, N, dtau)
v_mplus1, iterations = self.iterate(subdiagonal, diagonal, superdiagonal, v_mplus1, eex, max_policy_iter)
return r, v_mplus1
#vasicek modelのパラメーターを生成する.
def vasicek_params(self, r0, M, sigma, kappa, theta, T, prob, grid_struct_const=0.25, rs=None):
#金利の上限と下限. rsがNonneなら計算する. 所与としてできるということ. デフォルトは計算.
(r_min, r_max) = (rs[0], rs[-1]) if not rs is None else self.vasicek_limits(r0, sigma, kappa, theta, T, prob)
dt = T/float(M)
#グリッッドの時間幅はわかるが、価値の変化幅は分からないので.
N = self.calculate_N(grid_struct_const, dt, sigma, r_max, r_min)
dr = (r_max-r_min)/(N-1)
#上限は返さなくていい. コールに限れば.
return r_min, dr, N, dt
def calculate_N(self, max_structure_const, dt, sigma, r_max, r_min):
N = 0
while True:
N += 1
#grid_structure_interval = sigma^2*dt/dr^2
grid_structure_interval = dt*(sigma**2)/(((r_max-r_min)/float(N))**2)
#N=0だと一個の目が非常に大きい.
#少しずつ絞る.
#grid_structure_intervalは徐々に大きくなる.
if grid_structure_interval > max_structure_const:
break
return N
def vasicek_limits(self, r0, sigma, kappa, theta, T, prob=1e-6):
er = theta+(r0-theta)*math.exp(-kappa*T)
variance = (sigma**2)*T if kappa==0 else (sigma**2)/(2*kappa)*(1-math.exp(-2*kappa*T))
stdev = math.sqrt(variance)
#meanとstdを与えて正規分布を得る。probを与えれば、そのprobabilityに対応する値が返ってくる.
#今回は金利が返ってくる.
r_min = self.norminv(prob, er, stdev)
r_max = self.norminv(1-prob, er, stdev)
return r_min, r_max
#パラメーターはまとめて生成する.
def vasicek_diagonals(self, sigma, kappa, theta, r_min, dr, N, dtau):
#ようするにr.
#vaisicekで金利がどこに変化するかは分からないが、とりあえず自分の位置の金利はわかる.
rn = np.r_[0:N]*dr + r_min
subdiagonals = kappa*(theta-rn)*dtau/(2*dr) - 0.5*(sigma**2)*dtau/(dr**2)
diagonals = 1 + rn*dtau + sigma**2*dtau/(dr**2)
superdiagonals = -kappa*(theta-rn)*dtau/(2*dr) - 0.5*(sigma**2)*dtau/(dr**2)
# Implement boundary conditions.
if N > 0:
v_subd0 = subdiagonals[0]
superdiagonals[0] = superdiagonals[0] - subdiagonals[0]
diagonals[0] += 2*v_subd0
subdiagonals[0] = 0
if N > 1:
v_superd_last = superdiagonals[-1]
superdiagonals[-1] = superdiagonals[-1] - subdiagonals[-1]
diagonals[-1] += 2*v_superd_last
superdiagonals[-1] = 0
return subdiagonals, diagonals, superdiagonals
#行使するかしないか.
def check_exercise(self, V, eex):
return V > eex
def exercise_call_price(self, R, ratio, tau):
K = ratio*np.exp(-R*tau)
return K
def vasicek_policy_diagonals(self, subdiagonal, diagonal, superdiagonal, v_old, v_new, eex):
#行使されるかどうかTrue or False.
#すべての状態について.
has_early_exercise = self.check_exercise(v_new, eex)
#行使する状態について各iに関するdiagonalに0を入れる.
#これはf_i+1 = f_iとするためである.
subdiagonal[has_early_exercise] = 0
superdiagonal[has_early_exercise] = 0
policy = v_old/eex
policy_values = policy[has_early_exercise]
#行使された場合、価値はeex(K)となる。
#よって、bK = v_oldであり、その価値が最後まで引き継がれる.
#diagonal[has_early_exercise] = v_old/eex[has_early_exercise]
diagonal[has_early_exercise] = policy_values
#diagonal*eex = v_old
return subdiagonal, diagonal, superdiagonal
def iterate(self, subdiagonal, diagonal, superdiagonal, v_old, eex, max_policy_iter=10):
#このv_mplus1は不要.
#書いておくとわかりやすいためか.
v_mplus1 = v_old
v_m = v_old
#変化（行使）したか変化してないか. 変化したら1、してなければ0.
change = np.zeros(len(v_old))
prev_changes = np.zeros(len(v_old))
iterations = 0
while iterations <= max_policy_iter:
iterations += 1
#解く.
#v_oldは一期先の価格.
#これは繰り返すごとに変わる. diagonalsの変化によって.
v_mplus1 = self.tridiagonal_solve(subdiagonal, diagonal, superdiagonal, v_old)
#diagonalsを新しいv(v_mplus1)に応じて変える.
subdiagonal, diagonal, superdiagonal = self.vasicek_policy_diagonals(subdiagonal, diagonal, superdiagonal,
v_old, v_mplus1, eex)
#実際に計算上iterateしているのは上の二つ.
#ひたすらdiagonalsを更新している.
#v_oldは固定.
#diagonalsが更新されるのでv_mplus1は変わる.
#行使するかせざるか. True or False.
is_eex = self.check_exercise(v_mplus1, eex)
#行使するなら1にする.
change[is_eex] = 1
#2回目以上のiteration、すなわち1回目のiterationで行使が確認されたあと.
#変わっているところがないかチェック.
if iterations > 1:
change[v_mplus1 != v_m] = 1
is_no_more_eex = False if True in is_eex else True
#早期行使がされていないならやめる.
if is_no_more_eex:
break
#早期行使がなされた場合、された後の価値の変化を追う必要がある.
#価値に行使価格を代入する.
#ただし、これは計算には直接使わずcheckのために使う.
#最後の方に書いてもいい. ここに書くのはややややこしい.
v_mplus1[is_eex] = eex[is_eex]
#どの要素が変化したのかしてないのか.
#変化がなくなればすべて1になる（収束）
changes = (change == prev_changes)
#すべての値が変化していたらiterationをやめる.
is_no_further_changes = all((x == 1) for x in changes)
if is_no_further_changes:
break
prev_changes = change
v_m = v_mplus1
return v_mplus1, (iterations-1)
#Thomas algorism
def tridiagonal_solve(self, a, b, c, d):
nf = len(a) # Number of equations
#forループをしなくても各要素を演算.
#もともとがndarrayなので変化なし.
ac, bc, cc, dc = map(np.array, (a, b, c, d)) # Copy the array
#ac, bc, cc, dc = (a, b, c, d) # Copyでは無い場合.
for it in range(1, nf):
mc = ac[it]/bc[it-1]
bc[it] = bc[it] - mc*cc[it-1]
dc[it] = dc[it] - mc*dc[it-1]
xc = ac
xc[-1] = dc[-1]/bc[-1]
#状態を逆戻りして解く.
#最後にはdc[-1]/bc[-1]が入っている.
#それ以外に(dc[il]-cc[il]*xc[il+1])/bc[il]を入れる.
#xcは最初はaで、最後にrを返す.
for il in range(nf-2, -1, -1):
xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]
#ここはあっても無くてもいい.
del bc, cc, dc # Delete variables from memory
return xc
In [155]:
r0 = 0.05
R = 0.05
ratio = 0.95
sigma = 0.03
kappa = 0.15
theta = 0.05
prob = 1e-6
M = 250
max_policy_iter=10
grid_struct_interval = 0.25
rs = np.r_[0.0:2.0:0.1]
Vasicek = VasicekCZCB()
r, vals = Vasicek.vasicek_czcb_values(r0, R, ratio, 1., sigma, kappa, theta, M, prob,
max_policy_iter, grid_struct_interval, rs)
plt.figure(figsize=(10,5))
plt.title("Callable Zero Coupon Bond Values by r")
plt.plot(r, vals, label='1 yr')
for T in [5., 7., 10., 20.]:
r, vals = Vasicek.vasicek_czcb_values(r0, R, ratio, T, sigma, kappa, theta, M, prob,
max_policy_iter, grid_struct_interval, rs)
plt.plot(r, vals, label=str(T)+' yr', linestyle="--", marker=".")
plt.ylabel("Value ($)")
plt.xlabel("r")
plt.legend()
plt.grid(True)
plt.show()
上図はコール付きゼロクーポン債の、異なる満期と異なる金利に対する価値である。
1年満期のコール付きゼロクーポン債を考える。
この場合、行使価格はK（1年で割引）で固定されている.
あとはそれぞれの金利が与えられたときに、Kを上回るか上回らないかを考えればいい.
金利はvasicek modelで決まり、債権の価値は金利で決まる.
格子は縦軸が金利に、横軸が時間となる.
アメリカンなので早期行使を考えなくてはならない.
行使されたら価値はKになる.
In [ ]: