Optimal Fitting Criteria of SABR

Question

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()
```

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 ...')
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')
strike_spreads=SABR(swvolcube,a,b,calculation_date,'strike_spreads')
atm_strike_set=[f+i for i in strike_spreads]
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