C++ - Coding Black Scholes Formula

Question: Write a function to calculate the value of a call option using the Black-Scholes formula given by

where

To complete this task we will need to use the normal distribution function as demonstrated on the previous post. 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 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;
}

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

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.;

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

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

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;
// calculate the value of the cumulative normal distribution using Simpson's method
double normalDistribution(double x)
{
  if(x<-10.)return 0.; // return sensible values on limits
  if(x>10.)return 1.;
  // number of steps
  int N=1000;
  // range of integration
  double a=0,b=x;
  // local variables
  double s,h,sum=0.;
  // initialise the variables
  h=(b-a)/N;
  // add in the first few terms
  sum = sum + exp(-a*a/2.) + 4.*exp(-(a+h)*(a+h)/2.);
  // and the last one
  sum = sum + exp(-b*b/2.);
  // loop over terms 2 up to N-1
  for(int i=1;i<N/2;i++)
  {
    s = a + 2*i*h;
    sum = sum + 2.*exp(-s*s/2.);
    s = s + h;
    sum = sum + 4.*exp(-s*s/2.);
  }
  // complete the integral
  sum = 0.5 + h*sum/3./sqrt(8.*atan(1.));
  // return result
  return sum;
}
// 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;
}

leftright-investment