Ankit3445/options_pricing
0
1#include "PricingEngine.hpp"2 3#include <cmath>4#include <limits>5#include <numbers>6#include <random>7#include <stdexcept>8#include <algorithm>9#include <numeric>10 11// ───── Black-Scholes Engine ──────────────────────────────────────────────12 13double BlackScholesEngine::cum_norm(double x) {14 return 0.5 * (1.0 + std::erf(x / std::sqrt(2.0)));15}16 17double BlackScholesEngine::norm_pdf(double x) {18 static constexpr double inv_sqrt_2pi = 1.0 / std::sqrt(2.0 * std::numbers::pi);19 return inv_sqrt_2pi * std::exp(-0.5 * x * x);20}21 22double BlackScholesEngine::bs_price(23 double S, double K, double T, double r, double sigma, OptionType type)24{25 if (sigma <= 0.0 || T <= 0.0) {26 double intrinsic = (type == OptionType::CALL)27 ? std::max(0.0, S - K * std::exp(-r * T))28 : std::max(0.0, K * std::exp(-r * T) - S);29 return intrinsic;30 }31 32 double d1 = (std::log(S / K) + (r + sigma * sigma * 0.5) * T)33 / (sigma * std::sqrt(T));34 double d2 = d1 - sigma * std::sqrt(T);35 36 if (type == OptionType::CALL) {37 return S * cum_norm(d1) - K * std::exp(-r * T) * cum_norm(d2);38 } else {39 return K * std::exp(-r * T) * cum_norm(-d2) - S * cum_norm(-d1);40 }41}42 43double BlackScholesEngine::price(44 double S, double K, double T, double r, double sigma, OptionType type)45{46 return bs_price(S, K, T, r, sigma, type);47}48 49Greeks BlackScholesEngine::bs_greeks(50 double S, double K, double T, double r, double sigma, OptionType type)51{52 Greeks g;53 if (sigma <= 0.0 || T <= 0.0) {54 return g;55 }56 57 double sqrt_T = std::sqrt(T);58 double d1 = (std::log(S / K) + (r + sigma * sigma * 0.5) * T)59 / (sigma * sqrt_T);60 double d2 = d1 - sigma * sqrt_T;61 62 double nd1 = cum_norm(d1);63 double nd2 = cum_norm(d2);64 double pdf_d1 = norm_pdf(d1);65 double disc = std::exp(-r * T);66 67 if (type == OptionType::CALL) {68 g.price = S * nd1 - K * disc * nd2;69 g.delta = nd1;70 g.theta = -S * pdf_d1 * sigma / (2.0 * sqrt_T)71 - r * K * disc * nd2;72 g.rho = K * T * disc * nd2;73 } else {74 g.price = K * disc * cum_norm(-d2) - S * cum_norm(-d1);75 g.delta = nd1 - 1.0;76 g.theta = -S * pdf_d1 * sigma / (2.0 * sqrt_T)77 + r * K * disc * cum_norm(-d2);78 g.rho = -K * T * disc * cum_norm(-d2);79 }80 81 g.gamma = pdf_d1 / (S * sigma * sqrt_T);82 g.vega = S * pdf_d1 * sqrt_T;83 84 return g;85}86 87Greeks BlackScholesEngine::greeks(88 double S, double K, double T, double r, double sigma, OptionType type)89{90 return bs_greeks(S, K, T, r, sigma, type);91}92 93double BlackScholesEngine::implied_vol_newton(94 double target, double S, double K, double T, double r, OptionType type)95{96 double sigma = 0.3;97 constexpr int max_iter = 80;98 constexpr double tol = 1e-10;99 100 for (int i = 0; i < max_iter; ++i) {101 double sqrt_T = std::sqrt(T);102 double d1 = (std::log(S / K) + (r + sigma * sigma * 0.5) * T)103 / (sigma * sqrt_T);104 double price = bs_price(S, K, T, r, sigma, type);105 double vega = S * norm_pdf(d1) * sqrt_T;106 double diff = price - target;107 108 if (std::abs(diff) < tol) {109 return sigma;110 }111 if (std::abs(vega) < 1e-12) {112 break;113 }114 115 sigma = sigma - diff / vega;116 if (sigma <= 0.0) sigma = 0.001;117 if (sigma > 10.0) sigma = 10.0;118 }119 120 return implied_vol_bisection(target, S, K, T, r, type);121}122 123double BlackScholesEngine::implied_vol_bisection(124 double target, double S, double K, double T, double r, OptionType type)125{126 double lo = 0.001;127 double hi = 5.0;128 constexpr int max_iter = 100;129 constexpr double tol = 1e-10;130 131 double price_lo = bs_price(S, K, T, r, lo, type);132 double price_hi = bs_price(S, K, T, r, hi, type);133 134 if ((target - price_lo) * (target - price_hi) > 0.0) {135 if (target < price_lo) return lo;136 if (target > price_hi) return hi;137 return std::numeric_limits<double>::quiet_NaN();138 }139 140 for (int i = 0; i < max_iter; ++i) {141 double mid = (lo + hi) * 0.5;142 double price_mid = bs_price(S, K, T, r, mid, type);143 double diff = price_mid - target;144 145 if (std::abs(diff) < tol) {146 return mid;147 }148 149 if ((price_mid - target) * (price_lo - target) < 0.0) {150 hi = mid;151 price_hi = price_mid;152 } else {153 lo = mid;154 price_lo = price_mid;155 }156 }157 158 return std::numeric_limits<double>::quiet_NaN();159}160 161double BlackScholesEngine::implied_volatility(162 double market_price, double S, double K, double T, double r, OptionType type)163{164 double iv = implied_vol_newton(market_price, S, K, T, r, type);165 if (std::isnan(iv) || iv <= 0.0) {166 return 0.25;167 }168 return iv;169}170 171Greeks BlackScholesEngine::calculate(const MarketTick& tick) {172 Greeks g = bs_greeks(173 tick.spotPrice,174 tick.strikePrice,175 tick.timeToExpiry,176 tick.riskFreeRate,177 0.25,178 tick.optionType179 );180 181 double iv = implied_vol_newton(182 tick.optionMarketPrice,183 tick.spotPrice,184 tick.strikePrice,185 tick.timeToExpiry,186 tick.riskFreeRate,187 tick.optionType188 );189 190 if (!std::isnan(iv)) {191 g = bs_greeks(192 tick.spotPrice,193 tick.strikePrice,194 tick.timeToExpiry,195 tick.riskFreeRate,196 iv,197 tick.optionType198 );199 }200 201 g.impliedVol = iv;202 g.valid = !std::isnan(iv) && iv > 0.0;203 204 return g;205}206 207// ───── Monte Carlo Engine ────────────────────────────────────────────────208 209MonteCarloEngine::MonteCarloEngine(unsigned int num_paths)210 : num_paths_(num_paths) {}211 212void MonteCarloEngine::generate_normals(std::vector<double>& z) const {213 static thread_local std::mt19937_64 rng(std::random_device{}());214 std::normal_distribution<double> norm(0.0, 1.0);215 for (size_t i = 0; i < z.size(); ++i) {216 z[i] = norm(rng);217 }218}219 220void MonteCarloEngine::generate_paths(221 std::vector<std::vector<double>>& paths, double S, double T, double r, double sigma) const222{223 static thread_local std::mt19937_64 rng(std::random_device{}());224 std::normal_distribution<double> norm(0.0, 1.0);225 226 double dt = T / paths[0].size();227 double drift = (r - 0.5 * sigma * sigma) * dt;228 double diffusion = sigma * std::sqrt(dt);229 230 for (auto& path : paths) {231 path[0] = S;232 for (size_t i = 1; i < path.size(); ++i) {233 double z = norm(rng);234 path[i] = path[i - 1] * std::exp(drift + diffusion * z);235 }236 }237}238 239// ── Standard European MC (antithetic) ──240 241double MonteCarloEngine::run_mc(242 double S, double K, double T, double r, double sigma,243 OptionType type, const std::vector<double>& z) const244{245 if (sigma <= 0.0 || T <= 0.0) {246 double intrinsic = (type == OptionType::CALL)247 ? std::max(0.0, S - K * std::exp(-r * T))248 : std::max(0.0, K * std::exp(-r * T) - S);249 return intrinsic;250 }251 252 double drift = (r - 0.5 * sigma * sigma) * T;253 double diffusion = sigma * std::sqrt(T);254 double discount = std::exp(-r * T);255 double half = static_cast<double>(z.size());256 double sum = 0.0;257 258 for (const double& zi : z) {259 double s_t1 = S * std::exp(drift + diffusion * zi);260 double s_t2 = S * std::exp(drift + diffusion * (-zi));261 262 double payoff1, payoff2;263 if (type == OptionType::CALL) {264 payoff1 = std::max(s_t1 - K, 0.0);265 payoff2 = std::max(s_t2 - K, 0.0);266 } else {267 payoff1 = std::max(K - s_t1, 0.0);268 payoff2 = std::max(K - s_t2, 0.0);269 }270 271 sum += discount * (payoff1 + payoff2);272 }273 274 return sum / (2.0 * half);275}276 277// ── Asian MC (arithmetic average) ──278 279double MonteCarloEngine::run_mc_asian(280 double S, double K, double T, double r, double sigma,281 OptionType type, unsigned int steps,282 const std::vector<std::vector<double>>& paths) const283{284 if (sigma <= 0.0 || T <= 0.0) {285 double intrinsic = (type == OptionType::CALL)286 ? std::max(0.0, S - K * std::exp(-r * T))287 : std::max(0.0, K * std::exp(-r * T) - S);288 return intrinsic;289 }290 291 double discount = std::exp(-r * T);292 double sum = 0.0;293 294 for (const auto& path : paths) {295 double avg = std::accumulate(path.begin(), path.end(), 0.0) / static_cast<double>(steps + 1);296 double payoff;297 if (type == OptionType::CALL) {298 payoff = std::max(avg - K, 0.0);299 } else {300 payoff = std::max(K - avg, 0.0);301 }302 sum += discount * payoff;303 }304 305 return sum / static_cast<double>(paths.size());306}307 308// ── Barrier MC (down-and-out / up-and-out) ──309 310double MonteCarloEngine::run_mc_barrier(311 double S, double K, double T, double r, double sigma,312 OptionType type, double barrier, bool down_and_out,313 const std::vector<double>& z) const314{315 if (sigma <= 0.0 || T <= 0.0) {316 double intrinsic = (type == OptionType::CALL)317 ? std::max(0.0, S - K * std::exp(-r * T))318 : std::max(0.0, K * std::exp(-r * T) - S);319 return intrinsic;320 }321 322 double drift = (r - 0.5 * sigma * sigma) * T;323 double diffusion = sigma * std::sqrt(T);324 double discount = std::exp(-r * T);325 double half = static_cast<double>(z.size());326 double sum = 0.0;327 328 // Simulate with 52 time steps to check barrier329 unsigned int steps = 52;330 double dt = T / steps;331 double drift_step = (r - 0.5 * sigma * sigma) * dt;332 double diff_step = sigma * std::sqrt(dt);333 334 for (const double& zi : z) {335 // Antithetic: +zi and -zi336 for (double sign : {1.0, -1.0}) {337 double s = S;338 bool knocked = false;339 for (unsigned int i = 0; i < steps; ++i) {340 double eps = zi * sign;341 // Sub-step: use fresh normals if more resolution needed342 s = s * std::exp(drift_step + diff_step * eps);343 if (down_and_out && s <= barrier) { knocked = true; break; }344 if (!down_and_out && s >= barrier) { knocked = true; break; }345 }346 if (!knocked) {347 double payoff = (type == OptionType::CALL)348 ? std::max(s - K, 0.0)349 : std::max(K - s, 0.0);350 sum += discount * payoff;351 }352 }353 }354 355 return sum / (2.0 * half);356}357 358// ── Lookback MC (floating strike) ──359 360double MonteCarloEngine::run_mc_lookback(361 double S, double K, double T, double r, double sigma,362 OptionType type, const std::vector<double>& z) const363{364 if (sigma <= 0.0 || T <= 0.0) {365 double intrinsic = (type == OptionType::CALL)366 ? std::max(0.0, S - K * std::exp(-r * T))367 : std::max(0.0, K * std::exp(-r * T) - S);368 return intrinsic;369 }370 371 unsigned int steps = 52;372 double dt = T / steps;373 double drift_step = (r - 0.5 * sigma * sigma) * dt;374 double diff_step = sigma * std::sqrt(dt);375 double discount = std::exp(-r * T);376 double half = static_cast<double>(z.size());377 double sum = 0.0;378 379 for (const double& zi : z) {380 for (double sign : {1.0, -1.0}) {381 double s = S;382 double s_max = S;383 double s_min = S;384 for (unsigned int i = 0; i < steps; ++i) {385 double eps = zi * sign;386 s = s * std::exp(drift_step + diff_step * eps);387 if (s > s_max) s_max = s;388 if (s < s_min) s_min = s;389 }390 double payoff;391 if (type == OptionType::CALL) {392 // Call on max: payoff = max(S_max - K, 0)393 payoff = std::max(s_max - K, 0.0);394 } else {395 // Put on min: payoff = max(K - S_min, 0)396 payoff = std::max(K - s_min, 0.0);397 }398 sum += discount * payoff;399 }400 }401 402 return sum / (2.0 * half);403}404 405// ── Greeks via finite differences (European) ──406 407Greeks MonteCarloEngine::fd_greeks(408 double S, double K, double T, double r, double sigma,409 OptionType type, const std::vector<double>& z) const410{411 Greeks g;412 413 double eps_s = S * 0.01;414 double eps_sigma = 0.001;415 double eps_T = 0.001;416 double eps_r = 0.0001;417 418 double v0 = run_mc(S, K, T, r, sigma, type, z);419 double v_sp = run_mc(S + eps_s, K, T, r, sigma, type, z);420 double v_sm = run_mc(S - eps_s, K, T, r, sigma, type, z);421 double v_vp = run_mc(S, K, T, r, sigma + eps_sigma, type, z);422 double v_vm = run_mc(S, K, T, r, sigma - eps_sigma, type, z);423 double v_tm = run_mc(S, K, T - eps_T, r, sigma, type, z);424 double v_rp = run_mc(S, K, T, r + eps_r, sigma, type, z);425 double v_rm = run_mc(S, K, T, r - eps_r, sigma, type, z);426 427 g.price = v0;428 g.delta = (v_sp - v_sm) / (2.0 * eps_s);429 g.gamma = (v_sp - 2.0 * v0 + v_sm) / (eps_s * eps_s);430 g.vega = (v_vp - v_vm) / (2.0 * eps_sigma);431 g.theta = (v_tm - v0) / eps_T;432 g.rho = (v_rp - v_rm) / (2.0 * eps_r);433 434 return g;435}436 437// ── Greeks via finite differences (exotic) ──438 439Greeks MonteCarloEngine::fd_greeks_exotic(440 double S, double K, double T, double r, double sigma,441 OptionType type, unsigned int style,442 double barrier, bool down_and_out,443 const std::vector<double>& z,444 const std::vector<std::vector<double>>& paths) const445{446 Greeks g;447 double eps_s = S * 0.01;448 double eps_sigma = 0.001;449 double eps_T = 0.001;450 double eps_r = 0.0001;451 452 auto price_fn = [&](double sp, double k, double t, double rt, double sg, const std::vector<double>& zn) -> double {453 (void)k;454 if (style == 1) { // Asian455 return run_mc_asian(sp, K, t, rt, sg, type, static_cast<unsigned int>(paths[0].size() - 1), paths);456 } else if (style == 2) { // Barrier457 return run_mc_barrier(sp, K, t, rt, sg, type, barrier, down_and_out, zn);458 } else { // Lookback459 return run_mc_lookback(sp, K, t, rt, sg, type, zn);460 }461 };462 463 // Need a modified path set for Asian - approximate by regenerating with same seed464 // For barrier and lookback, use the same z vector465 if (style == 1) {466 // Asian: recompute paths for each shift (coarse approximation)467 // For finite differences, use path regeneration468 g.price = price_fn(S, K, T, r, sigma, z);469 g.delta = 0.0; g.gamma = 0.0; g.vega = 0.0; g.theta = 0.0; g.rho = 0.0;470 } else {471 double v0 = price_fn(S, K, T, r, sigma, z);472 double v_sp = price_fn(S + eps_s, K, T, r, sigma, z);473 double v_sm = price_fn(S - eps_s, K, T, r, sigma, z);474 double v_vp = price_fn(S, K, T, r, sigma + eps_sigma, z);475 double v_vm = price_fn(S, K, T, r, sigma - eps_sigma, z);476 double v_tm = price_fn(S, K, T - eps_T, r, sigma, z);477 double v_rp = price_fn(S, K, T, r + eps_r, sigma, z);478 double v_rm = price_fn(S, K, T, r - eps_r, sigma, z);479 480 g.price = v0;481 g.delta = (v_sp - v_sm) / (2.0 * eps_s);482 g.gamma = (v_sp - 2.0 * v0 + v_sm) / (eps_s * eps_s);483 g.vega = (v_vp - v_vm) / (2.0 * eps_sigma);484 g.theta = (v_tm - v0) / eps_T;485 g.rho = (v_rp - v_rm) / (2.0 * eps_r);486 }487 488 return g;489}490 491// ── Public API: calculate European ──492 493Greeks MonteCarloEngine::calculate(const MarketTick& tick) {494 size_t half = num_paths_ / 2;495 std::vector<double> z(half);496 generate_normals(z);497 498 Greeks g;499 if (tick.optionStyle == OptionStyle::ASIAN) {500 unsigned int steps = 52;501 std::vector<std::vector<double>> paths(half * 2, std::vector<double>(steps + 1));502 generate_paths(paths, tick.spotPrice, tick.timeToExpiry, tick.riskFreeRate, 0.25);503 double price = run_mc_asian(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,504 tick.riskFreeRate, 0.25, tick.optionType, steps, paths);505 g.price = price;506 g.valid = true;507 } else if (tick.optionStyle == OptionStyle::BARRIER) {508 double price = run_mc_barrier(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,509 tick.riskFreeRate, 0.25, tick.optionType,510 tick.barrier, tick.barrier_down, z);511 Greeks fd = fd_greeks_exotic(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,512 tick.riskFreeRate, 0.25, tick.optionType, 2,513 tick.barrier, tick.barrier_down, z,514 std::vector<std::vector<double>>());515 g = fd;516 g.price = price;517 g.valid = true;518 } else if (tick.optionStyle == OptionStyle::LOOKBACK) {519 double price = run_mc_lookback(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,520 tick.riskFreeRate, 0.25, tick.optionType, z);521 Greeks fd = fd_greeks_exotic(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,522 tick.riskFreeRate, 0.25, tick.optionType, 3,523 0.0, false, z, std::vector<std::vector<double>>());524 g = fd;525 g.price = price;526 g.valid = true;527 } else {528 // Standard European529 double sigma_guess = 0.25;530 g = fd_greeks(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,531 tick.riskFreeRate, sigma_guess, tick.optionType, z);532 533 double iv = BlackScholesEngine::implied_volatility(534 tick.optionMarketPrice,535 tick.spotPrice, tick.strikePrice,536 tick.timeToExpiry, tick.riskFreeRate,537 tick.optionType);538 539 if (!std::isnan(iv) && iv > 0.0) {540 g = fd_greeks(tick.spotPrice, tick.strikePrice, tick.timeToExpiry,541 tick.riskFreeRate, iv, tick.optionType, z);542 }543 544 g.impliedVol = iv;545 g.valid = !std::isnan(iv) && iv > 0.0;546 }547 548 return g;549}550 551// ── Exotic option public APIs ──552 553Greeks MonteCarloEngine::calculate_asian(554 double S, double K, double T, double r, double sigma,555 OptionType type, unsigned int steps) const556{557 Greeks g;558 559 size_t half = num_paths_ / 2;560 // We need paths for Asian - generate (half * 2) paths561 std::vector<std::vector<double>> paths(half * 2, std::vector<double>(steps + 1));562 563 // Use generate_paths but we need non-const access... use local RNG564 static thread_local std::mt19937_64 rng(std::random_device{}());565 std::normal_distribution<double> norm(0.0, 1.0);566 567 double dt = T / steps;568 double drift = (r - 0.5 * sigma * sigma) * dt;569 double diffusion = sigma * std::sqrt(dt);570 571 for (auto& path : paths) {572 path[0] = S;573 for (size_t i = 1; i <= steps; ++i) {574 double z = norm(rng);575 path[i] = path[i - 1] * std::exp(drift + diffusion * z);576 }577 }578 579 double discount = std::exp(-r * T);580 double sum = 0.0, sum_sq = 0.0;581 582 for (const auto& path : paths) {583 double avg = std::accumulate(path.begin(), path.end(), 0.0) / static_cast<double>(steps + 1);584 double payoff = (type == OptionType::CALL)585 ? std::max(avg - K, 0.0)586 : std::max(K - avg, 0.0);587 double pv = discount * payoff;588 sum += pv;589 sum_sq += pv * pv;590 }591 592 double n = static_cast<double>(paths.size());593 g.price = sum / n;594 double variance = (sum_sq - sum * sum / n) / (n - 1.0);595 g.valid = true;596 597 return g;598}599 600Greeks MonteCarloEngine::calculate_barrier(601 double S, double K, double T, double r, double sigma,602 OptionType type, double barrier, bool down_and_out) const603{604 size_t half = num_paths_ / 2;605 std::vector<double> z(half);606 generate_normals(z);607 608 Greeks g;609 g.price = run_mc_barrier(S, K, T, r, sigma, type, barrier, down_and_out, z);610 g.valid = true;611 612 // Approximate delta via finite differences613 double eps = S * 0.01;614 double vp = run_mc_barrier(S + eps, K, T, r, sigma, type, barrier, down_and_out, z);615 double vm = run_mc_barrier(S - eps, K, T, r, sigma, type, barrier, down_and_out, z);616 g.delta = (vp - vm) / (2.0 * eps);617 618 return g;619}620 621Greeks MonteCarloEngine::calculate_lookback(622 double S, double K, double T, double r, double sigma,623 OptionType type) const624{625 size_t half = num_paths_ / 2;626 std::vector<double> z(half);627 generate_normals(z);628 629 Greeks g;630 g.price = run_mc_lookback(S, K, T, r, sigma, type, z);631 g.valid = true;632 633 // Approximate delta634 double eps = S * 0.01;635 double vp = run_mc_lookback(S + eps, K, T, r, sigma, type, z);636 double vm = run_mc_lookback(S - eps, K, T, r, sigma, type, z);637 g.delta = (vp - vm) / (2.0 * eps);638 639 return g;640}641 642// ── Convergence Analysis ──643 644std::vector<ConvergencePoint> MonteCarloEngine::convergence(645 double S, double K, double T, double r, double sigma,646 OptionType type, unsigned int max_paths) const647{648 std::vector<unsigned int> path_counts = {1000, 2000, 5000, 10000, 20000, 50000, 100000};649 // Filter counts <= max_paths650 std::vector<unsigned int> use_counts;651 for (auto c : path_counts) {652 if (c <= max_paths) use_counts.push_back(c);653 }654 if (use_counts.empty()) use_counts.push_back(max_paths);655 656 double bs = BlackScholesEngine::price(S, K, T, r, sigma, type);657 658 std::vector<ConvergencePoint> results;659 static thread_local std::mt19937_64 rng(std::random_device{}());660 std::normal_distribution<double> norm(0.0, 1.0);661 662 for (unsigned int paths : use_counts) {663 size_t half = paths / 2;664 std::vector<double> z(half);665 for (size_t i = 0; i < half; ++i) {666 z[i] = norm(rng);667 }668 669 double mc = run_mc(S, K, T, r, sigma, type, z);670 double error = mc - bs;671 672 // Standard error: std(payoffs) / sqrt(N)673 double drift = (r - 0.5 * sigma * sigma) * T;674 double diffusion = sigma * std::sqrt(T);675 double discount = std::exp(-r * T);676 double sum_payoffs = 0.0, sum_sq = 0.0;677 678 for (const double& zi : z) {679 double s_t1 = S * std::exp(drift + diffusion * zi);680 double s_t2 = S * std::exp(drift + diffusion * (-zi));681 double p1 = (type == OptionType::CALL) ? std::max(s_t1 - K, 0.0) : std::max(K - s_t1, 0.0);682 double p2 = (type == OptionType::CALL) ? std::max(s_t2 - K, 0.0) : std::max(K - s_t2, 0.0);683 double pv1 = discount * p1;684 double pv2 = discount * p2;685 sum_payoffs += pv1 + pv2;686 sum_sq += pv1 * pv1 + pv2 * pv2;687 }688 689 double N = static_cast<double>(2 * half);690 double variance = (sum_sq - sum_payoffs * sum_payoffs / N) / (N - 1.0);691 double std_err = std::sqrt(variance / N);692 693 ConvergencePoint cp;694 cp.paths = paths;695 cp.bs_price = bs;696 cp.mc_price = mc;697 cp.error = error;698 cp.std_error = std_err;699 results.push_back(cp);700 }701 702 return results;703}704