European vanilla option pricing with C++ and analytic formulae

European vanilla option pricing with C++ and analytic formulae

In this article we will price a European vanilla option via the correct analytic solution of the Black-Scholes equation. We won't be concentrating on an extremely efficient or optimised implementation at this stage. Right now I just want to show you how the mathematical formulae correspond to the C++ code.

Black-Scholes Analytic Pricing Formula

The first stage in implementation is to briefly discuss the Black-Scholes analytic solution for the price of a vanilla call or put option. Consider the price of a European Vanilla Call, $C(S,t)$. $S$ is the underlying asset price, $K$ is the strike price, $r$ is the interest rate (or the "risk-free rate"), $T$ is the time to maturity and $\sigma$ is the (constant) volatility of the underlying asset $S$. $N$ is a function which will be described in detail below. The analytical formula for $C(S,t)$ is given by:

\begin{eqnarray*} C(S,t) = SN(d_1) - Ke^{-rT} N(d_2) \end{eqnarray*}

With $d_1$ and $d_2$ defined as follows:

\begin{eqnarray*} d_1 &=& \frac{log(S/K) + (r+\frac{\sigma^2}{2})T}{\sigma \sqrt{T} }\\ d_2 &=& d_1 - \sigma \sqrt{T} \end{eqnarray*}

Making use of put-call parity, we can also price a European vanilla put, $P(S,t)$, with the following formula:

\begin{eqnarray*} P(S,t) = Ke^{-rT} - S + C(S,t) = Ke^{-rT} - S + (SN(d_1) - Ke^{-rT} N(d_2)) \end{eqnarray*}

All that remains is to describe the function $N$, which is the cumulative distribution function of the standard normal distribution. The formula for $N$ is given by:

\begin{eqnarray*} N(x) = \frac{1}{\sqrt{2 \pi}} \int^x_{-\infty} e^{-t^{2} /2} dt \end{eqnarray*}

It would also help to have closed form solutions for the "Greeks". These are the sensitivities of the option price to the various underlying parameters. In order to calculate these sensitivities we need the formula for the probability density function of the standard normal distribution which is given below:

\begin{eqnarray*} f(x) = \frac{1}{\sqrt{2 \pi}} e^{-x^{2}/2} \end{eqnarray*}

Note that we won't be computing the "Greeks" in this article, we will leave them for later, but I wanted you to be aware of the statistical formulae which are necessary.

C++ Implementation

I'll assume you're familiar with a C++ IDE such as Microsoft Visual Studio/C++, Eclipse, Emacs or Vim. This code will not consist of any object oriented or generic programming techniques at this stage. Right now the goal is to help you understand a basic C++ implementation, without all of the additional object machinery to encapsulate away the mathematical functions. I've written out the program in full and then below I'll explain how each component works subsequently:

#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& v, const double& T) {
    return (log(S/K) + (r + (pow(-1,j-1))*0.5*v*v)*T)/(v*(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& v, const double& T) {
    return S * norm_cdf(d_j(1, S, K, r, v, T))-K*exp(-r*T) * norm_cdf(d_j(2, S, K, r, v, T));
}

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

int main(int argc, char **argv) {
    // First we create the parameter list
    double S = 100.0;  // Option price
    double K = 100.0;  // Strike price
    double r = 0.05;   // Risk-free rate (5%)
    double v = 0.2;    // Volatility of the underlying (20%)
    double T = 1.0;    // One year until expiry

    // Then we calculate the call/put values
    double call = call_price(S, K, r, v, T);
    double put = put_price(S, K, r, v, T);

    // Finally we output the parameters and prices
    std::cout << "Underlying:      " << S << std::endl;
    std::cout << "Strike:          " << K << std::endl;
    std::cout << "Risk-Free Rate:  " << r << std::endl;
    std::cout << "Volatility:      " << v << std::endl;
    std::cout << "Maturity:        " << T << std::endl;

    std::cout << "Call Price:      " << call << std::endl;
    std::cout << "Put Price:       " << put << std::endl;

    return 0;
}

Let's now take a look at the program step-by-step.

The first line is a pre-processor macro which tells the C++ compiler to make use of the C-standard mathematical constants. I've written an article on the different mathematical constants you can use in C++. Note that some compilers may not fully support these constants. Make sure you test them out first!

#define _USE_MATH_DEFINES

The next section imports the iostream and cmath libraries. This allows us to use the std::cout command to output code. In addition we now have access to the C mathematical functions such as exp, pow, log and sqrt:

#include <iostream>
#include <cmath>

Once we have the correct libraries imported we need to create the core statistics functions which make up most of the computation for the prices. Here is the standard normal probability density function. M_PI is the C-standard mathematical constant for $\pi$:

// 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);
}

The next statistics function is the approximation to the cumulative distribution function for the standard normal distribution. This approximation is found in Joshi. Note that this is a recursive function (i.e. it calls itself!):

// 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);
    }
}

The last remaining set of functions are directly related to the pricing of European vanilla calls and puts. We need to calculate the $d_j$ values first, otherwise we would be repeating ourselves by adding the function code directly to each of call_price and put_price:

// 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& v, const double& T) {
    return (log(S/K) + (r + (pow(-1,j-1))*0.5*v*v)*T)/(v*(pow(T,0.5)));
}

Now that we have the cumulative distribution function for the standard normal distribution (norm_cdf) coded up, as well as the $d_j$ function, we can calculate the closed-form solution for the European vanilla call price. You can see that it is a fairly simple formula, assuming that the aforementioned functions have already been implemented:

// 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& v, const double& T) {
    return S * norm_cdf(d_j(1, S, K, r, v, T))-K*exp(-r*T) * norm_cdf(d_j(2, S, K, r, v, T));
}

And similarly for the vanilla put (this formula can be derived easily via put-call parity):

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

Every C++ program must have a main program. This is where the functions described above are actually called and the values derived are output to the console. We make use of the iostream library to provide us with the std::cout and std::endl output handlers. Finally, we exit the program:

int main(int argc, char **argv) {
    // First we create the parameter list
    double S = 100.0;  // Option price
    double K = 100.0;  // Strike price
    double r = 0.05;   // Risk-free rate (5%)
    double v = 0.2;    // Volatility of the underlying (20%)
    double T = 1.0;    // One year until expiry

    // Then we calculate the call/put values
    double call = call_price(S, K, r, v, T);
    double put = put_price(S, K, r, v, T);

    // Finally we output the parameters and prices
    std::cout << "Underlying:      " << S << std::endl;
    std::cout << "Strike:          " << K << std::endl;
    std::cout << "Risk-Free Rate:  " << r << std::endl;
    std::cout << "Volatility:      " << v << std::endl;
    std::cout << "Maturity:        " << T << std::endl;

    std::cout << "Call Price:      " << call << std::endl;
    std::cout << "Put Price:       " << put << std::endl;

    return 0;
}

The output of the code on my Mac OSX system is as follows:

Underlying:      100
Strike:          100
Risk-Free Rate:  0.05
Volatility:      0.2
Maturity:        1
Call Price:      10.4506
Put Price:       5.57352

This code should give you a good idea of how closed form solutions to the Black-Scholes equations can be coded up in a procedural manner with C++. The next steps are to calculate the "Greeks" in the same vein, as closed form solutions exist, solutions for digital and power options, as well as a basic Monte Carlo pricer with which to validate against.