Robust Numerical Methods for Black–Scholes Implied Volatility
Article Quant Q&A · Author: opt
Summary
This discussion addresses failures of Newton–Raphson when solving for implied volatility from a Black–Scholes option price. Newton’s method can be fast near a root, but may fail to converge when its starting value or iteration steps are unsuitable. The code’s initial guess and update are presented as a practical example of an approach that can fail, though the responses do not diagnose that specific implementation in detail.
The answers recommend first expressing the option as out of the money, using put–call parity when needed, and point to a specialized method for robust Black implied-volatility calculation. They also contrast bracketing approaches, which preserve a root-containing interval but can be slow, with faster methods that may have convergence problems. Brent or related safeguarded methods combine bracketing with interpolation to improve reliability. A supplied Dekker implementation is explicitly offered with a warning to check it; the discussion gives no comparative benchmarks or comprehensive failure conditions.
Key ideas
- Newton–Raphson can converge quickly but is sensitive to the starting point and iteration behavior.
- Put–call parity can transform an in-the-money option into an out-of-the-money case for solving volatility.
- Bisection and related bracketing methods are reliable when the initial interval contains a root, though they may be slow.
- Brent and Dekker methods combine bracketing with interpolation for a balance of speed and robustness.
- The example root-search code is not accompanied by benchmarks and should be checked before use.
Tags
Full text
# What is an efficient method to find implied volatility?
# What is an efficient method to find implied volatility?
I have a code that finds the implied volatility using the Newton-Raphson method.
I set the number of trial to 1000 but sometimes it fails to converge and doesn't find the result.
Is there a better method to find the result? Are there any technical conditions in which this numerical method is expected to fail to converge to the solution?
Here is the C# code:
```
public double findIV(double S, double K, double r, double time, string type, double optionPrice)
{
int trial= 1000;
double ACCURACY = 1.0e-5;
double t_sqrt = Math.Sqrt(time);
double sigma = (optionPrice / S) / (0.398 * t_sqrt); // find initial value
for (int i = 0; i < trial; i++)
{
Option myCurrentOpt = new Option(type, S, K, time, r, 0, sigma); // create an Option object
double price = myCurrentOpt.BlackScholes();
double diff = optionPrice - price;
if (Math.Abs(diff) < ACCURACY)
return sigma;
double d1 = (Math.Log(S / K) + r * time) / (sigma * t_sqrt) + 0.5 * sigma * t_sqrt;
double vega = S * t_sqrt * ND(d1);
sigma = sigma + diff / vega;
}
throw new Exception("Failed to converge.");
}
public double ND(double X)
{
return (1.0 / Math.Sqrt(2.0 * Math.PI)) * Math.Exp(-0.5 * X * X);
}
```
## Answer by Mark Joshi (score 7, accepted)
https://quant.stackexchange.com/a/15232
Peter Jaeckel wrote a paper just on how to solve this problem:
By Implication (July 2006; Wilmott, pages 60-66, November 2006). Probably the most complicated trivial issue in financial mathematics: how to compute Black's implied volatility robustly, simply, efficiently, and fast
downloadable from jaeckel.org
In my experience the most important thing is to make sure that you are working with an option out of the money. If the option is in the money use put-call parity to transform to the other case.
## Answer by sigirisetti (score 3)
https://quant.stackexchange.com/a/15203
Bracketing methods such as Bisection and Regula Falsi are always known to converge but they are very slow.
Newton Raphson and secant methods are fast (quadratic convergence) but has convergence problems. Google for Newton Raphson convergence pitfalls. Classical ones such as"Trapped in local minima", "Diverge instead of converge" etc
Algorithms such as Dekker and Brent combine bracketing and quadratic convergence features and they are sure to converge and relatively faster as well.
## Answer by sigirisetti (score 3)
https://quant.stackexchange.com/a/15214
Below is the root search algorithm code I wrote in college. This is written in octave. It's simple to understand and re-write in C++.
Develop numerical methods algos as a separate module and integrate with your pricing and other code
I want to WARN you to re-check for bugs. It always converges for my objective functions
First function is Dekker method similar to Brent. Brent is more better. Second function is bracketing method. Either one can be used for root search
```
##
## Dekker root search.
##
## Eg call: dekker(@expMinusOne, -2, 3)
##
function c = dekker(f, lowerBound, upperBound)
MAX_ITER = 50;
iterations = 0;
a = lowerBound;
b = upperBound;
assert(validateBracket(f, a, b), "Given bounds doesn't bracket the root");
assert(!hasMultipleRoots(f, a, b), "Given bounds has mutiple roots");
fa = f(a);
fb = f(b);
#printf("fa = %f, fb = %f\n", fa, fb);
if(abs(fa) < abs(fb))
c = a;
d = b;
else
c = b;
d = a;
endif
linOps = 0;
do
s = linearInterpolation(f, c, d);
m = bisectionPoint(c, d);
a = c;
b = d;
# LI/LE and BI points are on same side of root
if(f(s) * f(m) >= 0)
# New root and lower bound are on same side of the root
# We can adjust one side of the bracket
if(f(c) * f(s) > 0)
if(abs(c - m) < abs(c - s) && linOps < 4)
c = s;
linOps++;
else
c = m;
linOps = 0;
endif
# New root and lower bound are on different sides of the root
# We can adjust both sides of the bracket
else
# Reduce the bracket.
if(abs(c - s) < abs(c - m))
c = s;
else
c = m;
endif
d = a;
endif
# LI/LE and BI points are on different sides of root.
# We can adjust both sides of the bracket
else
if(f(c) * f(s) > 0)
c = s;
d = m;
linOps++;
else
c = m;
d = s;
linOps = 0;
endif
endif
#printf("c = %f, d = %f\n", c, d);
# Check for Convergence
if(c == d)
#printf("DEKKER ROOT = %f in %d iterations.\n", c, iterations);
break;
endif
iterations++;
until (iterations == MAX_ITER)
if(iterations == MAX_ITER)
#printf("DEKKER ROOT = %f for maximum iterations %d.\n", c, MAX_ITER);
endif
endfunction
function s = linearInterpolation(f, a, b)
s = a - (b - a) * f(a)/(f(b) - f(a));
endfunction
function m = bisectionPoint(a, b)
m = (a + b)/2;
endfunction
function r = hasMultipleRoots(f, a, b)
if(f(a) == f(b))
r = true;
else
r = false;
endif
endfunction
function r = validateBracket(f, a, b)
if(f(a) * f(b) < 0)
r = true;
else
r = false;
endif
endfunction
##
## regulaFalsi( f, a, b, yAcc )
##
## rootsearching by regulaFalsi method
##
## f is a real numeric function s.t. f( a ) * f( b ) < 0.0
## xAcc convergence threshold
## nIter max number of iterations
function [c, x] = regulaFalsi( f, a, b, yAcc )
fb = f( b );
fa = f( a );
if ( fa * fb >= 0.0 ),
error(" f( a ) * f( b ) >= 0.0 " );
endif
# init loop variables
c = a;
iter = 1;
x = [];
do
cOld = c; # save previous value of c. Guard for infinite loop
c = a - fa * ( b - a ) / ( fb - fa ); # new point
fc = f( c );
if ( fc * fb < 0.0 ) # update interval
a = c;
fa = fc;
else
b = c;
fb = fc;
endif
x(iter,:)=[iter, c]; iter++; # unnecessary, just for display
until ( abs(fc) < yAcc || cOld == c ) # Exit criteria is on yAcc.
endfunction
```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.