#include <iostream>
#include <cmath>
#include <random>
using namespace std;
// Normal CDF using Error Function, erf
double norm_cdf(double x)
{
return 0.5 * (1 + erf(x / sqrt(2)));
}
// Enum for option type
enum OptionType { Call, Put };
// Black-Scholes formula for European options
// s0 = initial stock price, k = strike, r = interest rate
// v = implied volatility, t = time to expiry (years)
double blackscholes(OptionType type, double s0, double k, double r, double v, double t)
{
double d1 = (log(s0/k) + (r + 0.5*v*v) * t) / (v*sqrt(t));
double d2 = d1 - v * sqrt(t);
double sign = (type == Call) ? 1.0 : -1.0;
double nd1 = norm_cdf(sign * d1);
double nd2 = norm_cdf(sign * d2);
return sign * (s0 * nd1 - k * exp(-r*t) * nd2);
}
// Monte Carlo pricing of European option
double montecarlo(OptionType type, double s0, double k, double r, double v, double t, int nSims = 10000, bool show = true)
{
// Mersenne Twister RNG with fixed seed
mt19937 generator(123456);
normal_distribution<double> normsdist(0, 1);
double sign = (type == Call) ? 1.0 : -1.0;
double price = 0.0;
// Run Monte Carlo Simulations
for (int i = 1; i <= nSims; ++i)
{
// Generate Standard Normal Random Number, z
double z = normsdist(generator);
// Use GBM for Stock Price at Time t
double st = s0 * exp((r - 0.5*v*v) * t + v*sqrt(t)*z);
double payoff = max(sign * (st - k), 0.0);
// Martingale Price equals Discounted Payoff
price += payoff * exp(-r*t);
// Show Monte Carlo Intermmediate Results
if (show && i>0 && (i % 1000 == 0))
{
printf("%s%i%s%.6f\n", "nSims: ", i," price: ", price/i);
}
}
// Average Price
return price / nSims;
}
int main()
{
// European Call, S=100, K=100, r=1%, vol=10%, t=1 Year, Expected Price: 4.485236
double bs = blackscholes(Call, 100, 100, 0.01, 0.1, 1.0);
double mc = montecarlo(Call, 100, 100, 0.01, 0.1, 1.0, 10000);
printf("%s%.6f\n", "Monte Carlo Price: ", mc);
printf("%s%.6f\n", "Black-Scholes Price: ", bs);
return 0;
}
