# Fitting a volatility smile with pySABR -- Python implementation of SABR model

In order to model some volatility smiles I'm using the python's pySABR package.

I ran into a situation when I have two almost identical pieces of code for two different volatility smiles missing the ATM quotes and the pySABR can properly fit the ATM volatility in one case and can't in another.

This is the case when everything is working just fine:

import pysabr
import numpy as np
from pysabr import Hagan2002LognormalSABR as sabr
from pysabr import hagan_2002_lognormal_sabr as hagan2002

testStrikes = np.array([0.04, 0.06, 0.08, 0.10])
testVols = np.array([23.52, 16.24, 20.17, 26.19])
forward_3m_6m = (1/0.25) * (-1 + (1+0.0753*0.5) / (1+0.0747*0.25))

calibration = sabr(f = forward_3m_6m, shift = 0, t = 0.5, beta = 0.5).fit(testStrikes, testVols)
smile = []
test_smile = []
for strike in testStrikes:
smile.append(sabr(f = forward_3m_6m, shift = 0, t = 0.5, v_atm_n = 136.75/10000.00, beta = 0.5, rho = calibration[1], volvol = calibration[2]).lognormal_vol(strike) * 100.00)
test_smile.append(hagan2002.lognormal_vol(strike, forward_3m_6m, 0.5, calibration[0], beta, calibration[1], calibration[2]) * 100.00)

print(smile)
print(test_smile)
print(hagan2002.lognormal_vol(k = 0.0745136, f = forward_3m_6m, t = 0.5, alpha = calibration[0], beta = 0.5, rho = calibration[1], volvol = calibration[2]) * 100.00)


The output is

[23.52579980722943, 16.22619971530687, 20.186954023608315, 26.176954813043512]

[23.52656681608356, 16.227190950406076, 20.188104613955648, 26.178058303454062]

18.369296025878036


The difference between the first and the second volatility lists is that the first was built using the ATM normal volatility quote instead of the parameter alpha.

This is the code for the second case when I ran into a problem with different results for different methods and inadequate value of the resulting ATM lognormal volatility:

import pysabr
import numpy as np
from pysabr import Hagan2002LognormalSABR as LNsabr
from pysabr import hagan_2002_lognormal_sabr as hagan2002LN

strikes = np.array([0.05, 0.055, 0.06, 0.0650, 0.07, 0.08, 0.0850, 0.09, 0.095, 0.10])
LogNormalVols = np.array([18.90, 17.30, 16.34, 16.29, 17.19, 20.29, 21.89, 23.42, 24.84, 26.16])
vol_n_ATM = 141.01 / 10000.00
forward_3m_6m = (1 / 0.25) * (- 1 + ( 1+ 0.0756 * 0.5) / (1 + 0.0746 * 0.25))
discount_0m_6m = 0.0753

beta = 0.5
calibration_LN = LNsabr(forward_3m_6m, 0, 0.5, beta).fit(strikes, LogNormalVols)
modelVols_LN = []
test_LN = []
for strike in strikes:
modelVols_LN.append(LNsabr(forward_3m_6m, 0, 0.5, vol_n_ATM, beta, calibration_LN[1], calibration_LN[2]).lognormal_vol(strike) * 100.00)
test_LN.append(hagan2002LN.lognormal_vol(strike, f = forward_3m_6m, t = 0.5, alpha = calibration_LN[0], beta = beta, rho = calibration_LN[1], volvol = calibration_LN[2]) * 100.00)

print(modelVols_LN)
print(test_LN)
print(hagan2002LN.lognormal_vol(k = 0.0753, f = forward_3m_6m, t = 0.5, alpha = calibration_LN[0], beta = beta, rho = calibration_LN[1], volvol = calibration_LN[2]) * 100.00)


The output is

[20.14451199912703, 18.54257849585578, 17.499371075753768, 17.1951750118971, 17.677189811525658, 20.011518064279755, 21.365538860352675, 22.691679608999, 23.95616514380161, 25.148945356594965]

[68.433853990843, 67.98546902874503, 67.8392636042873, 67.91509150013229, 68.15026017672005, 68.91958756908917, 69.39185669606091, 69.89474554951227, 70.41454229587038, 70.94145252003263]

68.64132679555016


One can easily see that reproducing the volatility smile using all four of the calibrated parameters leads to a way higher lognormal volatilities than expected.

What am I doing wrong? Where is the error in the second code? Any help will be appreciated.

• what is the output of the fit? i.e. your calibration_LN object.
– will
Oct 2 at 10:58
• @will, the output is [0.18123299730498071, 0.5465080547109832, 0.8266797808547465], here the first element is alpha, second -- rho, third -- volvol. Note that both methods producing lognormal volatilities use the same calibration parameters. Oct 2 at 17:02
• Sorry, not what I meant. Normally, fitting functions output an object that tells you whether or not the fit was successful, the number or iterations, function evaluations, the fitting error, etc.
– will
Oct 3 at 12:11

## Edit: You need to specify the keyword arguments in your second example

The poor calibration of your second example comes from the fact, that you didn't define the keyword arguments in the LNsabr-object before the .fit() function.

calibration_LN = LNsabr(forward_3m_6m, 0, 0.5, beta).fit(strikes, LogNormalVols)


You only need to make the slight change of defining the keyword arguments:

calibration_LN = LNsabr(f = forward_3m_6m, shift = 0, t = 0.5, beta = beta).fit(strikes, LogNormalVols)


This results in a good fit:

I have extended my current appendix to contain a slightly rewritten version of your Python code provided above. You can run this yourself and see whether you get the same calibrated parameters. I am not completely apparent why the calibration becomes poor when not specifying the keyword arguments. I hope this resolves your issue. Please provide some feedback when you have time.

