ceph/src/boost/libs/math/example/autodiff_black_scholes.cpp

   1 //           Copyright Matthew Pulver 2018 - 2019.
   2 // Distributed under the Boost Software License, Version 1.0.
   3 //      (See accompanying file LICENSE_1_0.txt or copy at
   4 //           https://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <boost/math/differentiation/autodiff.hpp>
   7 #include <iostream>
   8 #include <stdexcept>
   9
  10 using namespace boost::math::constants;
  11 using namespace boost::math::differentiation;
  12
  13 // Equations and function/variable names are from
  14 // https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
  15
  16 // Standard normal probability density function
  17 template <typename X>
  18 X phi(X const& x) {
  19   return one_div_root_two_pi<X>() * exp(-0.5 * x * x);
  20 }
  21
  22 // Standard normal cumulative distribution function
  23 template <typename X>
  24 X Phi(X const& x) {
  25   return 0.5 * erfc(-one_div_root_two<X>() * x);
  26 }
  27
  28 enum class CP { call, put };
  29
  30 // Assume zero annual dividend yield (q=0).
  31 template <typename Price, typename Sigma, typename Tau, typename Rate>
  32 promote<Price, Sigma, Tau, Rate> black_scholes_option_price(CP cp,
  33                                                             double K,
  34                                                             Price const& S,
  35                                                             Sigma const& sigma,
  36                                                             Tau const& tau,
  37                                                             Rate const& r) {
  38   using namespace std;
  39   auto const d1 = (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
  40   auto const d2 = (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
  41   switch (cp) {
  42     case CP::call:
  43       return S * Phi(d1) - exp(-r * tau) * K * Phi(d2);
  44     case CP::put:
  45       return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1);
  46     default:
  47       throw std::runtime_error("Invalid CP value.");
  48   }
  49 }
  50
  51 int main() {
  52   double const K = 100.0;  // Strike price.
  53   auto const variables = make_ftuple<double, 3, 3, 1, 1>(105, 5, 30.0 / 365, 1.25 / 100);
  54   auto const& S = std::get<0>(variables);      // Stock price.
  55   auto const& sigma = std::get<1>(variables);  // Volatility.
  56   auto const& tau = std::get<2>(variables);    // Time to expiration in years. (30 days).
  57   auto const& r = std::get<3>(variables);      // Interest rate.
  58   auto const call_price = black_scholes_option_price(CP::call, K, S, sigma, tau, r);
  59   auto const put_price = black_scholes_option_price(CP::put, K, S, sigma, tau, r);
  60
  61   double const d1 = static_cast<double>((log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau)));
  62   double const d2 = static_cast<double>((log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau)));
  63   double const formula_call_delta = +Phi(+d1);
  64   double const formula_put_delta = -Phi(-d1);
  65   double const formula_vega = static_cast<double>(S * phi(d1) * sqrt(tau));
  66   double const formula_call_theta =
  67       static_cast<double>(-S * phi(d1) * sigma / (2 * sqrt(tau)) - r * K * exp(-r * tau) * Phi(+d2));
  68   double const formula_put_theta =
  69       static_cast<double>(-S * phi(d1) * sigma / (2 * sqrt(tau)) + r * K * exp(-r * tau) * Phi(-d2));
  70   double const formula_call_rho = static_cast<double>(+K * tau * exp(-r * tau) * Phi(+d2));
  71   double const formula_put_rho = static_cast<double>(-K * tau * exp(-r * tau) * Phi(-d2));
  72   double const formula_gamma = static_cast<double>(phi(d1) / (S * sigma * sqrt(tau)));
  73   double const formula_vanna = static_cast<double>(-phi(d1) * d2 / sigma);
  74   double const formula_charm =
  75       static_cast<double>(phi(d1) * (d2 * sigma * sqrt(tau) - 2 * r * tau) / (2 * tau * sigma * sqrt(tau)));
  76   double const formula_vomma = static_cast<double>(S * phi(d1) * sqrt(tau) * d1 * d2 / sigma);
  77   double const formula_veta = static_cast<double>(-S * phi(d1) * sqrt(tau) *
  78                                                   (r * d1 / (sigma * sqrt(tau)) - (1 + d1 * d2) / (2 * tau)));
  79   double const formula_speed =
  80       static_cast<double>(-phi(d1) * (d1 / (sigma * sqrt(tau)) + 1) / (S * S * sigma * sqrt(tau)));
  81   double const formula_zomma = static_cast<double>(phi(d1) * (d1 * d2 - 1) / (S * sigma * sigma * sqrt(tau)));
  82   double const formula_color =
  83       static_cast<double>(-phi(d1) / (2 * S * tau * sigma * sqrt(tau)) *
  84                           (1 + (2 * r * tau - d2 * sigma * sqrt(tau)) * d1 / (sigma * sqrt(tau))));
  85   double const formula_ultima =
  86       -formula_vega * static_cast<double>((d1 * d2 * (1 - d1 * d2) + d1 * d1 + d2 * d2) / (sigma * sigma));
  87
  88   std::cout << std::setprecision(std::numeric_limits<double>::digits10)
  89             << "autodiff black-scholes call price = " << call_price.derivative(0, 0, 0, 0) << '\n'
  90             << "autodiff black-scholes put  price = " << put_price.derivative(0, 0, 0, 0) << '\n'
  91             << "\n## First-order Greeks\n"
  92             << "autodiff call delta = " << call_price.derivative(1, 0, 0, 0) << '\n'
  93             << " formula call delta = " << formula_call_delta << '\n'
  94             << "autodiff call vega  = " << call_price.derivative(0, 1, 0, 0) << '\n'
  95             << " formula call vega  = " << formula_vega << '\n'
  96             << "autodiff call theta = " << -call_price.derivative(0, 0, 1, 0)
  97             << '\n'  // minus sign due to tau = T-time
  98             << " formula call theta = " << formula_call_theta << '\n'
  99             << "autodiff call rho   = " << call_price.derivative(0, 0, 0, 1) << '\n'
 100             << " formula call rho   = " << formula_call_rho << '\n'
 101             << '\n'
 102             << "autodiff put delta = " << put_price.derivative(1, 0, 0, 0) << '\n'
 103             << " formula put delta = " << formula_put_delta << '\n'
 104             << "autodiff put vega  = " << put_price.derivative(0, 1, 0, 0) << '\n'
 105             << " formula put vega  = " << formula_vega << '\n'
 106             << "autodiff put theta = " << -put_price.derivative(0, 0, 1, 0) << '\n'
 107             << " formula put theta = " << formula_put_theta << '\n'
 108             << "autodiff put rho   = " << put_price.derivative(0, 0, 0, 1) << '\n'
 109             << " formula put rho   = " << formula_put_rho << '\n'
 110             << "\n## Second-order Greeks\n"
 111             << "autodiff call gamma = " << call_price.derivative(2, 0, 0, 0) << '\n'
 112             << "autodiff put  gamma = " << put_price.derivative(2, 0, 0, 0) << '\n'
 113             << "      formula gamma = " << formula_gamma << '\n'
 114             << "autodiff call vanna = " << call_price.derivative(1, 1, 0, 0) << '\n'
 115             << "autodiff put  vanna = " << put_price.derivative(1, 1, 0, 0) << '\n'
 116             << "      formula vanna = " << formula_vanna << '\n'
 117             << "autodiff call charm = " << -call_price.derivative(1, 0, 1, 0) << '\n'
 118             << "autodiff put  charm = " << -put_price.derivative(1, 0, 1, 0) << '\n'
 119             << "      formula charm = " << formula_charm << '\n'
 120             << "autodiff call vomma = " << call_price.derivative(0, 2, 0, 0) << '\n'
 121             << "autodiff put  vomma = " << put_price.derivative(0, 2, 0, 0) << '\n'
 122             << "      formula vomma = " << formula_vomma << '\n'
 123             << "autodiff call veta = " << call_price.derivative(0, 1, 1, 0) << '\n'
 124             << "autodiff put  veta = " << put_price.derivative(0, 1, 1, 0) << '\n'
 125             << "      formula veta = " << formula_veta << '\n'
 126             << "\n## Third-order Greeks\n"
 127             << "autodiff call speed = " << call_price.derivative(3, 0, 0, 0) << '\n'
 128             << "autodiff put  speed = " << put_price.derivative(3, 0, 0, 0) << '\n'
 129             << "      formula speed = " << formula_speed << '\n'
 130             << "autodiff call zomma = " << call_price.derivative(2, 1, 0, 0) << '\n'
 131             << "autodiff put  zomma = " << put_price.derivative(2, 1, 0, 0) << '\n'
 132             << "      formula zomma = " << formula_zomma << '\n'
 133             << "autodiff call color = " << call_price.derivative(2, 0, 1, 0) << '\n'
 134             << "autodiff put  color = " << put_price.derivative(2, 0, 1, 0) << '\n'
 135             << "      formula color = " << formula_color << '\n'
 136             << "autodiff call ultima = " << call_price.derivative(0, 3, 0, 0) << '\n'
 137             << "autodiff put  ultima = " << put_price.derivative(0, 3, 0, 0) << '\n'
 138             << "      formula ultima = " << formula_ultima << '\n';
 139   return 0;
 140 }
 141 /*
 142 Output:
 143 autodiff black-scholes call price = 56.5136030677739
 144 autodiff black-scholes put  price = 51.4109161009333
 145
 146 ## First-order Greeks
 147 autodiff call delta = 0.773818444921273
 148  formula call delta = 0.773818444921274
 149 autodiff call vega  = 9.05493427705736
 150  formula call vega  = 9.05493427705736
 151 autodiff call theta = -275.73013426444
 152  formula call theta = -275.73013426444
 153 autodiff call rho   = 2.03320550539396
 154  formula call rho   = 2.03320550539396
 155
 156 autodiff put delta = -0.226181555078726
 157  formula put delta = -0.226181555078726
 158 autodiff put vega  = 9.05493427705736
 159  formula put vega  = 9.05493427705736
 160 autodiff put theta = -274.481417851526
 161  formula put theta = -274.481417851526
 162 autodiff put rho   = -6.17753255212599
 163  formula put rho   = -6.17753255212599
 164
 165 ## Second-order Greeks
 166 autodiff call gamma = 0.00199851912993254
 167 autodiff put  gamma = 0.00199851912993254
 168       formula gamma = 0.00199851912993254
 169 autodiff call vanna = 0.0410279463126531
 170 autodiff put  vanna = 0.0410279463126531
 171       formula vanna = 0.0410279463126531
 172 autodiff call charm = -1.2505564233679
 173 autodiff put  charm = -1.2505564233679
 174       formula charm = -1.2505564233679
 175 autodiff call vomma = -0.928114149313108
 176 autodiff put  vomma = -0.928114149313108
 177       formula vomma = -0.928114149313107
 178 autodiff call veta = 26.7947073115641
 179 autodiff put  veta = 26.7947073115641
 180       formula veta = 26.7947073115641
 181
 182 ## Third-order Greeks
 183 autodiff call speed = -2.90117322380992e-05
 184 autodiff put  speed = -2.90117322380992e-05
 185       formula speed = -2.90117322380992e-05
 186 autodiff call zomma = -0.000604548369901419
 187 autodiff put  zomma = -0.000604548369901419
 188       formula zomma = -0.000604548369901419
 189 autodiff call color = -0.0184014426606065
 190 autodiff put  color = -0.0184014426606065
 191       formula color = -0.0184014426606065
 192 autodiff call ultima = -0.0922426864775683
 193 autodiff put  ultima = -0.0922426864775683
 194       formula ultima = -0.0922426864775685
 195 **/