Skip to content

Instantly share code, notes, and snippets.

@nesteruk
Created June 6, 2020 11:26
Show Gist options
  • Save nesteruk/ec3f39e5c6d8c609c5cf6d84971fe5a5 to your computer and use it in GitHub Desktop.
Save nesteruk/ec3f39e5c6d8c609c5cf6d84971fe5a5 to your computer and use it in GitHub Desktop.
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment