|
| 1 | +#pragma once |
| 2 | + |
| 3 | +#include "lfmc/numerical_scheme.hpp" |
| 4 | +#include "lfmc/payoff.hpp" |
| 5 | +#include "lfmc/stochastic_process.hpp" |
| 6 | +#include "lfmc/types.hpp" |
| 7 | + |
| 8 | +#include <cassert> |
| 9 | +#include <cstdint> |
| 10 | +#include <functional> |
| 11 | +#include <memory> |
| 12 | +#include <random> |
| 13 | +#include <string> |
| 14 | +#include <thread> |
| 15 | +#include <utility> |
| 16 | +#include <vector> |
| 17 | + |
| 18 | +namespace lfmc { |
| 19 | + |
| 20 | +// A strategy sampler: given n samples and a seed, returns n iid samples of the |
| 21 | +// (possibly variance-reduced) payoff estimator. The seed guarantees RNG independence |
| 22 | +// between strategies and between exploration/exploitation phases. |
| 23 | +using SamplerFn = std::function<std::vector<double>(size_t n_samples, uint64_t seed)>; |
| 24 | + |
| 25 | +// Per-strategy statistics collected during the exploration phase |
| 26 | +struct StrategyStats { |
| 27 | + std::string name; |
| 28 | + double mean; |
| 29 | + double sample_variance; // unbiased (n-1 denominator) |
| 30 | + double wall_time_ms; |
| 31 | + size_t n_samples; |
| 32 | +}; |
| 33 | + |
| 34 | +// Full result of one ASVR run |
| 35 | +struct ASVRResult { |
| 36 | + double estimate; |
| 37 | + double estimated_variance; // Var(mu_ASVR), estimated from exploration data |
| 38 | + double estimated_stderr; // sqrt(estimated_variance) |
| 39 | + |
| 40 | + std::vector<StrategyStats> exploration_stats; // one entry per strategy |
| 41 | + std::vector<double> precision_weights; // w_k* used to allocate exploitation |
| 42 | + std::vector<size_t> exploitation_counts; // actual n_k samples per strategy |
| 43 | + |
| 44 | + // Diagnostics — useful for paper tables |
| 45 | + double plain_mc_variance_estimate; // sigma_0^2 / N (first strategy treated as plain MC) |
| 46 | + double variance_reduction_ratio; // plain_mc_variance / estimated_variance |
| 47 | + size_t n_exploration; |
| 48 | + size_t n_exploitation; |
| 49 | + size_t total_samples; |
| 50 | +}; |
| 51 | + |
| 52 | +struct ASVRConfig { |
| 53 | + // Fraction of total_samples used for exploration (split equally across K strategies) |
| 54 | + double exploration_fraction = 0.1; |
| 55 | + // Total thread budget; exploration uses up to K threads, exploitation uses the rest |
| 56 | + size_t n_threads = std::thread::hardware_concurrency(); |
| 57 | + // Minimum exploration samples per strategy regardless of exploration_fraction. |
| 58 | + // Needs to be large enough for a stable variance estimate (~100+ samples). |
| 59 | + size_t min_exploration_per_strategy = 100; |
| 60 | +}; |
| 61 | + |
| 62 | +// ─── ASVR Algorithm ────────────────────────────────────────────────────────── |
| 63 | +// |
| 64 | +// Phase 1 (Exploration): K strategies run in parallel, each generating |
| 65 | +// n_explore_each = max(min_per_strategy, alpha*N/K) samples. |
| 66 | +// Each strategy uses seed make_seed(k, 0) for full RNG independence. |
| 67 | +// |
| 68 | +// Weight computation: w_k* = (1/sigma_k^2) / sum_j(1/sigma_j^2) |
| 69 | +// (inverse-variance / "precision" weighting) |
| 70 | +// |
| 71 | +// Phase 2 (Exploitation): strategy k receives n_k = floor(w_k* * n_exploit) samples, |
| 72 | +// each run using seed make_seed(k, 1) — independent of the exploration phase. |
| 73 | +// Strategies with n_k > 0 run in parallel. |
| 74 | +// |
| 75 | +// Combination (unbiased by independence): |
| 76 | +// mu_explore = (1/K) * sum_k mu_k^explore (equal-weight pooled) |
| 77 | +// mu_exploit = sum_k w_k* * mu_k^exploit (precision-weighted) |
| 78 | +// mu_ASVR = alpha * mu_explore + (1-alpha) * mu_exploit |
| 79 | +// |
| 80 | +// The key unbiasedness guarantee: w_k* depends only on exploration samples; |
| 81 | +// exploitation samples are generated with fresh RNG seeds and are therefore |
| 82 | +// independent of the weight selection event. |
| 83 | + |
| 84 | +class AdaptiveVarianceReduction { |
| 85 | + public: |
| 86 | + static ASVRResult run(std::vector<std::pair<std::string, SamplerFn>> strategies, |
| 87 | + size_t total_samples, ASVRConfig config = {}); |
| 88 | + |
| 89 | + static StrategyStats run_strategy_batch(const std::string& name, const SamplerFn& sampler, |
| 90 | + size_t n_samples, uint64_t seed); |
| 91 | + |
| 92 | + private: |
| 93 | + static std::vector<double> compute_precision_weights(const std::vector<StrategyStats>& stats); |
| 94 | +}; |
| 95 | + |
| 96 | +// ─── Internal utilities ─────────────────────────────────────────────────────── |
| 97 | + |
| 98 | +namespace detail { |
| 99 | + |
| 100 | +// Splitmix64-based seed derivation: strategy k, phase p → unique uint64_t seed. |
| 101 | +// Two calls with the same (k, p) give the same seed; different (k, p) pairs are |
| 102 | +// statistically independent because they hit different streams of splitmix64. |
| 103 | +inline uint64_t make_seed(size_t k, size_t phase) noexcept { |
| 104 | + uint64_t x = (static_cast<uint64_t>(k) * 0x9e3779b97f4a7c15ULL) |
| 105 | + ^ (static_cast<uint64_t>(phase) * 0x6c62272e07bb0142ULL) |
| 106 | + ^ 0xdeadbeefcafe0000ULL; |
| 107 | + x ^= x >> 30; |
| 108 | + x *= 0xbf58476d1ce4e5b9ULL; |
| 109 | + x ^= x >> 27; |
| 110 | + x *= 0x94d049bb133111ebULL; |
| 111 | + x ^= x >> 31; |
| 112 | + return x; |
| 113 | +} |
| 114 | + |
| 115 | +// Generate `n` draws from N(0,1) |
| 116 | +inline std::vector<double> gen_normals(size_t n, std::mt19937_64& rng) { |
| 117 | + std::normal_distribution<double> dist{0.0, 1.0}; |
| 118 | + std::vector<double> v(n); |
| 119 | + for (auto& x : v) |
| 120 | + x = dist(rng); |
| 121 | + return v; |
| 122 | +} |
| 123 | + |
| 124 | +// Correct single-path generator (avoids the pre-allocate + push_back double-length bug |
| 125 | +// present in the existing PathGenerator) |
| 126 | +template <StochasticProcess SP, NumericalScheme<SP> NS> |
| 127 | +Path generate_path(const SP& process, const NS& scheme, const Normals& normals, size_t steps, |
| 128 | + double T) { |
| 129 | + const double dt = T / static_cast<double>(steps); |
| 130 | + Path path; |
| 131 | + path.reserve(steps + 1); |
| 132 | + double x = process.initial(); |
| 133 | + path.push_back(x); |
| 134 | + double t = 0.0; |
| 135 | + for (size_t i = 0; i < steps; ++i) { |
| 136 | + x = scheme.step(process, t, x, dt, normals[i]); |
| 137 | + path.push_back(x); |
| 138 | + t += dt; |
| 139 | + } |
| 140 | + return path; |
| 141 | +} |
| 142 | + |
| 143 | +// Unbiased sample mean and variance (Welford would be more numerically stable for large n, |
| 144 | +// but two-pass is fine for the batch sizes used here) |
| 145 | +inline std::pair<double, double> mean_variance(const std::vector<double>& samples) { |
| 146 | + assert(!samples.empty()); |
| 147 | + const double n = static_cast<double>(samples.size()); |
| 148 | + double sum = 0.0; |
| 149 | + for (double x : samples) |
| 150 | + sum += x; |
| 151 | + const double mean = sum / n; |
| 152 | + double sq = 0.0; |
| 153 | + for (double x : samples) { |
| 154 | + double d = x - mean; |
| 155 | + sq += d * d; |
| 156 | + } |
| 157 | + // n-1 denominator for unbiased variance; guard against n=1 |
| 158 | + const double variance = (samples.size() > 1) ? sq / (n - 1.0) : 0.0; |
| 159 | + return {mean, variance}; |
| 160 | +} |
| 161 | + |
| 162 | +} // namespace detail |
| 163 | + |
| 164 | +// ─── Strategy factory functions ────────────────────────────────────────────── |
| 165 | +// |
| 166 | +// Each factory captures the pricing parameters and returns a SamplerFn. |
| 167 | +// The SamplerFn is called with (n_samples, seed) and returns exactly n_samples |
| 168 | +// iid draws from the (possibly variance-reduced) estimator of E[payoff]. |
| 169 | + |
| 170 | +// Plain pseudo-random Monte Carlo — baseline strategy |
| 171 | +template <StochasticProcess SP, NumericalScheme<SP> NS> |
| 172 | +SamplerFn make_plain_mc_sampler( SP process, NS scheme, std::shared_ptr<Payoff> payoff, |
| 173 | + size_t steps, double T) { |
| 174 | + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector<double> { |
| 175 | + std::mt19937_64 rng{seed}; |
| 176 | + std::vector<double> samples; |
| 177 | + samples.reserve(n); |
| 178 | + for (size_t i = 0; i < n; ++i) { |
| 179 | + auto normals = detail::gen_normals(steps, rng); |
| 180 | + auto path = detail::generate_path(process, scheme, normals, steps, T); |
| 181 | + auto result = payoff->generate_payoffs({path}); |
| 182 | + if (result && !result->empty() && !(*result)[0].empty()) |
| 183 | + samples.push_back((*result)[0][0]); |
| 184 | + } |
| 185 | + return samples; |
| 186 | + }; |
| 187 | +} |
| 188 | + |
| 189 | +// Antithetic variates — each sample is (payoff(Z) + payoff(-Z)) / 2 |
| 190 | +// Returns n samples, each consuming one pair of paths. Variance is reduced |
| 191 | +// because the two paths are negatively correlated for monotone payoffs. |
| 192 | +template <StochasticProcess SP, NumericalScheme<SP> NS> |
| 193 | +SamplerFn make_antithetic_sampler( SP process, NS scheme, std::shared_ptr<Payoff> payoff, |
| 194 | + size_t steps, double T) { |
| 195 | + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector<double> { |
| 196 | + std::mt19937_64 rng{seed}; |
| 197 | + std::vector<double> samples; |
| 198 | + samples.reserve(n); |
| 199 | + for (size_t i = 0; i < n; ++i) { |
| 200 | + auto z = detail::gen_normals(steps, rng); |
| 201 | + Normals neg_z(steps); |
| 202 | + for (size_t j = 0; j < steps; ++j) |
| 203 | + neg_z[j] = -z[j]; |
| 204 | + |
| 205 | + auto path_pos = detail::generate_path(process, scheme, z, steps, T); |
| 206 | + auto path_neg = detail::generate_path(process, scheme, neg_z, steps, T); |
| 207 | + |
| 208 | + auto r_pos = payoff->generate_payoffs({path_pos}); |
| 209 | + auto r_neg = payoff->generate_payoffs({path_neg}); |
| 210 | + |
| 211 | + if (r_pos && r_neg && !r_pos->empty() && !r_neg->empty() && !(*r_pos)[0].empty() && |
| 212 | + !(*r_neg)[0].empty()) |
| 213 | + samples.push_back(0.5 * ((*r_pos)[0][0] + (*r_neg)[0][0])); |
| 214 | + } |
| 215 | + return samples; |
| 216 | + }; |
| 217 | +} |
| 218 | + |
| 219 | +// Control variates with in-batch OLS beta estimation. |
| 220 | +// target_payoff: the option price we're estimating |
| 221 | +// control_payoff: a payoff with known expectation `control_mean` |
| 222 | +// (e.g. EuropeanCall(0.0) gives S_T; E[S_T] = S0 * exp(mu*T)) |
| 223 | +// Returns: X_i - beta_hat * (Y_i - E[Y]) for each i |
| 224 | +template <StochasticProcess SP, NumericalScheme<SP> NS> |
| 225 | +SamplerFn make_control_variate_sampler( SP process, NS scheme, std::shared_ptr<Payoff> target, |
| 226 | + std::shared_ptr<Payoff> control, double control_mean, |
| 227 | + size_t steps, double T) { |
| 228 | + return [process, scheme, target, control, control_mean, steps, |
| 229 | + T](size_t n, uint64_t seed) -> std::vector<double> { |
| 230 | + std::mt19937_64 rng{seed}; |
| 231 | + std::vector<double> xs, ys; |
| 232 | + xs.reserve(n); |
| 233 | + ys.reserve(n); |
| 234 | + |
| 235 | + for (size_t i = 0; i < n; ++i) { |
| 236 | + auto z = detail::gen_normals(steps, rng); |
| 237 | + auto path = detail::generate_path(process, scheme, z, steps, T); |
| 238 | + auto rx = target->generate_payoffs({path}); |
| 239 | + auto ry = control->generate_payoffs({path}); |
| 240 | + if (rx && ry && !rx->empty() && !ry->empty() && !(*rx)[0].empty() && |
| 241 | + !(*ry)[0].empty()) { |
| 242 | + xs.push_back((*rx)[0][0]); |
| 243 | + ys.push_back((*ry)[0][0]); |
| 244 | + } |
| 245 | + } |
| 246 | + |
| 247 | + // OLS beta: minimises Var(X - beta*(Y - E[Y])) |
| 248 | + const size_t m = xs.size(); |
| 249 | + double mean_x = 0.0, mean_y = 0.0; |
| 250 | + for (size_t i = 0; i < m; ++i) { |
| 251 | + mean_x += xs[i]; |
| 252 | + mean_y += ys[i]; |
| 253 | + } |
| 254 | + mean_x /= static_cast<double>(m); |
| 255 | + mean_y /= static_cast<double>(m); |
| 256 | + |
| 257 | + double cov_xy = 0.0, var_y = 0.0; |
| 258 | + for (size_t i = 0; i < m; ++i) { |
| 259 | + cov_xy += (xs[i] - mean_x) * (ys[i] - mean_y); |
| 260 | + var_y += (ys[i] - mean_y) * (ys[i] - mean_y); |
| 261 | + } |
| 262 | + const double beta = (var_y > 0.0) ? cov_xy / var_y : 0.0; |
| 263 | + |
| 264 | + std::vector<double> samples(m); |
| 265 | + for (size_t i = 0; i < m; ++i) |
| 266 | + samples[i] = xs[i] - beta * (ys[i] - control_mean); |
| 267 | + return samples; |
| 268 | + }; |
| 269 | +} |
| 270 | + |
| 271 | +// Antithetic + control variates combined. |
| 272 | +// Each raw sample is the antithetic average; OLS beta is estimated from those averages. |
| 273 | +template <StochasticProcess SP, NumericalScheme<SP> NS> |
| 274 | +SamplerFn make_antithetic_cv_sampler( SP process, NS scheme, std::shared_ptr<Payoff> target, |
| 275 | + std::shared_ptr<Payoff> control, double control_mean, |
| 276 | + size_t steps, double T) { |
| 277 | + return [process, scheme, target, control, control_mean, steps, |
| 278 | + T](size_t n, uint64_t seed) -> std::vector<double> { |
| 279 | + std::mt19937_64 rng{seed}; |
| 280 | + std::vector<double> xs, ys; |
| 281 | + xs.reserve(n); |
| 282 | + ys.reserve(n); |
| 283 | + |
| 284 | + for (size_t i = 0; i < n; ++i) { |
| 285 | + auto z = detail::gen_normals(steps, rng); |
| 286 | + Normals neg_z(steps); |
| 287 | + for (size_t j = 0; j < steps; ++j) |
| 288 | + neg_z[j] = -z[j]; |
| 289 | + |
| 290 | + auto path_p = detail::generate_path(process, scheme, z, steps, T); |
| 291 | + auto path_n = detail::generate_path(process, scheme, neg_z, steps, T); |
| 292 | + |
| 293 | + auto xp = target->generate_payoffs({path_p}); |
| 294 | + auto xn = target->generate_payoffs({path_n}); |
| 295 | + auto yp = control->generate_payoffs({path_p}); |
| 296 | + auto yn = control->generate_payoffs({path_n}); |
| 297 | + |
| 298 | + if (xp && xn && yp && yn && !xp->empty() && !xn->empty() && !yp->empty() && |
| 299 | + !yn->empty() && !(*xp)[0].empty() && !(*xn)[0].empty() && !(*yp)[0].empty() && |
| 300 | + !(*yn)[0].empty()) { |
| 301 | + xs.push_back(0.5 * ((*xp)[0][0] + (*xn)[0][0])); |
| 302 | + ys.push_back(0.5 * ((*yp)[0][0] + (*yn)[0][0])); |
| 303 | + } |
| 304 | + } |
| 305 | + |
| 306 | + const size_t m = xs.size(); |
| 307 | + double mean_x = 0.0, mean_y = 0.0; |
| 308 | + for (size_t i = 0; i < m; ++i) { |
| 309 | + mean_x += xs[i]; |
| 310 | + mean_y += ys[i]; |
| 311 | + } |
| 312 | + mean_x /= static_cast<double>(m); |
| 313 | + mean_y /= static_cast<double>(m); |
| 314 | + |
| 315 | + double cov_xy = 0.0, var_y = 0.0; |
| 316 | + for (size_t i = 0; i < m; ++i) { |
| 317 | + cov_xy += (xs[i] - mean_x) * (ys[i] - mean_y); |
| 318 | + var_y += (ys[i] - mean_y) * (ys[i] - mean_y); |
| 319 | + } |
| 320 | + const double beta = (var_y > 0.0) ? cov_xy / var_y : 0.0; |
| 321 | + |
| 322 | + std::vector<double> samples(m); |
| 323 | + for (size_t i = 0; i < m; ++i) |
| 324 | + samples[i] = xs[i] - beta * (ys[i] - control_mean); |
| 325 | + return samples; |
| 326 | + }; |
| 327 | +} |
| 328 | + |
| 329 | +} // namespace lfmc |
0 commit comments