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):

\begin{eqnarray} g_0(x) = g(x_0) + (x - x_0) g'(x_0) \end{eqnarray}This can be rearranged for $x$:

\begin{eqnarray} x = \frac{y - g(x_0)}{g'(x_0)} + x_0 \end{eqnarray}This then leads to a recurrence relation for more accurate values of $x$:

\begin{eqnarray} x_{n+1} = \frac{y - g(x_n)}{g'(x_n)} + x_n \end{eqnarray}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:

\begin{eqnarray} \sigma_{n+1} = \frac{C_M - B(\sigma_n)}{\frac{\partial B}{\partial \sigma} (\sigma_n)} + \sigma_n \end{eqnarray}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.