Bartlett's delta gives wrong signs for calls and puts

Question

There is a paper by Bruce Bartlett introducing a modified delta for SABR model which accounts for the correlation between forward and volatility processes. The main result of the paper is that if $dF$ is a change in a forward rate $F$ then the average change in a SABR volatility parameter $\alpha$ is $\delta\alpha=\frac{\rho\nu}{F^\beta}dF$.

I took the following code computing Bartlett's delta from the SABR and SABR LIBOR Market Models in Practice book:

from scipy.stats import norm
import math

def haganLogNormalApprox (y, expiry , F_0 , alpha_0 , beta , nu , rho ):
    '''
    Function which returns the Black implied volatility ,
    computed using the Hagan et al. lognormal
    approximation .
    @var y: option strike
    @var expiry: option expiry (in years)
    @var F_0: forward interest rate
    @var alpha_0: SABR Alpha at t=0
    @var beta : SABR Beta
    @var rho: SABR Rho
    @var nu: SABR Nu
    '''
    one_beta = 1.0 - beta
    one_betasqr = one_beta * one_beta
    if F_0 != y:
        fK = F_0 * y
        fK_beta = math .pow(fK , one_beta / 2.0)
        log_fK = math .log(F_0 / y)
        z = nu / alpha_0 * fK_beta * log_fK
        x = math .log (( math .sqrt (1.0 - 2.0 * rho *
        z + z * z) + z - rho) / (1 - rho))
        sigma_l = (alpha_0 / fK_beta / (1.0 + one_betasqr /
        24.0 * log_fK * log_fK +
        math .pow( one_beta * log_fK , 4) / 1920.0) *
        (z / x))
        sigma_exp = ( one_betasqr / 24.0 * alpha_0 * alpha_0 /
        fK_beta / fK_beta + 0.25 * rho * beta *
        nu * alpha_0 / fK_beta +
        (2.0 - 3.0 * rho * rho) / 24.0 * nu * nu)
        sigma = sigma_l * ( 1.0 + sigma_exp * expiry)
    else:
        f_beta = math .pow(F_0 , one_beta)
        f_two_beta = math .pow(F_0 , (2.0 - 2.0 * beta ))
        sigma = (( alpha_0 / f_beta) * (1.0 +
        (( one_betasqr / 24.0) *
        ( alpha_0 * alpha_0 / f_two_beta ) +
        (0.25 * rho * beta * nu * alpha_0 / f_beta) +
        (2.0 - 3.0 * rho * rho) /
        24.0 * nu * nu) * expiry))
    return sigma

def dPlusBlack(F_0 , y, expiry , vol):
    '''
    Compute the d+ term appearing in the Black formula.
    @var F_0: forward rate at time 0
    @var y: option strike
    @var expiry: option expiry (in years)
    @var vol: Black implied volatility
    '''
    d_plus = ((math.log(F_0 / y) + 0.5 * vol * vol * expiry)
    / vol / math.sqrt(expiry ))
    return d_plus

def dMinusBlack(F_0 , y, expiry , vol):
    '''
    Compute the d- term appearing in the Black formula.
    @var F_0: forward rate at time 0
    @var y: option strike
    @var expiry: option expiry (in years)
    @var vol: Black implied volatility
    '''
    d_minus = (dPlusBlack(F_0 = F_0 , y = y, expiry = expiry ,
    vol = vol ) - vol * math.sqrt(expiry ))
    return d_minus

def black(F_0 , y, expiry , vol , isCall ):
    '''
    Compute the Black formula.
    @var F_0: forward rate at time 0
    @var y: option strike
    @var expiry: option expiry (in years)
    @var vol: Black implied volatility
    @var isCall: True or False
    '''
    option_value = 0
    if expiry * vol == 0.0:
        if isCall:
            option_value = max(F_0 - y, 0.0)
        else:
            option_value = max(y - F_0 , 0.0)
    else:
        d1 = dPlusBlack(F_0 = F_0 , y = y, expiry = expiry ,
        vol = vol)
        d2 = dMinusBlack(F_0 = F_0 , y = y, expiry = expiry ,
        vol = vol)
        if isCall:
            option_value = (F_0 * norm.cdf(d1) - y *
            norm.cdf(d2))
        else:
            option_value = (y * norm.cdf(-d2) - F_0 *
            norm.cdf(-d1))
    return option_value

