Bond Pricing with Vasicek Model

Introduction

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

The Model

\begin{equation} dr_{t}=\kappa\left(\mu-r_{t}\right)dt+\sigma dW_{t}. \end{equation}

Set up parameters


In [2]:
h = 1
kappa, mu, sigma = 1.5, .5, .1
theta = [kappa, mu, sigma]
u = np.linspace(0, 1e3, 1e2)

Characteristic function

Analytic solution

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

Numeric solution

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

Compare solutions

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


Compare characteristic functions

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


Compare transition densities

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)


CPU times: user 8.27 ms, sys: 7 µs, total: 8.27 ms
Wall time: 7.86 ms
CPU times: user 40.8 s, sys: 15.9 ms, total: 40.8 s
Wall time: 40.8 s

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


Bond Pricing

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$.

Analytic solution

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

Numeric solution

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

Comparison

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