# ATM cap prices in Vasicek model (Filipovic)

I am trying to replicate the ATM cap prices in table 7.1 (see bottom of this post) from Filipovic's book "Term Structure Models - A Graduate Course" which assume the Vasicek model and uses Black's formula. It is assumed that $$t_0 = 0$$ and the tenor is quarterly, i.e. $$\delta := T_i - T_{i-1} = 1/4$$ for $$i=1,...,119$$. Filipovic uses the paramerization $$dr(t) = (b-\beta r(t))dt + \sigma dW(t),$$

for the Vasicek model, while I've used the parameterization

$$dr(t) = a(b-r(t))dt + \sigma dW(t).$$

It is furthermore assumed that (my parameterization in parentheses) $$\beta = -0.86 \; (a = 0.86), \quad b/|\beta| = 0.09 \; (b = 0.09), \quad \sigma = 0.0148, \quad r(0) = 0.08.$$

Using the Affine Term Structure (ATS) I find that the ZCB prices can be found using $$P(t,T) = \exp\{-A(t,T) - B(t,T)\cdot r(t)\}$$, where

$$A(t,T) = (b - \frac{\sigma^2}{2\cdot a})\cdot (T-t - B(t,T)) - \frac{\sigma^2 \cdot B(t,T)^2}{4\cdot a},$$

$$B(t,T) = \frac{1}{a} (1-e^{-a \cdot (T-t)}).$$

Hence, I can calculate the simple forward (LIBOR) rates as $$F(t,T,T+\delta) = \frac{1}{\delta}\left(\frac{P(t,T)}{P(t,T+\delta)} - 1\right).$$

A cap is said to be ATM if its strike $$K$$ equals the swap rate $$R(0) = \frac{P(t,T_0) - P(t,T_n)}{\delta \sum_{i=1}^n P(t,T_i)}.$$

I use the Black's formula to price the caplets $$Cpl(t; T_{i-1}, T_i) = \delta P(t,T_i) \left[ F(t;T_{i-1}, T_i) N(d_+(i;t)) - K N(d_-(i;t))\right],$$ where $$d_\pm = \frac{\log\left(\frac{F(t,T_{i-1},T)}{K}\right) \pm \frac12 \sigma(t)^2 (T_{i-1}-t)}{\sigma(t) \sqrt{T_{i-1}-t}}$$ and $$N$$ denotes the Gaussian CDF and $$\sigma(t)$$ denotes the cap (implied) volatility.

Finally, I use that the price of the cap is the sum of the caplets.

I've tried price the first cap (maturity 1) in the table manually using the code below, but I do not get the correct price.

import numpy as np
from scipy.stats import norm

N = lambda x: norm(loc=0.0, scale=1.0).cdf(x)

class Vasicek:
def __init__(self, a, b, sigma):
self.a = a
self.b = b
self.sigma = sigma

def _calc_A(self, t):
B = self._calc_B(t)
return (self.b - self.sigma**2 / (2*self.a**2)) * (t - B) - \
self.sigma ** 2 * B ** 2 / (4 * self.a)

def _calc_B(self, t):
return  (1 - np.exp(-self.a * t)) / self.a

def calc_zcb(self, t, r0):
return np.exp(-self._calc_A(t) - self._calc_B(t) * r0)

# Parameters
a = np.array(0.86)
b = np.array(0.09)
sigma = np.array(0.0148)
r0 = np.array(0.08)
dt = np.array(0.25)

# Make model
mld = Vasicek(a=a, b=b, sigma=sigma)

# Zero coupon bonds
zcb_1 = mld.calc_zcb(t=np.array(0.25), r0=r0)
zcb_2 = mld.calc_zcb(t=np.array(0.50), r0=r0)
zcb_3 = mld.calc_zcb(t=np.array(0.75), r0=r0)
zcb_4 = mld.calc_zcb(t=np.array(1.00), r0=r0)

# Swap rate and simple forward rate
swap_rate = (zcb_1 - zcb_4) / (0.25 * (zcb_2 + zcb_3 + zcb_4))
fwd_1_2 = (zcb_1 / zcb_2 - 1) / 0.25
fwd_2_3 = (zcb_2 / zcb_3 - 1) / 0.25
fwd_3_4 = (zcb_3 / zcb_4 - 1) / 0.25

# Caplet prices
d1_1_2 = (np.log(fwd_1_2 / swap_rate) + 0.5 * sigma * 0.25) / (sigma * 0.25)
d2_1_2 = (np.log(fwd_1_2 / swap_rate) - 0.5 * sigma * 0.25) / (sigma * 0.25)
d1_2_3 = (np.log(fwd_2_3 / swap_rate) + 0.5 * sigma * 0.50) / (sigma * 0.50)
d2_2_3 = (np.log(fwd_2_3 / swap_rate) - 0.5 * sigma * 0.50) / (sigma * 0.50)
d1_3_4 = (np.log(fwd_3_4 / swap_rate) + 0.5 * sigma * 0.75) / (sigma * 0.75)
d2_3_4 = (np.log(fwd_3_4 / swap_rate) - 0.5 * sigma * 0.75) / (sigma * 0.75)

cpl_1_2 = 0.25 * zcb_2 * (fwd_1_2 * N(d1_1_2) - swap_rate * N(d2_1_2))
cpl_2_3 = 0.25 * zcb_3 * (fwd_2_3 * N(d1_2_3) - swap_rate * N(d2_2_3))
cpl_3_4 = 0.25 * zcb_4 * (fwd_3_4 * N(d1_3_4) - swap_rate * N(d2_3_4))

# Cap
cap = cpl_1_2 + cpl_2_3 + cpl_3_4

print(cap) # 0.011209781683249345