# Optimal Fitting Criteria of SABR

I was reading about SABR Model and curious about this. The process of fitting the SABR model involves finding values for the parameters α, β, ρ, ν that minimize the difference between model-implied option prices or volatilities and market-observed option prices or implied volatilities. This is typically done using numerical optimization techniques, such as the Levenberg-Marquardt algorithm or a least-squares optimization approach.

Can we (at least theoretically) use AIC/BIC criteria or else KL-Divergence for such optimisation problems? If yes can someone show me a python implementation of such? Or if not then can you please explain what are the problems and why in practice people doesn't use these?

import QuantLib as ql
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize
from scipy.optimize import  differential_evolution

moneyness=np.array([    120.00, 115.00, 114.00,     113.00,     112.00,     111.00 ,
110.00, 109.00 ,108.00,    107.00, 106.00, 105.00, 104.00, 103.00, 102.50, 102.00,
101.50, 101.00, 100.50, 100.00,     99.50,  99.00,  98.50,  98.00,  97.50,  97.00,
96.50,  96.00,  95.50,95.00,    94.50,  94.00,  93.50,  93.00,  92.00,  91.00,90.00 ])
moneyness=moneyness/100
fwd = 1.361794
strikes=fwd*moneyness
expiryTime = 30/365
marketVols = np.array([0.0302,  0.08827,    0.10019,    0.11206,    0.12343,    0.13305,    0.13968,
0.14325,  0.14431,    0.14337,    0.14077,    0.13706,    0.1325, 0.12664,
0.12355,  0.12024,    0.11694,    0.11371,    0.11062,    0.10772,    0.10526,
0.10307,  0.10118 ,0.09953,   0.09811,    0.0967, 0.09487,    0.09313,    0.09096,
0.08795,  0.08359,    0.07751,    0.07003,    0.06203,    0.04591,    0.02975,    0.01347 ])
marketVols=1-marketVols
marketVols[20]=1-10.97/100

def f(params):
params[0] = max(params[0], 1e-8)  # Avoid alpha going negative
params[1] = max(params[1], 1e-8)  # Avoid beta going negative
params[2] = max(params[2], 1e-8)  # Avoid nu going negative
params[3] = max(params[3], -0.999)  # Avoid rhp going < -1.0
params[3] = min(params[3], 0.999)  # Avoid rho going > 1.0

vols = np.array([
ql.sabrVolatility(strike, fwd, expiryTime, params[0],params[1],params[2],params[3])
#(Rate strike, Rate forward,  Time expiryTime,Real alpha,Real beta,Real nu,Real rho,

for strike in strikes
])
return ((vols - np.array(marketVols))**2 ).mean() **.5 #RSS Optimisation

bounds = [(0.000001,1000000),(0.000001,0.99) ,(0.00001, 1000000), (-0.99, 0.99)]
result = differential_evolution(f, bounds,tol=0.0000001)

params=result.x

