Stripping Cap Volatility into Forward Rate Volatilities

# Forward Rate Volatility Calculation - Caps


# Forward Rate Volatility Calculation - Caps












I am trying to calculate the forward rate volatilities from cap volatilities using Rebonato`s volatility model. Unfortunately, my approach always results in unrealistic forward rate vols. Furthermore, as I just have ATM Cap prices/vols, I cannot use the standard quantlib CapFloorVolTermSurface to construct my forward rate vols. Does anyone have a method, source, example in Python etc. for such a calculation?

This is the code I am using with some example data:

```
# Libraries
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
import matplotlib.pyplot as plt

# Input Parameters
cap_prices = [0.0, 0.0, 0.020011950921853195, 0.05074873822182878, 0.10336895891346086, 0.16588313145875605, 0.2533724574981232, 0.35174757325406103, 0.4701444936626773, 0.593081492686773, 0.7359526708319613, 0.877989696611585, 1.0260351088701913, 1.1600196162729612, 1.3184578716861968, 1.4829447531664126, 1.6496328932866386, 1.8065491133537495, 1.9820594340975828, 2.1490003466838026, 2.3225239850625234, 2.4917656955667002, 2.667050398665095, 2.8382619196071355, 3.0177156787135777, 3.189359811385567, 3.3707654678881247, 3.547975494866014, 3.730623411067357, 3.9023650797890843, 4.084721082037504, 4.261963691806112, 4.443314057979344, 4.613105915636499, 4.79282889899193, 4.967174224337874, 5.149200089709508, 5.310646042082241, 5.485811919729442, 5.659558511569442, 5.829614241194288, 5.994257618540332, 6.161495290757493, 6.323429110131764, 6.49006121048868, 6.647293059843137, 6.809224902486905, 6.96847370063799, 7.130463504868825, 7.281679545808297]
strikes = [0.004849792451814713, 0.005451685952049794, 0.0062512100930573805, 0.007193498719152967, 0.008261355914459145, 0.009416334324326169, 0.01055839873497409, 0.011689476404909173, 0.012754826632794203, 0.013766707730532313, 0.014703677769218961, 0.015590397988913039, 0.016469739452814425, 0.017311301978518785, 0.01804840121058648, 0.018723144744248053, 0.019403277949767234, 0.020055902491381006, 0.02061959943212093, 0.02115007947765411, 0.021641086946377633, 0.022081554325116624, 0.022475002908413484, 0.02285795251641165, 0.023222384445077595, 0.02355677277644328, 0.023849985574081144, 0.02409123389565899, 0.02427960601198556, 0.024453183806526894, 0.02461667740009192, 0.024769546952392436, 0.024911271573602314, 0.025041347690598202, 0.02515928730300895, 0.02526461610929479, 0.025356871485923523, 0.025435600304625477, 0.025501423796050745, 0.025562734092897478, 0.02562097641569784, 0.025675619242567206, 0.025726142034109705, 0.02577203490840359, 0.025812798169585514, 0.025847941687937193, 0.025876984128795716, 0.025899452027030454, 0.02591511856031686, 0.025928016839032297]
discount_factors = [0.99806973, 0.99565537, 0.99264992, 0.98876458, 0.98382339, 0.97767011, 0.97025005, 0.9617984, 0.95234658, 0.94217088, 0.9313256, 0.92001912, 0.90823235, 0.89578857, 0.88292001, 0.8701645, 0.85735653, 0.84405628, 0.83054739, 0.81741161, 0.80425875, 0.79119816, 0.77839357, 0.76587289, 0.75330698, 0.74080304, 0.72851779, 0.71661701, 0.70527458, 0.69455344, 0.68402033, 0.67362934, 0.66339853, 0.65334579, 0.64348887, 0.63384532, 0.62443248, 0.61526749, 0.6063673, 0.59773212, 0.58923836, 0.5808657, 0.57262467, 0.56452597, 0.55658046, 0.54879915, 0.54119311, 0.53377354, 0.52655172, 0.51953459, 0.51265042]
T = np.arange(0.5, 25.5, 0.5)

# Code for optimization

# Bacheliermodel for normal vol
def bachelier_caplet_price(F, K, sigma, T, df):
    d1 = (F - K) / (sigma * np.sqrt(T))
    caplet_price = df * ((F-K) * norm.cdf(d1) + sigma * np.sqrt(T) * norm.pdf(d1))
    return caplet_price

# Rebonato formula
def rebonato_volatility(T, alpha, beta, gamma, delta):
    return (alpha + beta * T) * np.exp(-gamma * T) + delta

# Objective Function
def objective_function(params, cap_prices, forward_swap_rates, discount_factors, T):
    alpha, beta, gamma, delta = params
    sigma = rebonato_volatility(T, alpha, beta, gamma, delta)
    model_prices = np.array([
        bachelier_caplet_price(forward_swap_rates[i], forward_swap_rates[i], T[i], sigma[i], discount_factors[i])
        for i in range(len(forward_swap_rates))
    ])
    return np.sum((model_prices - cap_prices) ** 2)

# Optimizer
initial_params = [0.1, 0.1, 0.1, 0.1]
result = minimize(objective_function, initial_params, args=(cap_prices, strikes, discount_factors, T), method='SLSQP')
optimized_params = result.x
alpha, beta, gamma, delta = optimized_params

