// ==========================================================================
// Author:  Yee Hsu
// Date:    10/10/2012
//
// More Reference ...
//
// http://en.wikipedia.org/wiki/Black%E2%80%93Scholes
// http://www.espenhaug.com/black_scholes.html
// http://www.numa.com/cgi-bin/numa/calc_op.pl
// 
// This is only partial implementation of the Black Scholes Model.
// There are other ones such as ...
//    - Black
//    - Binomial
//    - Binomial Smooth
//    - Binomial N
//    - Garman Kolhagen
//    - etc
// ==========================================================================
 
#include "stdafx.h"
#include "Black-Scholes.h"
#include "math.h"
 
#include <vector>
 
void BlackScholes::AssignValue(OptionProperty& opt)
{
      S = opt.price;
      K = opt.strike;
      T = opt.time;
      r = opt.interestrate;
      v = opt.volatility;
}
 
double BlackScholes::f(double x)
{
      double pi = 4.0*atan(1.0);
      return exp(-x*x*0.5)/sqrt(2*pi);
}
 
double BlackScholes::N(double x)
{
      return Boole(-10.0, x, 240);
}
 
// Boole's rule of approximation
double BlackScholes::Boole(double StartPoint, double EndPoint, int n)
{
      std::vector<double> X(n+1, 0.0);
      std::vector<double> Y(n+1, 0.0);
      double delta_x = (EndPoint - StartPoint)/double(n);
 
      for (int i=0; i<=n; i++)
      {
            X[i] = StartPoint + i*delta_x;
            Y[i] = f(X[i]);
      }
 
      double sum = 0;
      for (int t=0; t<=(n-1)/4; t++)
      {
            int ind = 4*t;
            sum += (1/45.0)*(14*Y[ind] + 64*Y[ind+1] + 24*Y[ind+2] + 64*Y[ind+3] + 14*Y[ind+4])*delta_x;
      }
      return sum;
}
 
// The cumulative normal distribution function
double BlackScholes::Calc_CND( double X )
{
      double L, K, w ;
      double const a1 = 0.31938153, a2 = -0.356563782, a3 = 1.781477937;
      double const a4 = -1.821255978, a5 = 1.330274429;
 
      L = fabs(X);
      K = 1.0 / (1.0 + 0.2316419 * L);
      w = 1.0 - 1.0 / sqrt(2 * PI) * exp(-L *L / 2) * (a1 * K + a2 * K *K + a3 * pow(K,3) + a4 * pow(K,4) + a5 * pow(K,5));
 
      if (X < 0 )
      {
            w= 1.0 - w;
      }
      return w;
}
 
/*
      The Black and Scholes (1973) Stock option formula
      This is the theoretical (or fair) value of the option,
      and should be compared with the actual trading price of the option in the market.
 
      This Black-Scholes calculator allows you to figure out the value of a European call or put option.
      The calculator uses the stock's current share price, the option strike price, time to expiration,
      risk-free interest rate, and volatility to derive the value of these options.
      The Black-Scholes calculation used by this tool assumes no dividend is paid on the stock.
*/
double BlackScholes::Calc_BlackScholes(OptionProperty& opt)
{
      AssignValue(opt);
 
      double d1=(log(S/K)+(r+v*v/2)*T)/(v*sqrt(T));
      double d2=d1-v*sqrt(T);
 
      if(opt.type == CALL)
            return S *Calc_CND(d1)-K * exp(-r*T)*Calc_CND(d2);
      else if(opt.type == PUT)
            return K * exp(-r * T) * Calc_CND(-d2) - S * Calc_CND(-d1);
}
 
// The Greeks measurements
/*
      Delta = change in OPTION price with respect to change in Stock price.
      if stock price changes, how sensitive will BS be. All Greeks are first derivatives. If
      delta is .6, if stock price changes so will the option price but by 60% of the stock price
      change. Delta is a sensitivity for the stock price. Delta is a percentage between 0 - 1.0.
      As it converges to 1, the stock price increases faster. As stock prices lower, option
      prices converges out of the money.
*/
double BlackScholes::Calc_Delta(OptionProperty& opt)
{
      AssignValue(opt);
 
      double d = (log(S/K) + T*(r + 0.5*v*v)) / (v*sqrt(T));
      if (opt.type == CALL)
            return N(d);
      else
            return N(d) - 1;
}
 
