C++ - Coding Black Scholes Formula

Question: Write a C++ code to calculate the option values for a call C(t = 0,S), and a put P(t = 0,S) using the Black-Scholes formula:

(2)

(3)

with

(4)

(5)

You can use the cumulative distribution function N(x) from earlier on. So rather than starting with an empty project copy and paste the code from the previous project into your new project, you can find a version of the code on my website (click here). Now we need to create a function for the call option, so we need to think about what arguments need to be supplied and what are the local variables to the function. Obviously we need to supply the asset price, current time both of which may vary, and the parameters for the function are the strike price, interest rate, volatility and maturity, and we want to return a real number as the answer. Inside the function, we will need to calculate the value of d1 and d2. This leads to a function definition of the form

double callOptionPrice(double S,double t,double X,double r,double sigma,double T)
{
  double d1;
  double d2;
  return 0;
}

download

You must place this function header in between the main function and the normalDistribution function, since we will be calling ”normalDistribution” inside ”callOptionPrice” and then ”callOptionPrice” inside ”main”. Next we can simply fill in the mathematical formulas as follows

double callOptionPrice(double S,double t,double X,double r,double sigma,double T)
{
  double d1=(log(S/X) + (r+sigma*sigma/2.)*(T-t))/(sigma*sqrt(T-t));
  double d2=(log(S/X) + (r-sigma*sigma/2.)*(T-t))/(sigma*sqrt(T-t));
  return normalDistribution(d1)*S - normalDistribution(d2)*X*exp(-r*(T-t));
}

int main()
{
  cout << "Call Option Price = " << callOptionPrice(1,0,1,0.05,0.2,1) << endl;
  return 0;
}

download

Now we need to test the function under different settings. There are obviously going to be problem in each of the following cases:

• S=0
• sigma=0
• t=T

since all will result in undefined mathematical values given the way that d1 and d2 are calculated. However, returning a value of plus or minus infinity for d1 and d2 does not result in an undefined value for the cumulative normal distribution since the function returns finite values from infinite limits. Using our knowledge of what happens to this function we can go through each case and decide what are the appropriate values to return.

Case S=0

In this situation we must make use of the boundary condition of the problem, which is that the option value is worthless if the asset value is zero. The only slight caveat here is that you should check whether S is smaller than a very small number rather than comparing it with zero. In code you add the following line to your function before the calculations.

if(S<1.e-14)return 0.;

download

By executing a return the function will stop and return the value without executing any more code.

Case sigma=0

This case is very similar to when t=T. Depending on the sign of the numerator in the calculation for d1 and d2 the function will either return zero or the asset minus the discounted strike price. In code this looks like

if(sigma<1.e-14)
  {
    if(S<X*exp(-r*(T-t)))return 0.;
    else return S-X*exp(-r*(T-t));
  }

download

Case t=T

Finally if t and T are almost equal then we are at maturity and we return the payoff. The code might look like

  if(fabs(T-t)<1.e-14)
  {
    if(S<X)return 0.;
    else return S-X;
  }

download

On adding each of these parts to the code you should test and validate each part using a range of parameters. Note here that another case that could cause problems is if t > T. There is no sensible value that can be returned in this case so if you add in a check for it you would be looking to exit the program if this happened or at least print a warning to the screen. Adding this in is left as an exercise. Your final code should look like

#include <iostream>
#include <cmath>
using namespace std;

double normalDistribution(double x)
{
  static const double RT2PI = sqrt(4.0*acos(0.0));
  static const double SPLIT = 10./sqrt(2);
  static const double a[] = {220.206867912376,221.213596169931,112.079291497871,33.912866078383,6.37396220353165,0.700383064443688,3.52624965998911e-02};
  static const double b[] = {440.413735824752,793.826512519948,637.333633378831,296.564248779674,86.7807322029461,16.064177579207,1.75566716318264,8.83883476483184e-02};
  
  const double z = fabs(x);
  double Nz = 0.0;
  
  // if z outside these limits then value effectively 0 or 1 for machine precision
  if(z<=37.0)
  {
    // NDash = N'(z) * sqrt{2\pi}
    const double NDash = exp(-z*z/2.0)/RT2PI;
    if(z<SPLIT)
    {
      const double Pz = (((((a[6]*z + a[5])*z + a[4])*z + a[3])*z + a[2])*z + a[1])*z + a[0];
      const double Qz = ((((((b[7]*z + b[6])*z + b[5])*z + b[4])*z + b[3])*z + b[2])*z + b[1])*z + b[0];
      Nz = RT2PI*NDash*Pz/Qz;
    }
    else
    {
      const double F4z = z + 1.0/(z + 2.0/(z + 3.0/(z + 4.0/(z + 13.0/20.0))));
      Nz = NDash/F4z;
    }
  }
  return x>=0.0 ? 1-Nz : Nz;
}

// return the value of a call option using the black scholes formula
double callOptionPrice(double S,double t,double X,double r,double sigma,double T)
{
  if(S<1.e-14)return 0.; // check if asset worthless
  if(sigma<1.e-14) // check if sigma zero
  {
    if(S<X*exp(-r*(T-t)))return 0.;
    else return S-X*exp(-r*(T-t));
  }
  if(fabs(T-t)<1.e-14) // check if we are at maturity
  {
    if(S<X)return 0.;
    else return S-X;
  }
  // calculate option price
  double d1=(log(S/X) + (r+sigma*sigma/2.)*(T-t))/(sigma*sqrt(T-t));
  double d2=(log(S/X) + (r-sigma*sigma/2.)*(T-t))/(sigma*sqrt(T-t));
  return normalDistribution(d1)*S - normalDistribution(d2)*X*exp(-r*(T-t));
}

int main()
{
  cout << "Call Option Price = " << callOptionPrice(1,0,1,0.05,0.2,1) << endl;
  return 0;
}

download