SABR Calibration: Fit the ATM Level Before the Smile

# Optimal Fitting Criteria of SABR


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

## Answer by user35980 (score 2, accepted)

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

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