BlackScholesSimd.cpp

#define _USE_MATH_DEFINES

#include <iostream>
#include <cmath>

#include <immintrin.h>

// Standard normal probability density function
double norm_pdf(const double& x) {
  return (1.0 / (pow(2 * M_PI, 0.5))) * exp(-0.5 * x * x);
}



__m256d norm_pdf(const __m256d& x)
{
  const static double mul = 1.0 / (pow(2 * M_PI, 0.5));
  const static __m256d mul_vector = {mul,mul,mul,mul};
  const static __m256d minus_half = {-0.5,-0.5,-0.5,-0.5};
  
  return _mm256_mul_pd(
    mul_vector,
    _mm256_exp_pd(
      _mm256_mul_pd(minus_half,
        _mm256_mul_pd(x, x)
      )
    )
  );
}

// An approximation to the cumulative distribution function
// for the standard normal distribution
// Note: This is a recursive function
double norm_cdf(const double& x) {
  double k = 1.0 / (1.0 + 0.2316419 * x);
  double k_sum = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + 1.330274429 * k))));

  if (x >= 0.0) {
    return (1.0 - (1.0 / (pow(2 * M_PI, 0.5))) * exp(-0.5 * x * x) * k_sum);
  }
  else {
    return 1.0 - norm_cdf(-x);
  }
}

// This calculates d_j, for j in {1,2}. This term appears in the closed
// form solution for the European call or put price
double d_j(const int& j, const double& S, const double& K, const double& r, const double& v, const double& T) {
  return (log(S / K) + (r + (pow(-1, j - 1)) * 0.5 * v * v) * T) / (v * (pow(T, 0.5)));
}

__m256d d1(const __m256d& s, const __m256d& k, 
  const __m256d& r, const __m256d& v, const __m256d& t)
{
  const static __m256d half = {0.5,0.5,0.5,0.5};

  return _mm256_mul_pd(
    _mm256_add_pd(
     _mm256_log_pd(
       _mm256_div_pd(s, k)
     ),
     _mm256_add_pd(r,
                   _mm256_mul_pd(half, _mm256_mul_pd(v, v))
     )
   ), _mm256_invsqrt_pd(
     _mm256_mul_pd(_mm256_mul_pd(v, v), t)
   ));
}

__m256d d2(const __m256d& s, const __m256d& k, 
  const __m256d& r, const __m256d& v, const __m256d& t)
{
  return _mm256_sub_pd(
    d1(s, k, r, v, t),
    _mm256_mul_pd(v, _mm256_sqrt_pd(t)));
}

__m256d call_price(const __m256d& s, const __m256d& k, 
  const __m256d& r, const __m256d& v, const __m256d& t)
{
  const auto D1 = d1(s, k, r, v, t);
  const auto D2 = d2(s, k, r, v, t);

  auto s_n_d1 = _mm256_mul_pd(s, _mm256_cdfnorm_pd(D1));

  auto e_min_rt = _mm256_exp_pd(
    _mm256_xor_pd(_mm256_mul_pd(r, t), _mm256_set1_pd(-0.0))
  );

  return _mm256_sub_pd(
    s_n_d1,
    _mm256_mul_pd(
      _mm256_mul_pd(k, e_min_rt),
      _mm256_cdfnorm_pd(D2)
    )
  );
}

__m256d put_price(const __m256d& s, const __m256d& k, 
  const __m256d& r, const __m256d& v, const __m256d& t)
{
  auto call = call_price(s, k, r, v, t);
  auto e_min_rt = _mm256_exp_pd(
    _mm256_xor_pd(_mm256_mul_pd(r, t), _mm256_set1_pd(-0.0))
  );

  return
    _mm256_add_pd(
      _mm256_sub_pd(_mm256_mul_pd(k, e_min_rt), s),
      call);
}

// Calculate the European vanilla call price based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double call_price(const double& S, const double& K, const double& r, const double& v, const double& T) {
  return S * norm_cdf(d_j(1, S, K, r, v, T))
  - K * exp(-r * T) * norm_cdf(d_j(2, S, K, r, v, T));
}

// Calculate the European vanilla put price based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double put_price(const double& S, const double& K, const double& r, const double& v, const double& T) {
  return -S * norm_cdf(-d_j(1, S, K, r, v, T))
  + K * exp(-r * T) * norm_cdf(-d_j(2, S, K, r, v, T));
}

int main(int argc, char** argv)
{
  using namespace std;

  // First we create the parameter list
  double S = 100.0;  // Option price
  double K = 100.0;  // Strike price
  double r = 0.05;   // Risk-free rate (5%)
  double v = 0.2;    // Volatility of the underlying (20%)
  double T = 1.0;    // One year until expiry

  // Then we calculate the call/put values
  double call = call_price(S, K, r, v, T);
  double put = put_price(S, K, r, v, T);

  // Finally we output the parameters and prices
  std::cout << "Underlying:      " << S << std::endl;
  std::cout << "Strike:          " << K << std::endl;
  std::cout << "Risk-Free Rate:  " << r << std::endl;
  std::cout << "Volatility:      " << v << std::endl;
  std::cout << "Maturity:        " << T << std::endl;

  std::cout << "Call Price:      " << call << std::endl;
  std::cout << "Put Price:       " << put << std::endl;

  __m256d ss = {100,100,100,100};
  __m256d kk = {100,100,100,100};
  __m256d rr = {0.05,0.05,0.05,0.05};
  __m256d vv = {0.2,0.2,0.2,0.2};
  __m256d tt = {1,1,1,1};

  auto call_scalar = call_price(S, K, r, v, T);
  auto put_scalar =  put_price(S, K, r, v, T);

  auto call_vector = call_price(ss, kk, rr, vv, tt);
  auto put_vector = put_price(ss, kk, rr, vv, tt);

  cout << "Scalar call = " << call_scalar <<
    ", put = " << put_scalar << "\n";

  cout << "Vector call = " << call_vector.m256d_f64[0] <<
    ", put = " << put_vector.m256d_f64[0] << "\n";

  return 0;
}