def computeFirstDerivative (v_u_plus_du , v_u_minus_du , du):
    '''
    Compute the first derivatve of a function using
    central difference
    @var v_u_plus_du: is the value of the function
    computed for a positive bump amount du
    @var v_u_minus_du : is the value of the function
    computed for a negative bump amount du
    @var du: bump amount
    '''
    first_derivative = (v_u_plus_du - v_u_minus_du ) / (2.0 * du)
    return first_derivative

def computeSABRDelta (y, expiry , F_0 , alpha_0 , beta , rho , nu , isCall):
    '''
    Compute the SABR delta.
    @var y: option strike
    @var expiry: option expiry (in years)
    @var F_0: forward interest rate
    @var alpha_0: SABR Alpha at t=0
    @var beta : SABR Beta
    @var rho: SABR Rho
    @var nu: SABR Nu
    @var isCall: True or False
    '''
    small_figure = 0.0001
    F_0_plus_h = F_0 + small_figure
    avg_alpha = (alpha_0 + (rho * nu /
    math .pow(F_0 , beta )) * small_figure )
    vol = haganLogNormalApprox (y, expiry , F_0_plus_h , avg_alpha ,
    beta , nu , rho)
    px_f_plus_h = black(F_0_plus_h , y, expiry , vol , isCall)
    F_0_minus_h = F_0 - small_figure
    avg_alpha = (alpha_0 + (rho * nu /
    math .pow(F_0 , beta )) * (-small_figure ))
    vol = haganLogNormalApprox (y, expiry , F_0_minus_h ,
    avg_alpha , beta ,
    nu , rho)
    px_f_minus_h = black(F_0_minus_h , y, expiry , vol , isCall)
    sabr_delta = computeFirstDerivative (px_f_plus_h ,px_f_minus_h ,
    small_figure )
    return sabr_delta

The code seems alright however I encountered a problem with wrong signs of deltas for several caplets (call option on a forward rate) and floorlets (put option on a forward rate) while working with SABR model calibrated to real surface. One would expect the delta of a call to be positive and the delta of a put to be negative which is violated in the following case

BartlettDeltaPut = computeSABRDelta(y=0.06, expiry=1.50, F_0=0.0962688131761622, 
                                 alpha_0=0.0895853076638471, beta=0.5, rho=0.235477576202461, nu=1.99479846430177, 
                                 isCall=False)
BartlettDeltaCall = computeSABRDelta(y=0.10, expiry=0.25, F_0=0.07942844548137806, 
                                 alpha_0=0.127693338654331, beta=0.5, rho=-0.473149790316068, nu=2.46284420168144, 
                                 isCall=True)

resulting in

0.21186868757223573
-0.0012938212806158644

In a contrary, the plain vanilla Black delta given by

import numpy as np

def Delta(k, f, t, v, isCall=True):

    d1 = (np.log(f/k) + v**2 * t/2) / (v * t**0.5)
    if isCall:
        delta = norm.cdf(d1)
    else:
        delta = norm.cdf(d1) - 1

    return delta

vol1 = haganLogNormalApprox(y=0.06, expiry=1.50, F_0=0.0962688131761622, 
                           alpha_0=0.0895853076638471, beta=0.5, nu=1.99479846430177, rho=0.235477576202461)
vol2 = haganLogNormalApprox(y=0.10, expiry=0.25, F_0=0.07942844548137806, 
                           alpha_0=0.127693338654331, beta=0.5, nu=2.46284420168144, rho=-0.473149790316068)
BlackDeltaPut = Delta(k=0.06, f=0.0962688131761622, t=1.50, v=vol1, isCall=False)
BlackDeltaCall = Delta(k=0.10, f=0.07942844548137806, t=0.25, v=vol2, isCall=True)

coupled with volatility values computed by Hagan et al. approximations from the code above would work just as expected producing negative delta for put and positive delta for call options:

-0.16385166669719764
0.1753400660949036

Why Bartlett's delta values don't make sense in this case? I looked through the code carefully and to the best of my knowledge it doesn't have any errors or typos.

Answer 1

