std::normal_distribution

.gitignore

@@ -6,4 +6,5 @@ compile_commands.json
 
 #ide
 .idea
-
+test.cpp
+a.out

examples/monte-carlo-calculator/main.cpp

@@ -1,48 +1,29 @@
 #include <aws/lambda-runtime/runtime.h>
 #include <iomanip>
 #include <sstream>
-#include <algorithm>    // Needed for the "max" function
+#include <algorithm>
 #include <cmath>
 #include <iostream>
+#include <random>
 #include <aws/core/utils/json/JsonSerializer.h>
 #include <aws/core/utils/memory/stl/SimpleStringStream.h>
 
 using namespace aws::lambda_runtime;
 using namespace Aws::Utils::Json;
 
-
-
-// A simple implementation of the Box-Muller algorithm, used to generate
-// gaussian random numbers - necessary for the Monte Carlo method below
-// Note that C++11 actually provides std::normal_distribution<> in 
-// the <random> library, which can be used instead of this function
-double gaussian_box_muller() {
-  double x = 0.0;
-  double y = 0.0;
-  double euclid_sq = 0.0;
-
-  // Continue generating two uniform random variables
-  // until the square of their "euclidean distance" 
-  // is less than unity
-  do {
-    x = 2.0 * rand() / static_cast<double>(RAND_MAX)-1;
-    y = 2.0 * rand() / static_cast<double>(RAND_MAX)-1;
-    euclid_sq = x*x + y*y;
-  } while (euclid_sq >= 1.0);
-
-  return x*sqrt(-2*log(euclid_sq)/euclid_sq);
-}
-
 // Pricing a European vanilla call option with a Monte Carlo method
 double monte_carlo_call_price(const int& num_sims, const double& S, const double& K, const double& r, const double& v, const double& T) {
-  double S_adjust = S * exp(T*(r-0.5*v*v));
-  double S_cur = 0.0;
+  double s_adjust = S * exp(T*(r-0.5*v*v));
+  double s_cur = 0.0;
   double payoff_sum = 0.0;
 
+  std::random_device rd {};
+  std::mt19937 prng {rd()};
+  std::normal_distribution<> d {0, 1};
+
   for (int i=0; i<num_sims; i++) {
-    double gauss_bm = gaussian_box_muller();
-    S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
-    payoff_sum += std::max(S_cur - K, 0.0);
+    s_cur = s_adjust * exp(sqrt(v*v*T)*d(prng));
+    payoff_sum += std::max(s_cur - K, 0.0);
   }
 
   return (payoff_sum / static_cast<double>(num_sims)) * exp(-r*T);