Skip to main content
Quantitative Finance

Return to Answer

Adding used code
Source Link
user35980
user35980
  • 1.5k
  • 7
  • 7

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

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
Source Link
user35980
user35980
  • 1.5k
  • 7
  • 7

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.