Bachelier Normal Implied Vol Python Calculation (Help) Jekel

Question

Writing a python script to calc Implied Normal Vol ; in line with Jekel article (Industry Standard).

https://jaeckel.000webhostapp.com/ImpliedNormalVolatility.pdf

They say they are using a Generalized Incomplete Gamma Function Inverse.

For a call:

F(x)=v/(K - F) -> find x that makes this true

Where F is Inverse Incomplete Gamma Function

And x = (K - F)/(T*sqrt(T) ; v is the value of a call

for that x, IV is =(K-F)/x*sqrt(T)

Example I am working with:

F=40

X=38

T=100/365

v=5.25

Vol= 20%

Using the equations I should be able to backout Vol of 20%

Scipy has upper and lower Incomplete Gamma Function Inverse in their special functions.

Lower: scipy.special.gammaincinv(a, y) : {a must be positive param}

Upper: scipy.special.gammainccinv(a, y) : {a must be positive param}

Implementation:

SIG= sympy.symbols('SIG')
F=40
T=100/365
K=38

def Objective(sig):
    SIG=sig
    return(special.gammaincinv(.5,((F-K)**2)/(2*T*SIG**2))+special.gammainccinv(.5,((F-K)**2)/(2*T*SIG**2))+5.25/(K-F))

x=optimize.brentq(Objective, -20.00,20.00, args=(), xtol=1.48e-8, rtol=1.48e-8, maxiter=1000, full_output=True)


IV=(K-F)/x*T**.5

Print(IV)

I know I am wrong, but Where am I going wrong / how do I fix it and use what I read in the article ?

Answer 1

Did you also post this on the Quantitative Finance Stack Exchange? You may get a better response there.

This is not my field, but it looks like your main problem is that brentq requires the passed Objective function to return values with opposite signs when passed the -20 and 20 arguments. However, this will not end up happening because according to the scipy docs, gammaincinv and gammainccinv always return a value between 0 and infinity.

I'm not sure how to fix this, unfortunately. Did you try implementing the analytic solution (rather than iterative root finding) in the second part of the paper?