Bartlett's delta as computed in your code is a simple finite difference (FD), also called bump and reprice, of the Black values. I do not think there is anything wrong here, besides the fact that you are not comparing like for like. Long story short, Bartlett's delta more accurately describes the change in the value of the option because it "incorporates" a Vega (or vanna /Ddelta Dvol) component. For certain conditions, this can result in different signs relative to Black deltas (which ignore the shape of the vol surface).

The simple Black delta is a "holding everything else constant" delta, meaning IVOL is fixed and stable. What you compute with the Bartlett delta is a FD value with different IVs for both shifts.

Using your haganLogNormalApprox function results in the following vols (I renamed some variables due to personal preference):

F_0=0.07942844548137806 
small_figure = 0.0001
y=0.10
expiry=0.25
alpha_0=0.127693338654331
beta=0.5
rho=-0.473149790316068
nu=2.46284420168144
F_pos = F_0 + small_figure
F_neg = F_0 - small_figure
avg_alpha_pos = (alpha_0 + (rho * nu / math.pow(F_0 , beta )) * small_figure )
avg_alpha_neg = (alpha_0 + (rho * nu / math.pow(F_0 , beta )) * (-small_figure ))
vol_0 = vol_0 = haganLogNormalApprox(y=0.10, expiry=0.25, F_0=F_0, 
                       alpha_0=0.127693338654331, beta=0.5, nu=2.46284420168144, rho=-0.473149790316068)
vol_neg = haganLogNormalApprox (y, expiry , F_neg ,avg_alpha_neg , beta ,nu , rho)
vol_pos = haganLogNormalApprox (y, expiry , F_pos ,avg_alpha_pos , beta ,nu , rho)

print(f' Vol Black = {vol_0}')
print(f' Vol Positive = {vol_pos}')
print(f' Vol Negative = {vol_neg}')
vol_avg = (vol_pos + vol_neg)/2
print(f' Vol Average (Centred) = {vol_avg}')

Result:

 Vol Black = 0.44137660374291093
 Vol Positive = 0.4396541985333539
 Vol Negative = 0.4431007329951456
 Vol Average (Centred) = 0.44137746576424974

I also rewrote Black because I find it easier to work with less convoluted functions:

import numpy as np
from scipy.stats import norm

def B(F,K,t, sigma, cp_flag):
    d1 = ((np.log(F/K)  + 0.5 * sigma **2 * t) / (sigma * np.sqrt (t)))
    d2 = d1 - sigma * np.sqrt(t) 
    opt = cp_flag*(F*norm.cdf(cp_flag*d1) - K*norm.cdf(cp_flag*d2))
    delta = cp_flag*norm.cdf(cp_flag*d1)
    return opt, delta

Now, let's plug in the numbers:

up = B(F_pos,0.10,0.25, vol_pos, 1)
down = B(F_neg,0.10,0.25, vol_neg , 1)
centred = B(F_0,0.10,0.25, vol_0 , 1)

print(f'Call Val Up = {up[0]}, Delta Up = {up[1]}')
print(f'Call Val Down = {down[0]}, Delta Down = {down[1]}')
print(f'Call Val Centred = {centred[0]}, Delta Centred = {centred[1]}')
print(f'Delta FD = {(up[0]-down[0])/(2.0*small_figure)}')
print(f'Delta FD centred = {(B(F_pos,0.10,0.25, vol_0, 1)[0]- B(F_neg,0.10,0.25, vol_0, 1)[0])/(2.0*small_figure)}')

to obtain

Call Val Up = 0.0015010173779386893, Delta Up = 0.1756511129543487
Call Val Down = 0.0015012761421948125, Delta Down = 0.17503196098493057
Call Val Centred = 0.0015011434512114622, Delta Centred = 0.1753400660949036
Delta FD = -0.0012938212806158644
Delta FD centred = 0.1753410636751336

As you can see, using a constant vol (the average of the up and down vol specifically) results in pretty much the same delta as in Black. What you plug into FD is a different IV for up and down, which is in line with the SABR model. I am using the formula from this answer, which are the rquations shown on P.91 of the article Managing Smile risk from Hagan et. al in the Willmott magazine. Re-written in Python, the code looks like this:

