Bisection Calibration for Implied Volatility in Deep ITM and OTM Options
Summary
The document explains a numerical difficulty in recovering Black–Scholes implied volatility for in-the-money and out-of-the-money calls. It presents Newton–Raphson as a method for an at-the-money call, where price errors are adjusted using vega, and notes that this approach can fail when vega is very small for options far from the money. Two sample option cases illustrate the issue, one out of the money and one in the money.
As an alternative, the included answer proposes bisection over a transformed volatility range: repeatedly test the midpoint, compare its model price with the observed price, and keep the half-interval consistent with that comparison. This avoids relying on vega and gives a bracketed search. The document provides no convergence analysis, stopping rule beyond a maximum iteration count, or discussion of price bounds and invalid inputs. Its examples use calls under a Black–Scholes setup, so implementation details and robustness for other instruments are not established.
Key ideas
- Newton–Raphson updates implied volatility using the option price error divided by vega.
- Vega can become very small for far in-the-money or out-of-the-money options, destabilizing Newton updates.
- Bisection searches a bounded volatility interval by retaining the side consistent with the price comparison.
- The examples concern Black–Scholes calls and do not establish behavior for other option models.
- A practical implementation needs an appropriate bracket and a defined stopping criterion.
Tags
Full text
# Computing implied volatilities of ITM and OTM options
# Computing implied volatilities of ITM and OTM options
For an ATM call the implied volatility can be computed by using the Newton-Raphson method:
```
import numpy as np
from scipy.stats import norm
def bs_d1(S,K,T,sigma,R):
return (np.log(S/K) + (r+sigma**2/2)*T)/(sigma*np.sqrt(T))
def bs_call(S,K,T,sigma,r):
d1 = bs_d1(S,K,T,sigma,r)
return S*norm.cdf(d1) - np.exp(-r*T)*K*norm.cdf(d1-sigma*np.sqrt(T))
def bs_vega(S,K,T,sigma,r):
d1 = bs_d1(S,K,T,sigma,r)
return S*np.sqrt(T)*norm.pdf(d1)
def sigma_i(c,S,K,T,r,eps=1e-3):
# first approximation
sigma = np.sqrt(2*np.pi/T)*c/S
while True:
diff = bs_call(S,K,T,sigma,r) - c
vega = bs_vega(S,K,T,sigma,r)
if np.abs(diff) < eps:
break
sigma = sigma - diff/vega
return sigma
```
However the method fails when I use OTM and ITM calls since Vega becomes very close to zero. One example:
```
# OTM
S = 100
K = 140
T = 0.25
r = 0.03
c = 0.0015581689540368482 # bs_call(S,K,T,0.2,r)
sigma = sigma_i(c,S,K,T,r)
```
Another example:
```
# ITM
S = 100
K = 80
T = 1.5
r = 0.03
c = 42.36652921873271 # bs_call(S,K,T,0.7,r)
sigma = sigma_i(c,S,K,T,r)
```
How to compute the implied volatility in these cases?
## Answer by Valometrics.com (score 2)
https://quant.stackexchange.com/a/50841
you can use this algoithm to calibrate volatility (the use of vega deos not work for far OTM and ITM options):
```
define vol1=0.0;
define vol2=1.99;
define z;
define z2;
for (i=0;i<max_iterations;i++) {
z=0.5*(vol1+vol2);
z2=z/(1.0-z);
if (price>bs_price(z2)) {
vol1=z;
vol2=vol2;
} else {
vol1=vol1;
vol2=z;
}
}
z=0.5*(vol1+vol2);
return z/(1.0-z);
```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.