In this article I want to discuss a practical application of the Black-Scholes model, design patterns and function objects in C++. In particular, we are going to consider the concept of Implied Volatility. In derivatives pricing, the implied volatility of an option is the value of the underlyings volatility (usually denoted by $\sigma$), which when input into an derivatives pricing model (such as with the Black-Scholes equation) provides a pricing value for the option which is equal to the currently observed market price of that option.

Motivation

The use of implied volatility is motivated by the fact that it is often more useful as a measure of relative value between options than the actual price of those options. This is because the option price is heavily dependent upon the price of the underlying asset. Options are often held in a delta neutral portfolio, which means that the portfolio is hedged for small moves in the price of the underlying asset. Given the need to continuously rebalance these portfolios, an option with a higher actual price can be cheaper in terms of its volatility if the underlying itself rises, as the underlying asset can be sold at a higher price.

Implied volatilities can in fact be considered a form of "prices", as their values have been determined by actual transactions in the marketplace, not on the basis of statistical estimates. Professional traders actually tend to quote values of options in terms of their "vol", rather than their actual market prices.

Our task for this article is to attempt to calculate implied volatilities using the Black-Scholes equation. If we let the market price of an option be given by $C_M$ and the Black-Scholes function by $B(S, K, r, T, \sigma)$, where $S$ is the underlying price, $K$ is the option strike, $r$ is the risk-free interest rate, $T$ is the time to maturity, then we are trying to find the value of $\sigma$ such that:

Root-Finding Algorithms

The first place to start is to consider whether an analytical inverse exists. Unfortunately, the Black-Scholes model does not lend itself to an analytical inverse and thus we are forced to carry out the computation using numerical methods. The common methods for achieving these are known as root finding algorithms. In this article we will look at Interval Bisection. In later articles we will consider Newton-Raphson and Brent's Method.

Interval Bisection Method

The first numerical algorithm considered is Interval Bisection. It makes use of a real analysis concept known as the intermediate value theorem. Since the Black-Scholes formula is continuous (and "well behaved", which for us means sufficiently smooth), we can make use of the theorem to help us find a volatility.

The theorem states, mathematically, that if $\exists m,n \in \mathbb{R}$, and $g(x):\mathbb{R}\rightarrow \mathbb{R}$, continuous, such that $g(m) < y$, $g(n) > y$, then $\exists p \in (m,n) \subset \mathbb{R}$ such that $g(p)=y$. For us, $g$ is the Black-Scholes function, $y$ is the current market price and $p$ is the volatility.

Heuristically, the theorem states that if a function is continuous and we know two (differing) values in the domain that generate an image interval, then there must exist a value in the domain between these two values that provides a mapped value lying within the image domain. This is exactly what we need in order to determine the implied volatility from the market prices.

Interval bisection is quite straightforward to understand. It is a "trial and error" algorithm. We pick the mid-point value, $c$, of an interval, and then either $g(c) = y$, $g(c) < y$ or $g(c) > y$. In the first instance the algorithm terminates. In the latter two cases, we subdivide the interval $(c,n)$ (respectively $(m,c)$) and find the corresponding midpoints. At each level of the recursion, the interval size is halved. In practice the algorithm is terminated using a tolerance value, $\epsilon$, once $|g(c) - y| < \epsilon$.

Interval bisection is not an efficient method for root-finding, but it is straightforward to implement, especially when we make use of functors and function templates.

Interval Bisection with Function Templates

So far we've made extensive use of object-orientation and, in particular, inheritance hierarchies as one of the main design patterns for solving quantitative finance problems. We're now going to use a very different set of tools, that of functors and function templates. We've discussed how functors work before and have used them to create PayOff classes. Now we are going to supplement their use by including function templates, which are a means of applying the template generic programming paradigm that we have previously applied to classes, to functions themselves.

The motivation for using functors and function templates is that we wish to create our interval bisection code in a reusable fashion. Ideally, the interval bisection code should be able to cope with a generic function, including those with their own set of (fixed) parameters. This is the situation we are in, because the Black-Scholes formula requires not only the volatility, but the strike, spot price, maturity time and interest rate.

Creating the Black-Scholes call option as a functor allows us to pass it to the interval bisection function with attached state, which in this case is the collection of extra arguments representing the option parameters. In fact, the interval bisection method doesn't even need to know about these parameters, as the only interface exposed to it is an option pricing operator() method, which accepts only a volatility.

Another reason to use function templatisation over inheritance is that we would be at the mercy of virtual functions. Virtual functions are relatively slow, as each time a function is called the vtable has to be queried, which is extra overhead. Further, virtual functions cannot be inlined, which is a major disadvantage. This leads to slower code. These problems also apply to function pointers as well, which is why we aren't making use of them either.

The approach we will follow is to make the interval_bisection function a template function, such that it can accept a generic function object, which in our case will be the Black-Scholes call pricing method. We'll now discuss the C++ implementation.

