0
$\begingroup$

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)
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

was missing g* in XBar equation2/3; now computes Value correctly 20.0000000001 Jaekel is a machine.

runfile('', wdir='')
0.11260414193445623
20.00000000000001
Improve this answer
$\endgroup$
1
  • $\begingroup$ Interesting I added g* to equation 2.3 yet the IV from the 2nd part (not quite sure what you're doing there) is slightly off from 20: SigmaBach = 20.00000000000001 and SigmaBach2 = 19.980963318763294. Seems the error was just using q2 in the second set of equations and I get 20.000000000008367 instead (or change q2 to q): jaeckel.org/ImpliedNormalVolatility.pdf - I'm missing something in the paper as ATM gives me divide by 0 errors... I think that's covered in section 1 though (not explicitly). $\endgroup$
    – Matt
    Commented Jan 4, 2022 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.