I will leave my old answer below for future readers to see.

## Old answer: There might be an error with the optimization routine in the pysabr package

The bad calibration might come from the optimization routine in the Pysabr package hitting a local minimum or similar.

I reimplemented your second example (which caused issues) in MATLAB using their inbuilt function of the same SABR lognormal polynomial expansion which is also provided in the pysabr package in Python (My code is provided in an appendix at the end of my answer.). The MATLAB documentation can be found here, and a clear-cut calibration example of the SABR model is documented here. In the latter, they walk through both ways of calibrating the SABR model: First, calibrating $$\alpha, \rho, \nu$$ directly and lastly calibrating $$\rho$$ and $$\nu$$ by implying $$\alpha$$ from the ATM volatility. Following the first calibration example in MATLABs documentation with your specified vols and strikes, I get the calibrated parameters to:

$$\left[\alpha, \rho, \nu\right] = \left[0.0504, \: 0.6407, \:0.8797\right],$$

where $$\nu$$ is the vol of vol parameter. If you insert these estimates in the calibration_LN list object, you will observe a very good fit for your SABR log-normal model, when calibrating $$\alpha$$ freely:

From this, I believe the optimization routine in the pysabr package might have some issues when calibrating the smile. I will update my answer if I find the reason why the package produces such a bad fit. As for now, this is the only help I can provide.

## Appendix:

### Matlab code

%Matlabs inbuilt functions uses settlementdate and exercisedate to
%calculate T internally. These are just chosen to get T=0.5:
Settle = '01-Jan-2020';
ExerciseDate = '01-Jul-2020';

%Checking if settle and exercise date becomes T=0.5.
T = yearfrac(Settle, ExerciseDate, 1);

Beta1 = 0.5;

%Following the matlab guide:
MarketStrikes = [0.05, 0.055, 0.06, 0.0650, 0.07, 0.08, 0.0850, 0.09, 0.095, 0.10]';
MarketVolatilities = [18.90, 17.30, 16.34, 16.29, 17.19, 20.29, 21.89, 23.42, 24.84, 26.16]'/100;

CurrentForwardValue = (1 / 0.25) * (- 1 + ( 1+ 0.0756 * 0.5) / (1 + 0.0746 * 0.25));

objFun = @(X) MarketVolatilities - ...
blackvolbysabr(X(1), Beta1, X(2), X(3), Settle, ...
ExerciseDate, CurrentForwardValue, MarketStrikes);

options_SABR = optimoptions(@lsqnonlin, 'Display', 'iter', ...
'MaxFunctionEvaluations', 5000, 'MaxIterations', 10000);

%Minimizer:
X = lsqnonlin(objFun, [0.5 0 1], [0 -1 0], [Inf 1 Inf], options_SABR);

Alpha1 = X(1);
Rho1 = X(2);
Nu1 = X(3);

Sabr_vol = blackvolbysabr(Alpha1, Beta1, Rho1, Nu1, Settle, ExerciseDate,
CurrentForwardValue, MarketStrikes);

[Alpha1, Rho1, Nu1]

plot(MarketStrikes/T, MarketVolatilities)
hold on
plot(MarketStrikes/T, Sabr_vol) %shifted to see the curvature
legend("MarketVolatilities", "SabrVolatilieis")


### Python code

import pysabr
import numpy as np
from pysabr import Hagan2002LognormalSABR as LNsabr
from pysabr import hagan_2002_lognormal_sabr as hagan2002LN
import matplotlib.pyplot as plt

strikes = np.array([0.05, 0.055, 0.06, 0.0650, 0.07, 0.08, 0.0850, 0.09,
0.095, 0.10])
LogNormalVols = np.array([18.90, 17.30, 16.34, 16.29, 17.19, 20.29, 21.89, 23.42, 24.84, 26.16])
vol_n_ATM = 141.01 / 10000.00
forward_3m_6m = (1 / 0.25) * (- 1 + ( 1+ 0.0756 * 0.5) / (1 + 0.0746 * 0.25))
discount_0m_6m = 0.0753

beta = 0.5
calibration_LN = LNsabr(f = forward_3m_6m, shift = 0, t = 0.5, beta = beta).fit(strikes, LogNormalVols)
Matlab_LN = [0.0504, 0.6407, 0.8797] #MATLAB calibrated parameters! WORKS!
modelVols_LN = []
test_LN = []
for strike in strikes:
modelVols_LN.append(LNsabr(forward_3m_6m, 0, 0.5, vol_n_ATM, beta, calibration_LN[1], calibration_LN[2]).lognormal_vol(strike) * 100.00)
test_LN.append(hagan2002LN.lognormal_vol(strike, f = forward_3m_6m, t = 0.5, alpha = calibration_LN[0], beta = beta, rho = calibration_LN[1], volvol = calibration_LN[2]) * 100.00)

print(np.around(modelVols_LN,3))
print(np.around(test_LN,3))
#print(hagan2002LN.lognormal_vol(k = 0.0753, f = forward_3m_6m, t = 0.5, alpha = calibration_LN[0], beta = beta, rho = calibration_LN[1], volvol = calibration_LN[2]) * 100.00)

plt.plot(strikes, modelVols_LN, "r--", strikes, test_LN)
print("------------------------------------------------------------------------------------------")
print("------------------------------------------------------------------------------------------")
print("Matlab calibrated parameters: %s. \nPysabr calibrated parameters: %s" % (Matlab_LN, calibration_LN))

plt.plot(strikes, modelVols_LN, "r--", strikes, test_LN)

• This answer solved all my problems. Thanks a lot -- a well deserved bounty! Oct 13 at 6:04
• @Hasek Glad I could help :-)
– Pleb
Oct 13 at 6:45