In the previous article on calculating implied volatility for options we made use of interval bisection to numerically solve for the implied volatility. In this article we are going to modify our code to make use of the Newton-Raphson process, which is more optimal for this problem domain than interval bisection.

For the previous calculation we made use of a function template to carry out the interval bisection. The template function accepted a function object (functor) of type T, $g$, which itself accepted a volatility parameter ($\sigma$) to provide an option price. The main design issue we have to contend with is that to use Newton-Raphson, we require a second function to represent $g'$, the derivative of $g$ (the vega of the option). There are multiple approaches to deal with this issue:

• Two functors - We can create a separate function object to represent the derivative of the function stored in the first functor. However, this method is not scalable if we wish to create second (or higher order) derivatives, or partial derivatives for any other variables. A derivative is also fundamentally a component of the original function, rather than a separate entity. Thus we won't be using this approach.
• Derivative method - It is possible to create a derivative method for the function object $g$. However, we still suffer from the same scaling issue in that we would need to create derivative methods for all orders or variables. We could create a generalised derivative function, but this would be a substantial amount of work for very little initial gain.
• Pointer to member function - To avoid the ugliness of the previous method, we can make use of a pointer to a member function. This method allows us to specify which function is to be called at the point of compilation, using templates. It is the approach taken in Joshi^[1] and we will follow it here.

Before considering the C++ implementation, we will briefly discuss how the Newton-Raphson root-finding algorithm works.

Newton-Raphson Method

Newton-Raphson is a more efficient algorithm for finding roots provided that some assumptions are met. In particular, $g$ must possess an easily calculated derivative. If the derivative is analytically calculable, then this aids efficiency even further. $g$ must also be 'well behaved', that is it cannot vary too wildly, else the derivative approximation made use of in the method becomes more innacurate.

In the case of the Black-Scholes formula, the derivative we seek is the derivative of the option price with respect to the volatility, $\frac{\partial B}{\partial \sigma}$. This is known as the vega of the option model. For a call option an analytical function exists for the vega, so we are in luck.

The idea of Newton-Raphson is to use the analytic derivative to make a linear estimate of where the solution should occur, which is much more accurate than the mid-point approach taken by Interval Bisection. Thus the starting approximation to $g$, $g_0$, is given by (where $x_0$ is our initial guess):

This can be rearranged for $x$:

This then leads to a recurrence relation for more accurate values of $x$:

As with interval bisection, the algorithm is terminated when $|g(x_n) - y| < \epsilon$, for some pre-specified tolerance $\epsilon$.

Rewritten with our financial notation (where I have suppressed the other parameters for the Black-Scholes formula, $B$), the relation becomes:

We can use this relation in order to find the implied volatility via a terminating criterion.

Pointer to a Member Function

The idea of a pointer to a member function is to restrict the usage of a function pointer to methods of a particular class. The scope resolution operator :: combined with the pointer dereference operator * is used to carry this out. In order to define such a function pointer we use the following syntax (where g_prime represents $g'$, the derivative of $g$):

double (T::*g_prime)(double) const

To invoke the method we use the following syntax:

(root_func.*g_prime)(x)

Where root_func is an object of type T, containing two methods g and g_prime, representing the value of the function and the derivative value of the function respectively that will be tested for roots.

One of the major benefits of this approach is that the calls to these function pointers can be inlined as they are treated as template parameters and so are evaluated at compile-time.

Implementing Newton-Raphson

In order to make use of the new approach, we need to add the explicit formula for a call option vega to our bs_prices.h header file and modify the original BlackScholesCall object we created in the previous tutorial to make use of it. Here is the added function to calculate the option vega (see the original tutorial for the full listing):

..
..

// Calculate the European vanilla call vega 'Greek' based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double call_vega(const double S, const double K, const double r, const double sigma, const double T) {
  return S * sqrt(T) * norm_pdf(d_j(1, S, K, r, sigma, T));
}

..
..

Here is the new header listing for black_scholes.h:

#ifndef __BLACK_SCHOLES_H
#define __BLACK_SCHOLES_H

class BlackScholesCall {
private:
  double S;  // Underlying asset price
  double K;  // Strike price
  double r;  // Risk-free rate
  double T;  // Time to maturity

public:
  BlackScholesCall(double _S, double _K, 
                   double _r, double _T);

  // This is the modified section. operator()
  // has been replaced with option_price and
  // we have added option_vega (both const)
  double option_price(double sigma) const;
  double option_vega(double sigma) const;
};
#endif

The source listing has also been changed to include the implementation of option_vega, given in black_scholes.cpp:

#ifndef __BLACK_SCHOLES_CPP
#define __BLACK_SCHOLES_CPP

#include "black_scholes.h"
#include "bs_prices.h"

BlackScholesCall::BlackScholesCall(double _S, double _K, 
                                   double _r, double _T) :
  S(_S), K(_K), r(_r), T(_T) {}

// Renamed from operator() to option_price()
double BlackScholesCall::option_price(double sigma) const {
  return call_price(S, K, r, sigma, T);
}

// New method added, which calls call_vega 
// to obtain the actual price
double BlackScholesCall::option_vega(double sigma) const {
  return call_vega(S, K, r, sigma, T);
}

#endif

The next stage is to create the Newton-Raphson solver itself. The function template will accept an object of type T (the functor) and two pointers to member functions (methods) of T, g and g_prime. Here is the listing for newton_raphson.h:

#ifndef __NEWTON_RAPHSON_H
#define __NEWTON_RAPHSON_H

#include <cmath>

template<typename T, 
         double (T::*g)(double) const,
         double (T::*g_prime)(double) const>
double newton_raphson(double y_target,       // Target y value
                      double init,           // Initial x value
                      double epsilon,        // Tolerance
                      const T& root_func) {  // Function objec
  // Set the initial option prices and volatility
  double y = (root_func.*g)(init);  // Initial option prices
  double x = init;                  // Initial volatility

  // While y and y_target are not similar enough
  // Take the vega of the option and recalculate
  // a new call price based on the best linear
  // approximation at that particular vol value
  while (fabs(y-y_target) > epsilon) {
    double d_x = (root_func.*g_prime)(x);
    x += (y_target-y)/d_x;
    y = (root_func.*g)(x);
  }
  return x;
}

#endif

Now we can create the main() function to wrap all of our code together:

#ifndef __MAIN_CPP
#define __MAIN_CPP

#include "black_scholes.h"
#include "newton_raphson.h"
#include <iostream>

int main(int argc, char **argv) {
  // First we create the parameter list
  double S = 100.0;  // Underlying spot price
  double K = 100.0;  // Strike price
  double r = 0.05;   // Risk-free rate (5%)
  double T = 1.0;    // One year until expiry
  double C_M = 10.5; // Option market price

  // Create the Black-Scholes Call functor
  BlackScholesCall bsc(S, K, r, T);

  // Newton Raphson parameters
  double init = 0.3;  // Our guess impl. vol of 30%
  double epsilon = 0.001;

  // Calculate the implied volatility
  double sigma = newton_raphson<BlackScholesCall, 
                                &BlackScholesCall::option_price,
                                &BlackScholesCall::option_vega>
    (C_M, init, epsilon, bsc);

  // Output the values
  std::cout << "Implied Vol: " << sigma << std::endl;

  return 0;
}

#endif

The output of the code is given by:

Implied Vol: 0.201317

This matches the implied volatility given in the previous article article on interval bisection. This process can be refined even further, making use of Brent's Method, however this will be the subject of a later article.

References

• [1] - Joshi, M.S. C++ Design Patterns and Derivatives Pricing, 2nd Ed, Cambridge University Press, 2008.