# Bachelier Normal Implied Vol Python Script Jaekel Not computing?

Encountering issue with Python Script for computing Nomral Implied Volatility from (Bacherlier). Using industry standard method, Jaekel-> https://jaeckel.000webhostapp.com/ImpliedNormalVolatility.pdf

Jaekel provides the numerical computation in the second part of the paper with machine accuracy, kind of incredible.

I have implemented Numerical method, it is not recovering the implied volatility of 20, please help. Code below.

# -*- coding: utf-8 -*-
"""
Created on Fri Jun  5 13:56:17 2020

@author: JS
"""

import scipy as sc
import numpy as np

#for Call priced at BachVol of 20
F=30
#not sure if Jaekel wants days to expirey or fraction->(days to expirey/365)
T=100/365
K=50
Sigma=20

#priced call following model equation Gives price of .112604141934456
V=Sigma*np.sqrt(T)*sc.stats.norm.pdf((F-K)/(Sigma*np.sqrt(T)))+(F-K)*sc.stats.norm.cdf((F-K)/(Sigma*np.sqrt(T)))
print(V)

#from equation (1.6) conditioned as stated for this example gives -.0055
PhiHatXstar= -V/20

#since -.0055 < -.00188203927  ; where PhiHatC=PhiHat(-9/4)= -.00188203927  apply (2.1)
#Assuming all the numbers in equation 2.2 for Eta are the same for all options, strikes and call values

#equation 2.1
g=1/(PhiHatXstar-.5)

#equation 2.2
Eta=(0.032114372355 - (g**2)*(0.016969777977 - (g**2)*(2.6207332461E-3-(9.6066952861E-5)*g**2)))/(1-(g**2)*(0.6635646938 - (g**2)*(0.14528712196 - 0.010472855461*g**2)))

#equation 2.3
XBar=((1/np.sqrt(2*np.pi))+Eta*g**2)

#equation 2.7 to solve XStar
PhiHatXBar=sc.stats.norm.cdf(XBar,0,1)+(sc.stats.norm.pdf(XBar,0,1)/XBar)

q=(PhiHatXBar-PhiHatXstar)/(sc.stats.norm.pdf(XBar,0,1))

XStar = XBar+ (3*q*(XBar**2)*(2-q*XBar*(2 + XBar**2)))/(6 + q*XBar*(-12 + XBar*(6*q + XBar*(-6 + q*XBar*(3 + XBar**2)))))

#Finall implied volatilty of Bachelier model should be 20 in this case
SigmaBach=abs(K-F)/abs(XStar*np.sqrt(T))

#Incorrect Result gives 151.0624.... not 20 as expected
print(SigmaBach)

#going the other way when -.0055 > -.00188203927  ; where PhiHatC=PhiHat(-9/4)= -.00188203927 Apply 2.4 & 2.5 instead

h=np.sqrt(-np.log(-PhiHatXstar))
XBar2=(9.4883409779-h*(9.6320903635-h*(0.58556997323 + 2.1464093351*h)))/(1-h*(0.65174820867 + h*(1.5120247828 + 6.6437847132E-5*h)))

PhiHatXBar2=sc.stats.norm.cdf(XBar2,0,1)+(sc.stats.norm.pdf(XBar2,0,1)/XBar2)

q2=(PhiHatXBar2-PhiHatXstar)/(sc.stats.norm.pdf(XBar2,0,1))

XStar2=XBar2+ (3*q*(XBar2**2)*(2-q*XBar2*(2 + XBar2**2)))/(6 + q*XBar2*(-12 + XBar2*(6*q + XBar2*(-6 + q*XBar2*(3 + XBar2**2)))))

SigmaBach2=abs(K-F)/abs(XStar2*np.sqrt(T))

#Also Incorrect Result 41.0308.... not 20  as expected
print(SigmaBach2)


runfile('', wdir='')