print(1-ql.sabrVolatility(1.279, fwd, expiryTime, params[0],params[1],params[2],params[3]))
newVols = [ql.sabrVolatility(strike, fwd, expiryTime, params[0],params[1],params[2],params[3]) for strike in strikes]
plt.plot(strikes, marketVols, marker='o', label="market")
plt.plot(strikes, newVols, marker='o', label="SABR")
plt.legend()
plt.show()
$$`$$
• Can you pls update your question to elaborate how the Python code in your post is related to your questions? Commented Aug 29, 2023 at 12:42
• I was using this code with some synthetic data to calibrate the SABR model using the optimisation method RSS. You can see I've put the comment in the code where I was doing the RSS optimisation under the function 'f'. Commented Aug 30, 2023 at 3:28

In practice (at least in the rates world), $$\beta$$ is preset and $$\alpha$$ is solved for to calibrate to the atm vols $$\sigma_{ATM}$$ (which are the most liquid and reliable of the market data available). For instance, in the case for normal vols and assuming a normal distribution of the fwds, $$\beta=0$$ then $$σ_{N,ATM}=α\left(1+\frac{2−3ρ^2}{24}ν^2T\right).$$ The $$\rho,\nu$$ parameters are then obtained via the sort of optimization routines you describe in order to incorporate the skew.

Edit: I have added a simple SABR calibration routine (employing RSS) I use to illustrate what I mean below.

def sabr_calibration(swvolcube,a,b,calculation_date,alphas,tolerance):
print('for '+str(a.normalized())+str(b.normalized())+' optimizing alpha to fit sabr to atm vols ...')
T=SABR(swvolcube,a,b,calculation_date,'T')
f=SABR(swvolcube,a,b,calculation_date,'f')
atm_sabr_vol=SABR(swvolcube,a,b,calculation_date,'atm_sabr_vol')
atm_vol=SABR(swvolcube,a,b,calculation_date,'atm_vol')
atm_smile=SABR(swvolcube,a,b,calculation_date,'atm_smile')
beta=SABR(swvolcube,a,b,calculation_date,'beta')
nu=SABR(swvolcube,a,b,calculation_date,'nu')
rho=SABR(swvolcube,a,b,calculation_date,'rho')
cubic0=-atm_vol*(f**(-beta))
cubic1=(1+((2-3*(rho**2))/24)*(nu**2)*T)
cubic2=(rho*beta*nu*T)/(4*(f**(1-beta)))
cubic3=(beta*(beta-2)*T)/(24*(f**(2-2*beta)))
coeff=[cubic3,cubic2,cubic1,cubic0]
roots=[np.roots(coeff)[i] for i in range(0,len(np.roots(coeff)))]
positive_roots=[i for i in roots if i>0]
positive_roots.sort()
root=positive_roots[0]
dummy_alpha = alphas.loc[str(a.normalized()).lower(), str(b.normalized()).lower()]
dummy_nu = nus.loc[str(a.normalized()).lower(), str(b.normalized()).lower()]
dummy_rho = rhos.loc[str(a.normalized()).lower(), str(b.normalized()).lower()]
if abs(atm_vol-atm_sabr_vol)>tolerance:
dummy_alpha.setValue(root)
sabr_calibration(swvolcube, a, b, calculation_date,alphas,tolerance)
else:
print('atm vol = ' + str(atm_vol))
print('atm sabr vol = '+str(atm_sabr_vol))
print('optimized alpha = '+str(dummy_alpha.value()))
print('error is = ' + str(abs(atm_vol - atm_sabr_vol)))
print('calibrating nu and rho ...')
params = np.array([nu, rho])
def calib(params):
vols = np.array([
ql.sabrVolatility(strike, f, T, dummy_alpha.value(), beta, *params, vol_type)
for strike in atm_strike_set
])
return ((vols - np.array(atm_smile)) ** 2).mean() ** .5
cons = (
{'type': 'ineq', 'fun': lambda x: 0.999 + x[1]},
{'type': 'ineq', 'fun': lambda x: 0.999 - x[1]},
{'type': 'ineq', 'fun': lambda x: x[0] - 1e-15},
)
result = minimize(calib, params,constraints=cons)
new_params = result['x']
nu = new_params[0]
rho = new_params[1]
dummy_nu.setValue(nu)
dummy_rho.setValue(rho)
atm_sabr_vol = SABR(swvolcube, a, b, calculation_date, 'atm_sabr_vol')
atm_vol = SABR(swvolcube, a, b, calculation_date, 'atm_vol')
if abs(atm_vol - atm_sabr_vol) > tolerance:
print('re-calibrating ...')
sabr_calibration(swvolcube, a, b, calculation_date, alphas, tolerance)

return dummy_alpha.value(),nu,rho,beta
• Thank you for you detailed post and the approximate close form for computing the sigma. I'm actually more interested into knowing about the optimisation techniques, like is RSS (linear/non-linear) the best one or we could utilise AIC/BIC or KL-Divergence to get better optimisation result. Are those better/worse and if yes then what's the logic behind it? And if someone could guide me to implement it in some language (python). Commented Aug 30, 2023 at 3:33
• In my experience RSS is sufficient to obtain $\rho,\nu$ once $\alpha$ is solved for. I'm unfamiliar with the other optimization techniques you mention. The point I wanted to stress was there is a hierarchical importance in the parameter calibration methodology i.e. ATM vols should be an exact match and must be calibrated first - $\rho,\nu$ have little bearing on these. Their fitting is hence of secondary importance. Commented Aug 30, 2023 at 6:26
• I understand, thank you for taking time and answering this. Can you please refer me to some more guides (academic papers, blogs, maybe step by step guide to implement in market(python code) anything) where I can read and know about this more? Commented Aug 30, 2023 at 7:16
• I found uk.mathworks.com/help/fininst/calibrating-the-sabr-model.html very useful when learning about this stuff. It's not python but is very practical and hands-on nonetheless. Commented Aug 30, 2023 at 7:21
• Nice answer: +1 Commented Aug 30, 2023 at 7:50