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

[Hasek]

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
from pysabr import Hagan2002LognormalSABR as LNsabr
import numpy as np

from pysabr import hagan_2002_lognormal_sabr as hagan2002
from pysabr import hagan_2002_normal_sabr as hagan2002normal
from pysabr import black

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 = LNsabr(f = forward_3m_6m, shift = 0, t = 0.5, beta = 0.5).fit(testStrikes, testVols)
smile = []
smile2 = []
for strike in testStrikes:
    smile.append(LNsabr(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)
    smile2.append(hagan2002.lognormal_vol(strike, forward_3m_6m, 0.5, calibration[0], 0.5, calibration[1], calibration[2]) * 100.00)

print(smile)
print(smile2)
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)

The output is

[23.525799837769064, 16.226199750352116, 20.186954010467034, 26.176954772322432]
[23.526566860448394, 16.227191003316406, 20.18810462154904, 26.178058282618206]
18.369296047068513

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, forward_3m_6m, 0.5, calibration_LN[0], 0.5, calibration_LN[1], calibration_LN[2]) * 100.00)

print(modelVols_LN)
print(test_LN)
print(hagan2002LN.lognormal_vol(0.0753, forward_3m_6m, 0.5, calibration_LN[0], 0.5, calibration_LN[1], calibration_LN[2]) * 100.00)

The output is

[20.14451195328051, 18.542578468604265, 17.499371066739126, 17.195175016467612, 17.677189818996556, 20.011518053296353, 21.365538837348016, 22.691679574639426, 23.956165099043584, 25.148945302368077]
[68.43385412708137, 67.98546916428838, 67.83926373469444, 67.91509162220223, 68.15026028834066, 68.91958765670267, 69.39185677119379, 69.89474561228646, 70.41454234659224, 70.94145255911624]
68.52050168384525

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.