C++ Implementation

We've previously considered analytical pricing of European call options in some depth. However, I've reproduced and modified the code here for completeness (bs_prices.h). If you are interested in the details of the algorithms or approximations, check out the aforementioned link:

#ifndef __BS_PRICES_H
#define __BS_PRICES_H

#define _USE_MATH_DEFINES

#include <iostream>
#include <cmath>

// Standard normal probability density function
double norm_pdf(const double x) {
  return (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x);
}

// An approximation to the cumulative distribution function
// for the standard normal distribution
// Note: This is a recursive function
double norm_cdf(const double x) {
  double k = 1.0/(1.0 + 0.2316419*x);
  double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k))));

  if (x >= 0.0) {
    return (1.0 - (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x) * k_sum);
  } else {
    return 1.0 - norm_cdf(-x);
  }
}

// This calculates d_j, for j in {1,2}. This term appears in the closed
// form solution for the European call or put price
double d_j(const int j, const double S, const double K, const double r, const double sigma, const double T) {
  return (log(S/K) + (r + (pow(-1,j-1))*0.5*sigma*sigma)*T)/(sigma*(pow(T,0.5)));
}

// Calculate the European vanilla call price based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double call_price(const double S, const double K, const double r, const double sigma, const double T) {
  return S * norm_cdf(d_j(1, S, K, r, sigma, T))-K*exp(-r*T) * norm_cdf(d_j(2, S, K, r, sigma, T));
}

#endif

The following is a listing for black_scholes.h, which contains a basic function object (functor) for handling calculation of options prices when provided with a volatility, $\sigma$. Notice that all of the Black-Scholes model option paramaters are stored as private variables, with the exception of the volatility $\sigma$. This is because the function call operator() method takes it as a parameter. This method will eventually be (repeatedly) called by the interval_bisection function:

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

  double operator()(double sigma) const;
};
#endif

The following is a listing for black_scholes.cpp, which contains the source implementation for the Black-Scholes options class. This class simply provides a wrapper around the analytical price for the call option, but crucially specifies the needed parameters, in such a way that the interval_bisection function avoids having to concern itself with them. Notice that $S$, $K$, $r$ and $T$ are referencing private member data, while $\sigma$ is being passed from the method call as a parameter:

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

double BlackScholesCall::operator()(double sigma) const {
  return call_price(S, K, r, sigma, T);
}

#endif

The following is a listing for interval_bisection.h, which contains the function template for carrying out interval bisection. As with classes, we need to add the template<typename T> syntax to the signature of the function in order to tell the compiler it should be expecting a generic type as one of its parameters. The function body implicitly calls operator() of the function object $g$, so any object passed to it must define operator() for a sole parameter.

The remainder of the code carries out the Interval Bisection algorithm for finding a root of the generic function $g$:

#ifndef __INTERVAL_BISECTION_H
#define __INTERVAL_BISECTION_H

#include <cmath>

// Creating a function template
// Trying to find an x such that |g(x) - y| < epsilon,
// starting with the interval (m, n).
template<typename T>
double interval_bisection(double y_target,  // Target y value
                          double m,         // Left interval value
                          double n,         // Right interval value
                          double epsilon,   // Tolerance
                          T g) {            // Function object of type T, named g

  // Create the initial x mid-point value
  // Find the mapped y value of g(x)
  double x = 0.5 * (m + n);
  double y = g(x);

  // While the difference between y and the y_target
  // value is greater than epsilon, keep subdividing
  // the interval into successively smaller halves
  // and re-evaluate the new y.
  do {
    if (y < y_target) {
      m = x;
    }

    if (y > y_target) {
      n = x;
    }

    x = 0.5 * (m + n);
    y = g(x);
  } while (fabs(y-y_target) > epsilon);

  return x;
}

#endif

The final listing is for the main implementation (main.cpp). I've hardcoded some option parameters (but you could easily modify this to input them from the command line), so that the implied volatility can be calculated. Firstly we create a BlackScholesCall instance and pass it the required parameters (without the volatility, at this stage). Then we define the parameters for the interval bisection itself, and finally pass the interval_bisection the bsc function object as a generic parameter. This then calculates the implied volatility of the option for us:

#ifndef __MAIN_CPP
#define __MAIN_CPP

#include "black_scholes.h"
#include "interval_bisection.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);

  // Interval Bisection parameters
  double low_vol = 0.05;
  double high_vol = 0.35;
  double episilon = 0.001;

  // Calculate the implied volatility
  double sigma = interval_bisection(C_M, low_vol, high_vol, episilon, bsc);

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

  return 0;
}
#endif

The output of the code is as follows:

Implied Vol: 0.201318

Thus we obtain an implied volatility of 20.1% for this particular call option. The object-oriented and generic aspects of the above code lend themselves naturally to extension and re-use. In subsequent articles we will create a Newton-Raphson solver as well as a Brent's Method solver to calculate the same values, but in a significantly more optimal manner.