# Calculation of forward rate vols
forward_rate_vols = np.array([rebonato_volatility(t, alpha, beta, gamma, delta) for t in T])
```

The results I get are obviously wrong, as they should be closer to the implied cap vols and not just increasing.

The implied cap vols look like this:

```
implied_cap_vols = {'1': 0.20199999999999999, '2': 0.366, '3': 0.496, '4': 0.609, '5': 0.693, '6': 0.7509999999999999, '7': 0.7709999999999999, '8': 0.7909999999999999, '9': 0.802, '10': 0.807, '12': 0.8059999999999999, '15': 0.794, '20': 0.7609999999999999, '25': 0.723}
```

Thank you very much for any help on this!

## Answer by Hasek (score 2)

https://quant.stackexchange.com/a/78454

What you are asking for is called volatility stripping.

A cap is a series of consecutive caplets which are call options on forward rates. Caplets are not traded individually on the market so there are no quoted prices or volatilities for them and thus you need to somehow infer them from cap prices or volatilities. Note that since there are several caplets in a cap, there is no unique solution.

There are several common ways to approach this problem. The easiest one is probably to assume a piecewise-constant volatility term structure:

- sort the caps by maturity in ascending order

- calculate price differences for all consecutive caps

- partition the caplets in price differences by maturities

- for each caplet partition solve for the volatility which allows to reprice the price difference by some root finding algorithm like bisection or Newton's method

This approach should be pretty stable and easy to implement.

I suggest you to take a look at OpenGamma Caplet Stripping paper in case you're looking for something more sophisticated. It isn't easy to say what is going wrong with your numerical optimization without seeing your code and input data. I'm not entirely sure either what did you mean by a parametric volatility model. Are you talking about something like Rebonato's volaitlity parametrization for the LIBOR Market Model $$\sigma^{model}(t) = (a + b\cdot t)\cdot e^{-c\cdot t} + d$$ optimized to give the best solution for the following MSE problem? $$\min\limits_{a,b,c,d}\sum\limits_{K}\left(\sigma^{model}(K,t,a,b,c,d) - \sigma^{market}(K,t)\right)^2$$

I suggest you to start with the QuantLib built-in functionality. Here is an example of how you can strip a given cap volatility surface.

```
import QuantLib as ql
import numpy as np
import matplotlib.pyplot as plt

calc_date = ql.Date(23, 2, 2024)
ql.Settings.instance().evaluationDate = calc_date

dates = [calc_date, ql.Date(25, 6, 2024), ql.Date(24, 2, 2025), 
         ql.Date(24, 2, 2026), ql.Date(24, 2, 2027), ql.Date(24, 2, 2028), 
         ql.Date(24, 2, 2029), ql.Date(24, 2, 2030)]
yields = [0.07, 0.07, 0.065, 0.05, 0.052, 0.055, 0.057, 0.063]

dcc = ql.Actual365Fixed()
calendar = ql.NullCalendar()
interpolation = ql.Linear()
compounding = ql.Compounded
compounding_frequency = ql.Annual
term_structure = ql.ZeroCurve(dates, yields, dcc, calendar, interpolation, 
                              compounding, compounding_frequency)
ts_handle = ql.YieldTermStructureHandle(term_structure)

start_date = calc_date
end_date = start_date + ql.Period(6, ql.Years)
period = ql.Period(3, ql.Months)
bdc = ql.ModifiedFollowing
rule = ql.DateGeneration.Forward
end_of_month = False
schedule = ql.Schedule(start_date, end_date, period, calendar, bdc, bdc, 
                       rule, end_of_month)
ts_handle = ql.YieldTermStructureHandle(term_structure)

ibor_index = ql.USDLibor(ql.Period(3, ql.Months), ts_handle)
ibor_index.addFixing(ql.Date(23, 2, 2024), 0.07)
notional = 1000000
ibor_leg = ql.IborLeg([notional], schedule, ibor_index)

strikes = [0.06, 0.07, 0.08]
expiries = [ql.Period(i, ql.Years) for i in range(1, 6)]
vols = ql.Matrix(len(expiries), len(strikes))
data = [[47.27, 55.47, 64.07, 58.71, 59.90],
        [45.65, 54.15, 61.47, 55.55, 58.23],
        [46.69, 52.65, 59.32, 54.29, 57.40]]
for i in range(vols.rows()):
    for j in range(vols.columns()):
        vols[i][j] = data[j][i] / 100.0

capfloor_surf = ql.CapFloorTermVolSurface(2, calendar, bdc,
                                          expiries, strikes, vols, dcc)

optionlet_surf = ql.OptionletStripper1(capfloor_surf, ibor_index, 
                                       ql.nullDouble(), 1e-6, 100)
ovs_handle = ql.OptionletVolatilityStructureHandle(
    ql.StrippedOptionletAdapter(optionlet_surf))

tenors = np.arange(1/360, 4.75, 0.25)
strike = 0.07
optionlet_vols = [ovs_handle.volatility(t, strike) for t in tenors]
vols = [100 * vol for vol in optionlet_vols]
plt.plot(tenors, vols, "-", 
         label="Caplet {0}% Stripped Vols".format(round(strike * 100)))
plt.legend(loc="lower center")
```

This will be the result for that (completely random) data