# Bond Pricing with Vasicek Model

Author:<br>
Stanislav Khrapov<br>
<a href="mailto:khrapovs@gmail.com">khrapovs@gmail.com</a><br>
http://sites.google.com/site/khrapovs/<br>

## Introduction

The following code is the example of adapting methodology of<br>
<a href = "http://onlinelibrary.wiley.com/doi/10.1111/1468-0262.00164/abstract">Duffie, D., Pan, J., & Singleton, K. J. (2000). Transform Analysis and Asset Pricing for Affine Jump-Diffusions. Econometrica, 68(6), 1343–1376.</a>
<br>
to the Vasicek model of interest rates.

Set up environment


```
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


```
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
<p>
\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
<p>
\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*}


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


```
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:<p>
\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*}
<p>
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.


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


```
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)$


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


```
a = alpha_num(u, h, theta)
b = beta_num(u, h, theta)
sol_integrate = np.array([a, b])
```

### Compare solutions

Plot them.


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


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


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


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


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


```
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:
<p>
\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}
<p>
But in order to compute the bond price we obviously need to set $u=0$ in this transform:
<p>
\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}
<p>
In case of Vasicek model we need to solve the following system of ODE
<p>
\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
<p>
\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.


```
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.
<p>
\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}


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


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


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


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


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


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


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

baae5af6e65a4ec165129909f6ed9b7b00b44f1b
ipynb
Jupyter Notebook
Vasicek.ipynb
khrapovs/finmetrix-code
f278df1c15a225385846c2f0d7a6700c5737e901
MIT
2015-07-03T16:34:29.000Z
2019-05-09T13:10:26.000Z
Vasicek.ipynb
khrapovs/finmetrix-code
f278df1c15a225385846c2f0d7a6700c5737e901
MIT
Vasicek.ipynb
khrapovs/finmetrix-code
f278df1c15a225385846c2f0d7a6700c5737e901
MIT
2016-04-01T05:33:44.000Z
2020-07-12T06:58:25.000Z
Qwen/Qwen-72B
1. YES
2. YES