def vol(β,α, ρ, ν, t_ex, f, K):
        A = α /(((f*K)**((1-β)/2))*(1+((1-β)**2)/24*np.emath.logn(2,(f/K))+ ((1-β)**4)/1920*np.emath.logn(4,(f/K))))
        B = 1+(((1-β)**2)/24*(α**2/(f*K)**(1-β))+(1/4)*α*β*ρ*ν/((f*K)**((1-β)/2))+(2-3*ρ**2)/24*ν**2)*t_ex
        z = ν/α*(f*K)**((1-β)/2)*np.log(f/K)
        χ_z = np.log((np.sqrt(1-2*ρ*z+z**2)+z-ρ)/(1-ρ))
        atm = α/(f**(1-β))*(1+(((1-β)**2)/24*(α**2/(f*K)**(1-β))+(1/4)*α*β*ρ*ν/((f*K)**((1-β)/2))+(2-3*ρ**2)/24*ν**2)*t_ex)
        cond = f==K
        return atm if cond else A*z/χ_z*B

If we plot this as a function of F, it is clear that it makes sense for IV to be above Black IV if F is shocked down, and below Black IV is F is shocked up, even if you would not adjust $α$ (which makes the effect a tad more pronounced):

β, α, ρ, ν, t_ex, f = 0.5, 0.127693338654331, -0.473149790316068, 2.46284420168144, 0.25, 0.07942844548137806
K=0.1

plt.plot( np.arange(0.05,0.14,0.003), [vol(β, alpha_0, ρ, ν, t_ex, F, K) for F in np.arange(0.05,0.14,0.003)], label = 'SABR Vols')
font = {'family':'serif','color':'blue','size':20}
plt.xlabel("Forward",fontdict = font)
plt.ylabel("IVOL",fontdict = font)
plt.title("SABR Model",fontdict = font)
plt.vlines(F_0, 0, vol(β, alpha_0, ρ, ν, t_ex, F_0, K), colors ="r", label = 'Current Forward')
plt.legend(loc = 'lower right')
plt.show()

The red line corresponds to the IV at the current forward. Despite being a different implementation, the resulting IV is almost identical to your implemention and also results in a similar negative delta when using a FD method with adjusted IV.

After all, the entire reason for Bartlett (2006) providing a refined delta under the SABR model is to account for the effect of vol. If you move along the vol surface when the underlying changes, you will have different market values, not only due to changes in the underlying, but also due to changes in IV.

However, it was shown that for a portfolio that is both delta and vega hedged, the original SABR Greeks given by Hagan et al. (2002) provide essentially the same result as Bartlett’s new SABR Greeks. See for example Hedging under SABR model by Bruce Bartlett in Wilmott Feb 2006. This claim is also easy to verify if we add vega to the Black formula:

def B(F,K,t, sigma, cp_flag):
    d1 = ((np.log(F/K)  + 0.5 * sigma **2 * t) / (sigma * np.sqrt (t)))
    d2 = d1 - sigma * np.sqrt(t) 
    opt = cp_flag*(F*norm.cdf(cp_flag*d1) - K*norm.cdf(cp_flag*d2))
    delta = cp_flag*norm.cdf(cp_flag*d1)
    vega = F*norm.pdf(d1)*np.sqrt(t)
    return opt, delta , vega

If you compute the new option value for a $F_{pos}$ value of the underlying, you get (using the new vol according to the vol surface):

res = B(F_pos,0.10,0.25, vol_pos , 1)
print(f'Call New = {res[0]}')

Call New = 0.0015010173779386893

Using a Taylor approximation with Black delta would yield:

print(f'Value according to Black Delta = {centred[0] + centred[1]*(F_pos-F_0)}')

Value according to Black Delta = 0.001518677457820953

However, this is too high. Using Bartlett's delta (which is in this case negative) gives

print(f"Value according to Bartlett's Delta = {centred[0] + (up[0]-down[0])/(2.0*small_figure)*(F_pos-F_0)}")

Value according to Bartlett's Delta = 0.0015010140690834006

which is almost identical to the value obtained by using Black delta and Vega together (and very close to the actual value of Black had you computed it outright without Greeks).

print(f'Value according to Black Delta and Vega = {centred[0] + centred[1]*(F_pos-F_0) + centred[2]*(vol_pos-vol_0)}')

Value according to Black Delta and Vega = 0.0015010228771048103