Calibrating a One-Factor Gaussian HJM Model to Forward Rates and Caps

# Calibrate 1-factor Gaussian HJM model on forward rates and ATM caps prices


# Calibrate 1-factor Gaussian HJM model on forward rates and ATM caps prices












I'm trying to solve the following problem as a part of the Interest Rate Models course

The algorithm that I'm following is

- derive simple rates from the given forward rates via $L(0, T_i) = \frac{(1+\frac{1}{2}F(0,T_{i-1},T_i))(1+L(0,T_{i-1})T_{i-1})-1}{T_i}$

- from the list of simple rates we can build the discount curve (calculate prices of zero-coupon bonds with corresponding maturities): $P(0,T_i) = \frac{1}{1+L(0,T_i)T_i}$

- calculate forward swap rates from estimated discount curve as $R_{swap}(T_0, T_n) = \frac{P(0,T_0)-P(0,T_n)}{0.5\sum_{i=1}^n P(0,T_i)}$, since given caps are ATM their respective strikes correspond to forward swap rates of the same maturity

- after we have strike rates for caps we are ready to invert the Black cap price formula $Cap(0,T_n) = \sum_{i=1}^n Cpl(T_{i-1}, T_i)$ (the price of the cap is the sum of the prices of all the caplets) with $Cpl(T_{i-1}, T_i) = 0.5P(0,T_i)(F(0,T_{i-1},T_i)Ф(d_{1,i})-\kappa Ф(d_{2,i}))$ and $d_{1,i}=\frac{\ln(F(0,T_{i-1},T_i)/\kappa)+\frac{\sigma^2}{2}T_{i-1}}{\sigma\sqrt{T_{i-1}}}$, $d_{2,i} = d_{1,i} - \sigma\sqrt{T_{i-1}}$, where $Ф$ is the standard normal cumulative distribution function, in order to obtain implied volatilities from the given market prices for four caps

- from implied volatilities we can calculate Black cap vegas which are just the derivatives of the price w.r.t. volatility parameter $\frac{\partial Cpl(T_{i-1},T_i)}{\partial\sigma}=0.5P(0,T_i)F(0,T_{i-1},T_i)\sqrt{T_{i-1}}\phi(d_{1,i})$ and $\frac{\partial Cap(0,T_n)}{\partial\sigma}=\sum_{i=1}^n\frac{\partial Cpl(T_{i-1},T_i)}{\partial\sigma}$ (the Black cap vega is the sum of all Black caplet vegas), where $\phi$ is the standard normal probability density function

- finally we can susbtitute the 1-factor Gaussian HJM volatility $\sigma(T) = \nu\exp^{-\beta T}$ into the Black cap price formula and solve the following minimization problem to calibrate the model parameters $\nu$ and $\beta$: $\min\limits_{\theta}\sum\limits_{n=1}^4\frac{1}{\left(\frac{\partial Cap_n}{\partial\sigma}\right)^2}\left(\hat{Cap_n^\theta}-Cap_n\right)^2$, where $\theta=(\nu,\beta)$ is the minimization parameter

This is my Python code which implements the aforementioned calibration algorithm:

