Black-Scholes-Merton Options Pricing Engine in C++

Hi again everybody,
here's also a C++ class for pricing options premiums (Black-Scholes-Merton).

Code:

/*
  BSM.cpp - Black-Scholes-Merton options pricing engine

  2016-01-26-Tu: v0.99: init
  Author: U.M. in Germany (user botpro at www.elitetrader.com)
 
  Compile:
    g++ -Wall -O2 -std=c++11 BSM.cpp -o BSM.exe
 
  Run:
    Linux/Unix: ./BSM.exe
    Windows:    BSM.exe

    Tested on 64bit Linux only

  Example runs on Linux:
    ./BSM.exe 100 100 30 42
       Spot=100.00 Strike=100 Vola%=30.00 ExpDays=42.00(t=0.16667) AnnDriftPct%=0.00 AnnDividPct%=0.00 --> Call=4.88297 Put=4.88297

    ./BSM.exe 100 100 30 63
       Spot=100.00 Strike=100 Vola%=30.00 ExpDays=63.00(t=0.25000) AnnDriftPct%=0.00 AnnDividPct%=0.00 --> Call=5.97853 Put=5.97853

    ./BSM.exe 100 100 30 252
       Spot=100.00 Strike=100 Vola%=30.00 ExpDays=252.00(t=1.00000) AnnDriftPct%=0.00 AnnDividPct%=0.00 --> Call=11.92354 Put=11.92354
   
    ./BSM.exe 100 100 30 504
       Spot=100.00 Strike=100 Vola%=30.00 ExpDays=504.00(t=2.00000) AnnDriftPct%=0.00 AnnDividPct%=0.00 --> Call=16.79960 Put=16.79960

  Misc:
    - each invocation of the calc() function calculates Call.Value and Put.Value
    - uses Year=252 trading days, this can be changed via constructor parameter
    - extracted for standalone use from my BlackScholes class named TCBS
*/

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>

using namespace std;


#define ZEROADJ  1e-18


double normal_cdf(const double value)
{
   return 0.5 * erfc(-value * M_SQRT1_2);
}


class BSM
  { // Black-Scholes-Merton options pricing engine

    public:
      double S, K, s, t, r, q, dbDaysInYear;
      double u,  st, rt, qt, ut;
      double z,  z0;
      double z1, z2;  // ie. d1, d2 in orig BS
      double p1, p2;

      struct TS
        {
          double Value;    // the calculated premium will be stored here
         
       // The Greeks (omitted here):
       // double Delta, Vega, Theta, Rho, Gamma, Vanna, Charm, Speed, Zomma, Color, DVegaDTime, Vomma,
       //        DualDelta, DualGamma, DualTheta, Ultima, Lambda, ForwardDelta, DriftlessDelta, Cross;
         
          TS()
            {
              Value = 0.0;
              //...
            }
        };
      TS Call, Put;

      BSM()
        {
          S = K = s = t = r = q = dbDaysInYear = u = st = rt = qt = ut = 0;
          z = z0 = z1 = z2 = p1 = p2 = 0;
        }
       
    public:
      void calc(const double AdbSpot_S, const double AdbStrike_K, const double AdbAnnVolaPct, const double AdbExpireDays,
                const double AdbAnnDriftPct = 0.0, const double AdbAnnDividPct = 0.0, const double AdbDaysInYear = 252.0)
        {
          dbDaysInYear = AdbDaysInYear;
       
          S = AdbSpot_S;
          K = AdbStrike_K;
          s = AdbAnnVolaPct  / 100.0;
         
          t = AdbExpireDays  / dbDaysInYear;
          t = max(t, ZEROADJ);    // safety
         
          r = AdbAnnDriftPct / 100.0;
          q = AdbAnnDividPct / 100.0;

          st = s * sqrt(t);
          rt = r * t;
          qt = q * t;
          u  = r - q;
          ut = u * t;

#if 0
          // orig BS d1, d2:
          z1 = (log(S / K) + (u + s * s / 2.0) * t) / st;
          z2 = z1 - st;
#else
          // my unpublished method using z1, z2 (makes Ito's lemma unnecessary): 100% identical to d1, d2 of orig BS:
          z  = st / 2.0;
          z0 = (log(S / K) + ut) / st;
          z1 = z0 + z;
          z2 = z0 - z;
#endif

          // calc the corrosponding probabilities p1, p2 belonging to z1, z2:
          p1 = normal_cdf(z1);
          p2 = normal_cdf(z2);

          // calc Call premium:
          const double exp_mqt = exp(-qt), exp_mrt = exp(-rt);
          Call.Value = S * exp_mqt * p1      - K * exp_mrt * p2;
          if (Call.Value < 0.0) Call.Value = ZEROADJ;   // safety

          // calc Put premium: via Put/Call Parity:
          Put.Value = Call.Value + K * exp_mrt - S * exp_mqt;
          if (Put.Value < 0.0) Put.Value = ZEROADJ;     // safety
        }
  };


int main(int argc, char* argv[])
  {
    if (argc < 5) { printf("%s dbSpot dbStrike dbAnnVola%% dbExpDays\n", argv[0]); return 1; }

    const double dbSpot       = atof(argv[1]);
    const double dbStrike     = atof(argv[2]);
    const double dbAnnVolaPct = atof(argv[3]);
    const double dbExpDays    = atof(argv[4]);
 
    BSM B;

    B.calc(dbSpot, dbStrike, dbAnnVolaPct, dbExpDays);
    printf("Spot=%.2f Strike=%.0f Vola%%=%.2f ExpDays=%.2f(t=%.5f) AnnDriftPct%%=%.2f AnnDividPct%%=%.2f --> Call=%.5f Put=%.5f\n",
            B.S,      B.K,        B.s * 100.0, B.t * B.dbDaysInYear, B.t,  B.r * 100.0, B.q * 100.0,     B.Call.Value, B.Put.Value);

    return 0;   
  }

 

in

double normal_cdf(const double value)
{
return 0.5 * erfc(-value * M_SQRT1_2);
}

M_SQRT1_2 is not defined

where is erfc declared?

are we missing header file? This is VS2013

Thank you

 

The above one works in GNU g++.

Try by replacing the above one it with this one below.
Your C++ compiler should understand C++11. The last Microsoft VS I used was VS2008, so I don't know if VS2013 has the erf() function available, but I think it should have, so the following should work:

// Returns the probability of [-inf,x] of a gaussian distribution
double normal_cdf(const double x, const double mu = 0.0, const double sigma = 1.0)
{
return 0.5 * (1.0 + erf((x - mu) / (sigma * sqrt(2.0))));
}

 

Here's an updated BSM.cpp with the above changes:

 

Great, it works,thank you

 

Ok, here's hopefully the final version, no changes in code, just some cosmetics, and added one more example.