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 ])
fwd = 1.361794
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 ])

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)


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")
  • $\begingroup$ Can you pls update your question to elaborate how the Python code in your post is related to your questions? $\endgroup$
    – Alper
    Aug 29, 2023 at 12:42
  • $\begingroup$ 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'. $\endgroup$
    – Starlord22
    Aug 30, 2023 at 3:28

1 Answer 1


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 ...')
atm_strike_set=[f+i for i in strike_spreads]
roots=[np.roots(coeff)[i] for i in range(0,len(np.roots(coeff)))]
positive_roots=[i for i in roots if i>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:
    sabr_calibration(swvolcube, a, b, calculation_date,alphas,tolerance)
    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]
    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
  • $\begingroup$ 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). $\endgroup$
    – Starlord22
    Aug 30, 2023 at 3:33
  • $\begingroup$ 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. $\endgroup$
    – user35980
    Aug 30, 2023 at 6:26
  • $\begingroup$ 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? $\endgroup$
    – Starlord22
    Aug 30, 2023 at 7:16
  • $\begingroup$ 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. $\endgroup$
    – user35980
    Aug 30, 2023 at 7:21
  • $\begingroup$ Nice answer: +1 $\endgroup$ Aug 30, 2023 at 7:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.