# 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);
}

• A few more details, please: what original data do you have and what function are you trying to set to 0?
– user59
Oct 29, 2014 at 16:16
• Just added the C# code.
– opt
Oct 29, 2014 at 21:22
• You can have a look at this answer, which points out to a couple of alternatives.
– SRKX
Oct 30, 2014 at 1:29

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

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.

• Why out-of-the-money is ideal? Oct 31, 2014 at 13:08
• the intrinsic value of the in the money tends to swamp the numerics. Oct 31, 2014 at 21:07
• He actually wrote a more recent paper ("Let's be rational") in which he "show how Black's volatility can be implied from option prices with as little as two iterations to maximum attainable precision on standard (64 bit floating point) hardware for all possible inputs."
– Bram
Nov 3, 2014 at 18:57
• @MarkJoshi One must note that as the OP is apparently coding in c#, she will encounter issues related to stackoverflow.com/questions/3580389/… if she only adapts Jäckel's c++ code. (c# not having DBL_EPSILON is already a bad sign though.) Handling machine precision in c# to ensure max 2 iterations for any inputs as Jäckel did in c++ is non trivial. Just mutating Jäckels code to c# gives really unacceptable errors, for instance. Apr 2, 2019 at 7:28
• the website is not findable anymore :( Feb 10, 2020 at 16:39

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.

• Almost all programming languages and numerical systems should have some zero finder usually based on Brent. These will likely be faster and more stable than anything one might build yourself. Oct 29, 2014 at 16:52
• Thanks for the comments. Could you please address me towards some solution implemented in C#/C++ based on the Dekker and Brent algorithm?
– opt
Oct 29, 2014 at 20:55
• Is the bisection method really that slow? If you bracket between volatility 0 and 1, as few as 20 steps will get you within 10^-6 of the true volatility.
– user59
Oct 29, 2014 at 23:16
• Speed is always relative :) Every milli second counts in building pricing or enterprise solutions. If you are building number of Vol surfaces from different inputs then probably you want to reduce total time. Oct 30, 2014 at 3:13
• @Opt: I suggest you to separate your root search algorithms from other code. Root search algos are useful in other places as well so build them separately and plug them where ever required Oct 30, 2014 at 3:40

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