Debugging Jaekel’s Bachelier Implied Volatility Approximation
Summary
The post demonstrates a numerical implementation problem in recovering normal, or Bachelier, implied volatility from a call option price. The author prices a call under the Bachelier model with stated forward, strike, expiry, and volatility inputs, then applies an approximation method described in Jaekel’s paper. The initial implementation returns volatility estimates that differ substantially from the input, prompting a search for the coding or formula error.
A reply identifies a missing factor in the intermediate expression used to construct the approximation’s starting value. After that correction, the reported calculation recovers the original volatility to high numerical precision. The example is a useful illustration of checking an implied-volatility inversion against a price generated from known model inputs. It does not provide a general implementation, independent validation across option moneyness or expiry, or discussion of numerical edge cases, so the reported success should be understood as a single worked example.
Key ideas
- Bachelier option prices can be inverted to obtain normal implied volatility using a numerical approximation.
- A missing factor in an intermediate approximation step caused the example to return incorrect volatility estimates.
- Correcting that term made the inversion recover the volatility used to generate the sample price.
- Testing an inversion against a model price with known inputs is a useful implementation check.
- The post does not assess performance across a broad range of option inputs or numerical edge cases.
Tags
Full text
# Bachelier Normal Implied Vol Python Script Jaekel Not computing?
# 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)
```
## Answer by JBerstein (score 1)
https://quant.stackexchange.com/a/54684
was missing g* in XBar equation2/3; now computes Value correctly 20.0000000001 Jaekel is a machine.
```
runfile('', wdir='')
0.11260414193445623
20.00000000000001
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.