7
$\begingroup$

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);
    }
Improve this question
$\endgroup$
3
  • $\begingroup$ A few more details, please: what original data do you have and what function are you trying to set to 0? $\endgroup$
    – user59
    Oct 29, 2014 at 16:16
  • $\begingroup$ Just added the C# code. $\endgroup$
    – opt
    Oct 29, 2014 at 21:22
  • $\begingroup$ You can have a look at this answer, which points out to a couple of alternatives. $\endgroup$
    – SRKX
    Oct 30, 2014 at 1:29

3 Answers 3

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Why out-of-the-money is ideal? $\endgroup$
    – SmallChess
    Oct 31, 2014 at 13:08
  • 1
    $\begingroup$ the intrinsic value of the in the money tends to swamp the numerics. $\endgroup$
    – Mark Joshi
    Oct 31, 2014 at 21:07
  • 1
    $\begingroup$ 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." $\endgroup$
    – Bram
    Nov 3, 2014 at 18:57
  • 1
    $\begingroup$ @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. $\endgroup$
    – Olórin
    Apr 2, 2019 at 7:28
  • $\begingroup$ the website is not findable anymore :( $\endgroup$
    – Marine Galantin
    Feb 10, 2020 at 16:39
3
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – rhaskett
    Oct 29, 2014 at 16:52
  • $\begingroup$ Thanks for the comments. Could you please address me towards some solution implemented in C#/C++ based on the Dekker and Brent algorithm? $\endgroup$
    – opt
    Oct 29, 2014 at 20:55
  • $\begingroup$ 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. $\endgroup$
    – user59
    Oct 29, 2014 at 23:16
  • $\begingroup$ 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. $\endgroup$
    – sigirisetti
    Oct 30, 2014 at 3:13
  • $\begingroup$ @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 $\endgroup$
    – sigirisetti
    Oct 30, 2014 at 3:40
3
$\begingroup$

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
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.