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

$\endgroup$
  • $\begingroup$ Why out-of-the-money is ideal? $\endgroup$ – SmallChess Oct 31 '14 at 13:08
  • 1
    $\begingroup$ the intrinsic value of the in the money tends to swamp the numerics. $\endgroup$ – Mark Joshi Oct 31 '14 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 '14 at 18:57
  • $\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$ – Olorin Apr 2 at 7:28
2
$\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.

$\endgroup$
  • $\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 '14 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 '14 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$ – barrycarter Oct 29 '14 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 '14 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 '14 at 3:40
2
$\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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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