/*
      Gamma=Change ins DELTA with respect to change in stock price. Gamma is second derivative of
      stock price but first derivative of delta. Longer term of gamma of call option vs time to
      expire also indicates lower delta. Delta peaks at the money when the option is at the
      money. This is where the delta is changing the most rapidly. After the peak, the delta
      converges to 1 but becoming more stable. Delta converges more stable as it converges from 0
      and the rate of change is slower.
*/
double BlackScholes::Calc_Gamma(OptionProperty& opt)
{
      AssignValue(opt);
 
      double d = (log(S/K) + T*(r + 0.5*v*v)) / (v*sqrt(T));
      return f(d) / S / v / sqrt(T);
}
 
/*
      Vega=Change in OPTION price with respect to CHANGE in volatility. This is significant since
      the sensitivity is high. The higher vega means the option price is change (or more
      sensitive) to changes in volatility. The longer the term, the more sensitive as vega grows
      as the option life grows. Gamma peaks at the in the money.
*/
double BlackScholes::Calc_Vega(OptionProperty& opt)
{
      AssignValue(opt);
 
      double d = (log(S/K) + T*(r + 0.5*v*v)) / (v*sqrt(T));
      return S*f(d)*sqrt(T);
}
 
/*
      Rho=Change in OPTION price with respect to change in interest rate. It is less sensitive as
      compare to other inputs of volatility. Option price is more sensitive to changes in
      volatility versus interest rate.
*/
double BlackScholes::Calc_Rho(OptionProperty& opt)
{
      AssignValue(opt);
 
      double d = (log(S/K) + T*(r + 0.5*v*v)) / (v*sqrt(T));
      if (opt.type == CALL)
            return T*K*exp(-r*T)*N(d - v*sqrt(T));
      else
            return -T*K*exp(-r*T)*N(v*sqrt(T) - d);
}
 
/*
      Theta=Change in OPTION price with respect to change in Term (T). This is known as time
      decay. This is always negative. As time passes, the option and gets closer to expiration,
      it becomes less valueable. Time erodes the value of the option.
*/
double BlackScholes::Calc_Theta(OptionProperty& opt)
{
      AssignValue(opt);
 
      double d = (log(S/K) + T*(r + 0.5*v*v)) / (v*sqrt(T));
      if (opt.type == CALL)
            return -S*f(d)*v/2/sqrt(T) - r*K*exp(-r*T)*N(d - v*sqrt(T));
      else
            return -S*f(d)*v/2/sqrt(T) + r*K*exp(-r*T)*N(v*sqrt(T) - d);
}
 
/*
int _tmain(int argc, _TCHAR* argv[])
{
      IDerivativeModel* pDev = new BlackScholes();
      OptionProperty opt;
 
      // Option properties
      opt.price = 60.0;             // stock price
      opt.strike = 65.0;                  // strike
      opt.time = 0.25;              // 3 month
      opt.interestrate = 0.08;      // 8%
      opt.volatility = 0.3;         // 30%
      opt.type = CALL;
 
      // Derivative Measurements   
      printf("Black Scholes Delta = %f\n", pDev->Calc_Delta(opt));
      printf("Black Scholes Gamma = %f\n", pDev->Calc_Gamma(opt));
      printf("Black Scholes Vega  = %f\n", pDev->Calc_Vega(opt));
      printf("Black Scholes Rho   = %f\n", pDev->Calc_Rho(opt));
      printf("Black Scholes Theta = %f\n", pDev->Calc_Theta(opt));
 
      // Black-Scholes Specific
      BlackScholes* pbs = static_cast<BlackScholes*>(pDev);
      printf("Black Scholes Value = %f\n", pbs->Calc_BlackScholes(opt));
 
      return 0;
}
*/