```
import numpy as np
import math
import scipy
from scipy import optimize
from scipy import stats

Time = [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]
ForwardRates = [0.06, 0.08, 0.09, 0.10, 0.10, 0.10, 0.09, 0.09]
CapPrices = [0.002, 0.008, 0.012, 0.016]
SimpleRates = []
DiscountFactors = []

SimpleRates.append(ForwardRates[0])
DiscountFactors.append(1/(1+Time[0]*SimpleRates[0]))

for i in range(1, len(Time)):
    SimpleRates.append(((1+0.5*ForwardRates[i])*(1+Time[i-1]*SimpleRates[i-1])-1)/Time[i])
    DiscountFactors.append(1/(1+Time[i]*SimpleRates[i]))

SwapRates = [None]
for k in range(1, len(Time)):
    s = sum(DiscountFactors[i] for i in range(1, k + 1))
    SwapRates.append((DiscountFactors[0] - DiscountFactors[k])/(0.5*s))

def d1(i, strike, sigma):
    return (np.log(ForwardRates[i]/strike)+0.5*math.pow(sigma, 2)*Time[i-1])/(sigma*math.sqrt(Time[i-1]))

def d2(i, strike, sigma):
    return (np.log(ForwardRates[i]/strike)+0.5*math.pow(sigma, 2)*Time[i-1])/(sigma*math.sqrt(Time[i-1])) - sigma*math.sqrt(Time[i-1])

def caplet(begin, end, strike, sigma):
    return 0.5*DiscountFactors[end]*(ForwardRates[end]*scipy.stats.norm.cdf(d1(end, strike, sigma))-strike*scipy.stats.norm.cdf(d2(end, strike, sigma)))

def cap(end, sigma):
    strike = SwapRates[end]
    price = 0
    for j in range(1, end + 1):
        price += caplet(j - 1, j, strike, sigma)
    return price

ImpliedVolatilities = []
for j in range(1, 8, 2):
    ImpliedVolatilities.append(scipy.optimize.bisect(lambda sigma: cap(j, sigma) - CapPrices[int(0.5*(j-1))], 0.001, 0.5))

def capletVega(begin, end, strike, sigma):
    return 0.5*DiscountFactors[end]*ForwardRates[end]*math.sqrt(Time[begin])*scipy.stats.norm.pdf(d1(end, strike, sigma))

def capVega(end, sigma):
    strike = SwapRates[end]
    vega = 0
    for j in range(1, end + 1):
        vega += capletVega(j - 1, j, strike, sigma)
    return vega

vegas = []
for j in range(1, 8, 2):
    vegas.append(capVega(j, ImpliedVolatilities[int(0.5*(j-1))]))

def modelVolatility(v, b, time):
    return v*math.exp(-b*time)

def model_d1(i, strike, v, b):
    return (np.log(ForwardRates[i]/strike)+0.5*math.pow(modelVolatility(v, b, i), 2)*Time[i-1])/(modelVolatility(v, b, i)*math.sqrt(Time[i-1]))

def model_d2(i, strike, v, b):
    return (np.log(ForwardRates[i]/strike)+0.5*math.pow(modelVolatility(v, b, i), 2)*Time[i-1])/(modelVolatility(v, b, i)*math.sqrt(Time[i-1])) - modelVolatility(v, b, i)*math.sqrt(Time[i-1])

def modelCaplet(begin, end, strike, v, b):
    return 0.5*DiscountFactors[end]*(ForwardRates[end]*scipy.stats.norm.cdf(model_d1(end, strike, v, b))-strike*scipy.stats.norm.cdf(model_d2(end, strike, v, b)))

def modelCap(end, v, b):
    strike = SwapRates[end]
    price = 0
    for j in range(1, end + 1):
        price += modelCaplet(j - 1, j, strike, v, b)
    return price

expression = lambda theta: sum([math.pow(math.pow(vegas[int(0.5*(j-1))], 2),-1)*math.pow(modelCap(j, theta[0], theta[1]) - CapPrices[int(0.5*(j-1))], 2) for j in range(1, 8, 2)])
print("Implied volatilities are ")
print(ImpliedVolatilities)
print("Black vegas are ")
print(vegas)

theta = scipy.optimize.minimize(expression, [0.2, 0.2])
print(theta)

print(scipy.optimize.least_squares(expression, [0.2, 0.2]))
```

The implied volatilities that I got are

```
[0.19000807387027085, 0.13909831531679448, 0.09679878262167403, 0.09194895048610853]
```

The Black cap vegas are

```
[0.010510044022142772, 0.03762906503532581, 0.07393034105012006, 0.11751273646625979]
```

Running the scipy.optimize.minimize function gives me $\nu = 0.24239177, \beta = 0.25328877$ and running the scipy.optimize.least_squares gives $\nu = 0.24220154, \beta = 0.25285692$, however the automatic grader on Coursera does not accept neither $\nu=0.24$ nor $\beta=0.25$ as correct answers.

I would be extremely grateful if someone will find time to go through my algorithm (is my understanding of what I'm supposed to do in order to calibrate model correct?) and code (are there any errors in implementation?). Are my values for implied volatilies and vegas correct as an intermediate step? Any help will be appreciated.