Author:
Stanislav Khrapov
khrapovs@gmail.com
http://sites.google.com/site/khrapovs/
The following code is the example of adapting methodology of
Duffie, D., Pan, J., & Singleton, K. J. (2000). Transform Analysis and Asset Pricing for Affine Jump-Diffusions. Econometrica, 68(6), 1343–1376.
to the Vasicek model of interest rates.
Set up environment
In [1]:
import scipy.integrate as si
import numpy as np
import matplotlib.pylab as plt
import functools
Set up parameters
In [2]:
h = 1
kappa, mu, sigma = 1.5, .5, .1
theta = [kappa, mu, sigma]
u = np.linspace(0, 1e3, 1e2)
The analytic solution for the characteristic function is exponentially affine
\begin{equation} E_{t}\left[\exp\left(iur_{T}\right)\right]=\exp\left(\alpha_{t}\left(u\right)+\beta_{t}\left(u\right)r_{t}\right), \end{equation} with
\begin{eqnarray*} \beta_{t}\left(u\right) & = & iue^{-\kappa h}\\ \alpha_{t}\left(u\right) & = & \mu iu\left(1-e^{-\kappa h}\right)-\frac{\sigma^{2}}{4\kappa}u^{2}\left(1-e^{-2\kappa h}\right) \end{eqnarray*}
In [3]:
def alpha(u, h, theta):
kappa, mu, sigma = theta
e = np.exp( - kappa * h )
return 1j * u * mu * (1 - e) - u ** 2 * sigma ** 2 * (1 - e ** 2) / kappa / 4
def beta(u, h, theta):
kappa, mu, sigma = theta
return 1j * u * np.exp( - kappa * h )
Evaluate these functions on the grid.
In [4]:
a = alpha(u, h, theta)
b = beta(u, h, theta)
sol_analytic = np.array([a, b])
We need to solve the following system of ODE:
\begin{eqnarray*} \dot{\alpha}_{t} & = & -\kappa\mu\beta_{t}-\frac{1}{2}\sigma^{2}\beta_{t}^{2},\\ \dot{\beta}_{t} & = & \kappa\beta_{t}, \end{eqnarray*}
with terminal conditions $\beta_{T}=iu$ and $\alpha_{T}=0$.
The following function is the right hand side of the ODE system. Since the solution is actually in terms of distance $h$ to the terminal date and not in calendar time $t$, the left-hand side derivatives go with minus.
In [5]:
def df(t, x, theta):
kappa, mu, sigma = theta
da = - kappa * mu * x[1] - .5 * sigma ** 2 * x[1] ** 2
db = kappa * x[1]
return [-da, -db]
This function solves the system of complex valued ODE.
In [6]:
def f(df, u, h, theta):
dt = h / 1e1
df_args = functools.partial(df, theta = theta)
r = si.complex_ode(df_args)
sol = []
for v in u:
x0 = [0., v * 1j]
r.set_initial_value(x0, 0)
while r.successful() and r.t < h:
r.integrate(r.t + dt)
sol.append(r.y)
return np.array(sol).T
Separate solution into functions $\alpha_t(u)$ and $\beta_t(u)$
In [7]:
alpha_num = lambda u, h, theta: f(df, u, h, theta)[0]
beta_num = lambda u, h, theta: f(df, u, h, theta)[1]
Evaluate numerical solutions at specific values.
In [8]:
a = alpha_num(u, h, theta)
b = beta_num(u, h, theta)
sol_integrate = np.array([a, b])
Plot them.
In [9]:
fig, axes = plt.subplots(2, 2, figsize = (8, 6), sharex = True)
for sol in [sol_analytic, sol_integrate]:
axes[0,0].plot(u, sol[0].real)
axes[1,0].plot(u, sol[0].imag)
axes[0,1].plot(u, sol[1].real)
axes[1,1].plot(u, sol[1].imag)
axes[0,0].set_title('Alpha real')
axes[1,0].set_title('Alpha imag')
axes[0,1].set_title('Beta real')
axes[1,1].set_title('Beta imag')
axes[1,0].set_xlabel('u')
axes[1,1].set_xlabel('u')
axes[0,1].legend(['Analytic','Numeric'])
plt.show()
Define conditional characteristic functions.
In [10]:
rt = 1.5
psi_analytic = lambda u: np.exp( alpha(u, h, theta) + beta(u, h, theta) * rt )
psi_numeric = lambda u: np.exp( alpha_num(u, h, theta) + beta_num(u, h, theta) * rt )
Plot them.
In [11]:
x = np.linspace(0, 60, 1e2)
psi_a = psi_analytic(x)
psi_n = psi_numeric(x)
fig, axes = plt.subplots(2, 1, figsize = (8, 6), sharex = True)
for psi in [psi_a, psi_n]:
axes[0].plot(x, psi.real)
axes[1].plot(x, psi.imag)
axes[1].set_xlabel('u')
axes[0].set_title('Real')
axes[1].set_title('Imag')
axes[0].legend(['Analytic','Numeric'])
plt.show()
Fast Fourier inverse of charcteristic function. Returns density.
In [12]:
def CFinverse(psi, A = -1e2, B = 1e2, N = 1e5):
eta = (N - 1) / N * 2 * np.pi / (B - A)
lmbd = (B - A) / (N - 1)
k = np.arange(0, N)
x = A + lmbd * k
v = eta * k
y = psi(v) * np.exp(- 1j * A * v) * eta / np.pi
f = np.fft.fft(y)
return x, f.real
Compute densities. Note how much time was spent on each operation.
In [13]:
N = 1e4
%time x, density_analytic = CFinverse(psi_analytic, N = N)
%time x, density_numeric = CFinverse(psi_numeric, N = N)
Plot transition densities.
In [14]:
plt.plot(x, density_analytic)
plt.plot(x, density_numeric)
plt.xlim([0, 2])
plt.legend(['Analytic','Numeric'])
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()
Thanks to Duffie, Pan, & Singleton (2000, Econometrica) we know how to compute the following transform for AJD models:
\begin{equation} \mathbb{T}^{\chi}\left(u,Y_{t},t,T\right)=E_{t}\left[\exp\left(-\int_{t}^{T}r_{s}ds\right)\exp\left(ur_{T}\right)\right] =\exp\left(\alpha_{t}(u)+\beta_{t}(u)r_{t}\right). \end{equation}
But in order to compute the bond price we obviously need to set $u=0$ in this transform:
\begin{equation} B_{t,T}=E_{t}\left[\exp\left(-\int_{t}^{T}r_{s}ds\right)\right]=\mathbb{T}^{\chi}\left(0,r_{t},t,T\right)=\exp\left(\alpha_{t}\left(0\right)+\beta_{t}\left(0\right)r_{t}\right), \end{equation}
In case of Vasicek model we need to solve the following system of ODE
\begin{eqnarray*} \dot{\beta}_{t} & = & 1+\kappa\beta_{t},\\ \dot{\alpha}_{t} & = & -\kappa\mu\beta_{t}-\frac{1}{2}\sigma^{2}\beta_{t}^{2}, \end{eqnarray*} with terminal conditions $\beta_{T}=u$ and $\alpha_{T}=0$.
The analytic solution to this problem is known to be
\begin{eqnarray*} \beta_{t}\left(0\right) & = & \frac{1}{\kappa}\left(e^{-\kappa h}-1\right),\\ \alpha_{t}\left(0\right) & = & -\left(\mu-\frac{\sigma^{2}}{2\kappa^{2}}\right)\left(\beta_{t}\left(0\right)+h\right)-\frac{\sigma^{2}}{4\kappa}\beta_{t}^{2}\left(0\right). \end{eqnarray*} It is coded below.
In [15]:
def beta(h, theta):
kappa, mu, sigma = theta
return ( np.exp( - kappa * h ) - 1 ) / kappa
def alpha(h, theta):
kappa, mu, sigma = theta
a = - (mu - .5 * sigma ** 2 / kappa ** 2) * (beta(h, theta) + h)
a -= .25 * sigma ** 2 / kappa * beta(h, theta) ** 2
return a
Most of the time we are actually interested in yields, rather than in bond prices directly.
\begin{equation} y_{t,T}= -\frac{1}{h}\log B_{t,T}=-\frac{1}{h}\alpha_{t}\left(0\right)-\frac{1}{h}\beta_{t}\left(0\right)r_{t} =A\left(h\right) + B\left(h\right)r_t \end{equation}
In [16]:
A_exact = lambda h, theta: - alpha(h, theta) / h
B_exact = lambda h, theta: - beta(h, theta) / h
We can compute these for the interval of $H$ days in case of analytic solution.
In [17]:
H = np.linspace(.1, 10, 1e2)
a = A_exact(H, theta)
b = B_exact(H, theta)
sol_analytic = np.array([a, b])
For numerical solution we need to set up the system of ODEs.
In [18]:
def df(t, x, theta):
kappa, mu, sigma = theta
da = - kappa * mu * x[1] - .5 * sigma ** 2 * x[1] ** 2
db = 1 + kappa * x[1]
return [-da, -db]
The solver is adjusted since for bond pricign we only need the value of functions $\alpha_t(u)$ and $\beta_t(u)$ at zero but for variety of $h$.
In [19]:
def f(df, H, theta):
df_args = functools.partial(df, theta = theta)
r = si.ode(df_args)
sol = []
for h in H:
dt = h / 1e2
x0 = [0, 0]
r.set_initial_value(x0, 0)
while r.successful() and r.t < h:
r.integrate(r.t + dt)
sol.append(r.y)
return np.array(sol).T
The solutions are written as function below.
In [20]:
alpha_num = lambda H, theta: f(df, H, theta)[0]
beta_num = lambda H, theta: f(df, H, theta)[1]
A_numeric = lambda H, theta: - alpha_num(H, theta) / H
B_numeric = lambda H, theta: - beta_num(H, theta) / H
We can evaluate these at a specific time interval.
In [21]:
a = A_numeric(H, theta)
b = B_numeric(H, theta)
sol_integrate = np.array([a, b])
Now plot the solutions and the resulting yield.
In [22]:
fig, axes = plt.subplots(3, 1, figsize = (8, 6), sharex = True)
for sol in [sol_analytic, sol_integrate]:
axes[0].plot(H, sol[0])
axes[1].plot(H, sol[1])
axes[2].plot(H, sol[0] + sol[1] * rt)
axes[0].set_title('A')
axes[1].set_title('B')
axes[2].set_title('Yields conditional on $r_0=$ ' + str(rt))
axes[2].legend(['Analytic','Numeric'])
axes[2].set_xlabel('H')
axes[2].set_ylabel('%')
plt.show()