Correcting Vega in Newton–Raphson Implied Volatility Estimation
Summary
This exchange diagnoses why a Newton–Raphson routine for Black–Scholes implied volatility diverges. The original code divides the option-price error by gamma, but implied volatility is the quantity being solved for, so the derivative needed for the update is vega. Gamma measures sensitivity to the underlying price and is not the correct slope for this iteration. The accepted answer replaces gamma with vega, defined for the option model used, while retaining an iterative volatility update.
The example shows the original volatility estimate turning negative and then exploding after the incorrect update; the suggested correction addresses that derivative mismatch. The response does not provide convergence results or discuss safeguards. Newton–Raphson can still be sensitive to its starting value and to cases with very small vega, so the code shown should not be taken as a robust solver without additional bounds or fallback methods.
Key ideas
- Newton–Raphson for implied volatility must use vega, the option price derivative with respect to volatility.
- Gamma measures curvature with respect to the underlying price and cannot serve as the volatility update denominator.
- The example's incorrect derivative leads to negative volatility estimates and numerical divergence.
- The proposed correction does not address convergence safeguards such as volatility bounds or fallback solvers.
Tags
Full text
# estimate implied volatility using newton-raphson in python
# estimate implied volatility using newton-raphson in python
I am trying to calculate the implied volatility using newton-raphson in python, but the value diverges instead of converge. What is wrong with the code?
```
s = stock price
k = strike
t = time to maturity
rf = risk free interest
cp = +/-1 call/put
price = option price
def newtonRap(cp, price, s, k, t, rf):
v = sqrt(2*pi/t)*price/s
print "initial volatility: ",v
for i in range(1, 10):
d1 = (log(s/k)+(rf+0.5*pow(v,2))*t)/(v*sqrt(t))
d2 = d1 - v*sqrt(t)
gamma = norm.pdf(d1)/(s*v*sqrt(t))
price0 = cp*s*norm.cdf(cp*d1) - cp*k*exp(-rf*t)*norm.cdf(cp*d2)
v = v - (price0 - price)/gamma
print "price, gamma, volatility\n",(price0, gamma, v)
if abs(price0 - price) < 1e-10 :
break
return v
v = newtonRap(cp=1, price = 1.52, s=23.95, k=24, t=71.0/365, rf=0.05)
print v
```
output:
```
initial volatility: 0.36069926906
price, gamma, volatility
(1.6055072570611344, 0.10385864414094476, -0.46260492259786345)
price, gamma, volatility
(-1.8488599102758396, -0.080851497020229368, -42.129859511995726)
price, gamma, volatility
(-23.767706818545953, -1.6137689013848907e-22, -1.5669967860233743e+23)
price, gamma, volatility
(-23.767706818545953, -0.0, -inf)
RuntimeWarning: divide by zero encountered in double_scalars
v = v - (price0 - price)/gamma
RuntimeWarning: invalid value encountered in double_scalars
d1 = (log(s/k)+(rf+0.5*pow(v,2))*t)/(v*sqrt(t))
price, gamma, volatility
(nan, nan, nan)
price, gamma, volatility
(nan, nan, nan)
```
## Answer by user2686641 (score -1)
https://quant.stackexchange.com/a/11412
```
def newtonRap(cp, price, s, k, t, rf):
v = sqrt(2*pi/t)*price/s
print "initial volatility: ",v
for i in range(1, 100):
d1 = (log(s/k)+(rf+0.5*pow(v,2))*t)/(v*sqrt(t))
d2 = d1 - v*sqrt(t)
vega = s*norm.pdf(d1)*sqrt(t)
price0 = cp*s*norm.cdf(cp*d1) - cp*k*exp(-rf*t)*norm.cdf(cp*d2)
v = v - (price0 - price)/vega
print "price, vega, volatility\n",(price0, vega, v)
if abs(price0 - price) < 1e-25 :
break
return v
```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.