From 5b15818952d3e88965512dbb2401f6cefa678fa5 Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Fri, 20 Feb 2026 00:42:07 -0500 Subject: [PATCH 01/12] temp archive --- include/lfmc/{ => archive}/Estimator.hpp | 0 include/lfmc/{ => archive}/Manager.hpp | 0 include/lfmc/{ => archive}/NumericalScheme.hpp | 0 include/lfmc/{ => archive}/PathGenerator.hpp | 0 include/lfmc/{ => archive}/Payoff.hpp | 0 include/lfmc/{ => archive}/RandomGenerator.hpp | 0 include/lfmc/{ => archive}/StochasticProcess.hpp | 0 include/lfmc/{ => archive}/VarianceReductionStrategy.hpp | 0 include/lfmc/{ => archive}/types.hpp | 0 9 files changed, 0 insertions(+), 0 deletions(-) rename include/lfmc/{ => archive}/Estimator.hpp (100%) rename include/lfmc/{ => archive}/Manager.hpp (100%) rename include/lfmc/{ => archive}/NumericalScheme.hpp (100%) rename include/lfmc/{ => archive}/PathGenerator.hpp (100%) rename include/lfmc/{ => archive}/Payoff.hpp (100%) rename include/lfmc/{ => archive}/RandomGenerator.hpp (100%) rename include/lfmc/{ => archive}/StochasticProcess.hpp (100%) rename include/lfmc/{ => archive}/VarianceReductionStrategy.hpp (100%) rename include/lfmc/{ => archive}/types.hpp (100%) diff --git a/include/lfmc/Estimator.hpp b/include/lfmc/archive/Estimator.hpp similarity index 100% rename from include/lfmc/Estimator.hpp rename to include/lfmc/archive/Estimator.hpp diff --git a/include/lfmc/Manager.hpp b/include/lfmc/archive/Manager.hpp similarity index 100% rename from include/lfmc/Manager.hpp rename to include/lfmc/archive/Manager.hpp diff --git a/include/lfmc/NumericalScheme.hpp b/include/lfmc/archive/NumericalScheme.hpp similarity index 100% rename from include/lfmc/NumericalScheme.hpp rename to include/lfmc/archive/NumericalScheme.hpp diff --git a/include/lfmc/PathGenerator.hpp b/include/lfmc/archive/PathGenerator.hpp similarity index 100% rename from include/lfmc/PathGenerator.hpp rename to include/lfmc/archive/PathGenerator.hpp diff --git a/include/lfmc/Payoff.hpp b/include/lfmc/archive/Payoff.hpp similarity index 100% rename from include/lfmc/Payoff.hpp rename to include/lfmc/archive/Payoff.hpp diff --git a/include/lfmc/RandomGenerator.hpp b/include/lfmc/archive/RandomGenerator.hpp similarity index 100% rename from include/lfmc/RandomGenerator.hpp rename to include/lfmc/archive/RandomGenerator.hpp diff --git a/include/lfmc/StochasticProcess.hpp b/include/lfmc/archive/StochasticProcess.hpp similarity index 100% rename from include/lfmc/StochasticProcess.hpp rename to include/lfmc/archive/StochasticProcess.hpp diff --git a/include/lfmc/VarianceReductionStrategy.hpp b/include/lfmc/archive/VarianceReductionStrategy.hpp similarity index 100% rename from include/lfmc/VarianceReductionStrategy.hpp rename to include/lfmc/archive/VarianceReductionStrategy.hpp diff --git a/include/lfmc/types.hpp b/include/lfmc/archive/types.hpp similarity index 100% rename from include/lfmc/types.hpp rename to include/lfmc/archive/types.hpp From 9e3ba8254b63143c7d0259438026c380543c5ea4 Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Fri, 20 Feb 2026 01:33:06 -0500 Subject: [PATCH 02/12] rng engine and path engine --- include/lfmc/Engine.hpp | 1 + include/lfmc/Estimator.hpp | 1 + include/lfmc/NumericalScheme.hpp | 42 +++++++++++ include/lfmc/PathGenerator.hpp | 41 +++++++++++ include/lfmc/Payoff.hpp | 1 + include/lfmc/Pipeline.hpp | 1 + include/lfmc/RandomSource.hpp | 28 +++++++ include/lfmc/StochasticProcess.hpp | 20 +++++ include/lfmc/archive/NumericalScheme.hpp | 86 ---------------------- include/lfmc/archive/RandomGenerator.hpp | 2 +- include/lfmc/archive/StochasticProcess.hpp | 55 -------------- include/lfmc/archive/types.hpp | 22 ------ include/lfmc/types.hpp | 10 +++ src/CMakeLists.txt | 1 + src/EulerMaruyama.cpp | 17 +++++ src/GeometricBrownianMotion.cpp | 25 +++++++ src/RandomSource.cpp | 21 ++++++ 17 files changed, 210 insertions(+), 164 deletions(-) create mode 100644 include/lfmc/Engine.hpp create mode 100644 include/lfmc/Estimator.hpp create mode 100644 include/lfmc/NumericalScheme.hpp create mode 100644 include/lfmc/PathGenerator.hpp create mode 100644 include/lfmc/Payoff.hpp create mode 100644 include/lfmc/Pipeline.hpp create mode 100644 include/lfmc/RandomSource.hpp create mode 100644 include/lfmc/StochasticProcess.hpp delete mode 100644 include/lfmc/archive/NumericalScheme.hpp delete mode 100644 include/lfmc/archive/StochasticProcess.hpp delete mode 100644 include/lfmc/archive/types.hpp create mode 100644 include/lfmc/types.hpp create mode 100644 src/EulerMaruyama.cpp create mode 100644 src/GeometricBrownianMotion.cpp create mode 100644 src/RandomSource.cpp diff --git a/include/lfmc/Engine.hpp b/include/lfmc/Engine.hpp new file mode 100644 index 0000000..6f70f09 --- /dev/null +++ b/include/lfmc/Engine.hpp @@ -0,0 +1 @@ +#pragma once diff --git a/include/lfmc/Estimator.hpp b/include/lfmc/Estimator.hpp new file mode 100644 index 0000000..6f70f09 --- /dev/null +++ b/include/lfmc/Estimator.hpp @@ -0,0 +1 @@ +#pragma once diff --git a/include/lfmc/NumericalScheme.hpp b/include/lfmc/NumericalScheme.hpp new file mode 100644 index 0000000..40509b5 --- /dev/null +++ b/include/lfmc/NumericalScheme.hpp @@ -0,0 +1,42 @@ +#pragma once + +#include + +namespace lfmc { + +template +concept NumericalScheme = + requires(S const& s, P const& p, double t, double x, double dt, double z) { + { s.step(p, t, x, dt, z) } -> std::convertible_to; + }; + +/** + * @brief Exact simulation for Geometric Brownian Motion. + * + * Uses the closed-form solution: + * X_T = X_0 * exp((mu - 0.5*sigma^2)*T + sigma*sqrt(T)*Z) + * + * This is faster and more accurate than Euler-Maruyama for GBM. + */ +// class GBMExact { +// public: +// explicit GBMExact(GeometricBrownianMotion gbm) : gbm_(gbm) {} +// +// /** +// * @brief Simulate terminal value using exact solution. +// * @param x0 Initial value. +// * @param T Time to maturity. +// * @param z Standard normal random variable. +// * @return Terminal value X_T. +// */ +// double simulate_terminal(double x0, double T, double z) const noexcept { +// double drift_adjusted = (gbm_.mu - 0.5 * gbm_.sigma * gbm_.sigma) * T; +// double diffusion_term = gbm_.sigma * std::sqrt(T) * z; +// return x0 * std::exp(drift_adjusted + diffusion_term); +// } +// +// private: +// GeometricBrownianMotion gbm_; +// }; + +} // namespace lfmc diff --git a/include/lfmc/PathGenerator.hpp b/include/lfmc/PathGenerator.hpp new file mode 100644 index 0000000..256e55f --- /dev/null +++ b/include/lfmc/PathGenerator.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include "lfmc/NumericalScheme.hpp" +#include "lfmc/StochasticProcess.hpp" +#include "lfmc/types.hpp" + +namespace lfmc { + +template + requires NumericalScheme +class PathGenerator { + private: + Process process_; + Scheme scheme_; + double T_; + size_t steps_; + + public: + PathGenerator(Process process, Scheme scheme, double T, size_t steps) + : process_(std::move(process)), scheme_(std::move(scheme)), T_(T), steps_(steps) {} + + Path generate(const Normals& normals) const { + const double dt = T_ / steps_; + + Path path(steps_ + 1); + + double t = 0.0; + double x = process_.initial(); + path.push_back(x); + + for (size_t i = 0; i < steps_; ++i) { + x = scheme_.step(process_, t, x, dt, normals[i]); + path.push_back(x); + t += dt; + } + + return path; + } +}; + +} // namespace lfmc diff --git a/include/lfmc/Payoff.hpp b/include/lfmc/Payoff.hpp new file mode 100644 index 0000000..6f70f09 --- /dev/null +++ b/include/lfmc/Payoff.hpp @@ -0,0 +1 @@ +#pragma once diff --git a/include/lfmc/Pipeline.hpp b/include/lfmc/Pipeline.hpp new file mode 100644 index 0000000..6f70f09 --- /dev/null +++ b/include/lfmc/Pipeline.hpp @@ -0,0 +1 @@ +#pragma once diff --git a/include/lfmc/RandomSource.hpp b/include/lfmc/RandomSource.hpp new file mode 100644 index 0000000..46a58fd --- /dev/null +++ b/include/lfmc/RandomSource.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include "lfmc/types.hpp" + +#include + +namespace lfmc { + +class RandomSource { + public: + virtual ~RandomSource() = default; + virtual Normals generate(size_t n) = 0; +}; + +class PseudoRandomSource : public RandomSource { + private: + std::mt19937 rng_; + std::normal_distribution dist_; + + public: + PseudoRandomSource(unsigned seed = std::random_device{}()); + + Normals generate(size_t n) override; + + void seed(unsigned seed); +}; + +} // namespace lfmc diff --git a/include/lfmc/StochasticProcess.hpp b/include/lfmc/StochasticProcess.hpp new file mode 100644 index 0000000..a33e3ff --- /dev/null +++ b/include/lfmc/StochasticProcess.hpp @@ -0,0 +1,20 @@ +#pragma once + +#include + +/** + * @file StochasticProcess.hpp + * @brief Defines the StochasticProcess concept for stochastic differential equations (SDEs) and + * provides some implementations. + */ + +namespace lfmc { + +template +concept StochasticProcess = requires(P const& p, double t, double x) { + { p.initial() } -> std::convertible_to; + { p.drift(t, x) } -> std::convertible_to; + { p.diffusion(t, x) } -> std::convertible_to; +}; + +} // namespace lfmc diff --git a/include/lfmc/archive/NumericalScheme.hpp b/include/lfmc/archive/NumericalScheme.hpp deleted file mode 100644 index ac9614b..0000000 --- a/include/lfmc/archive/NumericalScheme.hpp +++ /dev/null @@ -1,86 +0,0 @@ -#pragma once - -#include "StochasticProcess.hpp" - -#include -#include - -/** - * @file NumericalScheme.hpp - * @brief Defines the NumericalScheme concept for numerical methods solving SDEs and - * provides some implementations. - */ - -namespace lfmc { - -/** - * @brief Concept for numerical schemes solving SDEs. - * - * A NumericalScheme must implement a step function that computes the next state - * given the current state, time step, and a standard normal random variable. - * - * @tparam S Numerical scheme type. - * @tparam P Stochastic process type. - * - * Requires: - * - S must have a method: - * double step(const P& process, double x_current, double dt, double z) const noexcept; - * where: - * - process: Stochastic process defining drift and diffusion. - * - x_current: Current state. - * - dt: Time step size. - * - z: Standard normal random variable N(0,1). - * The method returns the next state X_{t+dt}. - */ -template -concept NumericalScheme = - StochasticProcess

&& requires(S const& s, P const& p, double x, double dt, double z) { - { s.step(p, x, dt, z) } -> std::same_as; - }; - -template struct EulerMaruyama { - /** - * @brief Compute the next state using Euler-Maruyama. - * @param Stochastic process defining drift and diffusion. - * @param x_current Current state. - * @param dt Time step size. - * @param z Standard normal random variable N(0,1). - * @return Next state X_{t+dt}. - */ - double step(const P& process, double x_current, double dt, double z) const noexcept { - double drift = process.drift(x_current); - double diffusion = process.diffusion(x_current); - return x_current + drift * dt + diffusion * std::sqrt(dt) * z; - } -}; - -/** - * @brief Exact simulation for Geometric Brownian Motion. - * - * Uses the closed-form solution: - * X_T = X_0 * exp((mu - 0.5*sigma^2)*T + sigma*sqrt(T)*Z) - * - * This is faster and more accurate than Euler-Maruyama for GBM. - */ -// class GBMExact { -// public: -// explicit GBMExact(GeometricBrownianMotion gbm) : gbm_(gbm) {} -// -// /** -// * @brief Simulate terminal value using exact solution. -// * @param x0 Initial value. -// * @param T Time to maturity. -// * @param z Standard normal random variable. -// * @return Terminal value X_T. -// */ -// double simulate_terminal(double x0, double T, double z) const noexcept { -// double drift_adjusted = (gbm_.mu - 0.5 * gbm_.sigma * gbm_.sigma) * T; -// double diffusion_term = gbm_.sigma * std::sqrt(T) * z; -// return x0 * std::exp(drift_adjusted + diffusion_term); -// } -// -// private: -// GeometricBrownianMotion gbm_; -// }; - -} // namespace lfmc diff --git a/include/lfmc/archive/RandomGenerator.hpp b/include/lfmc/archive/RandomGenerator.hpp index eae17d9..5a20f25 100644 --- a/include/lfmc/archive/RandomGenerator.hpp +++ b/include/lfmc/archive/RandomGenerator.hpp @@ -1,6 +1,6 @@ #pragma once -#include "types.hpp" +#include "lfmc/types.hpp" #include diff --git a/include/lfmc/archive/StochasticProcess.hpp b/include/lfmc/archive/StochasticProcess.hpp deleted file mode 100644 index 8c0acb1..0000000 --- a/include/lfmc/archive/StochasticProcess.hpp +++ /dev/null @@ -1,55 +0,0 @@ -#pragma once - -#include -#include - -/** - * @file StochasticProcess.hpp - * @brief Defines the StochasticProcess concept for stochastic differential equations (SDEs) and - * provides some implementations. - */ - -namespace lfmc { - -template -concept StochasticProcess = requires(P const& p, double x) { - { p.drift(x) } -> std::same_as; - { p.diffusion(x) } -> std::same_as; -}; - -struct GeometricBrownianMotion { - double mu; - double sigma; - - /** - * @brief Constructor to initialize GBM parameters. - * @param driftCoefficient The drift coefficient (mu). - * @param diffusionCoefficient The diffusion coefficient (sigma/volatility). - */ - GeometricBrownianMotion(double driftCoefficient, double diffusionCoefficient) - : mu(driftCoefficient), sigma(diffusionCoefficient) { - if (sigma < 0.0) { - throw std::invalid_argument("Diffusion coefficient (sigma) must be non-negative"); - } - } - - /** - * @brief Compute drift term at state x. - * @param x The current state variable. - * @return The drift term mu * x. - */ - double drift(double x) const noexcept { - return mu * x; - } - - /** - * @brief Compute diffusion term at state x. - * @param x The current state variable. - * @return The diffusion term sigma * x. - */ - double diffusion(double x) const noexcept { - return sigma * x; - } -}; - -} // namespace lfmc diff --git a/include/lfmc/archive/types.hpp b/include/lfmc/archive/types.hpp deleted file mode 100644 index 2aba43b..0000000 --- a/include/lfmc/archive/types.hpp +++ /dev/null @@ -1,22 +0,0 @@ -#pragma once - -#include - -namespace lfmc { - -struct State { - double initialValue; - double timeToMaturity; - size_t steps; -}; - -using Path = std::vector; -using Normals = std::vector; - -struct ManagerConfig { - // TODO temp - size_t numNoVarianceReductionSimulations; - size_t numAntitheticVariatesSimulations; -}; - -} // namespace lfmc diff --git a/include/lfmc/types.hpp b/include/lfmc/types.hpp new file mode 100644 index 0000000..df4c691 --- /dev/null +++ b/include/lfmc/types.hpp @@ -0,0 +1,10 @@ +#pragma once + +#include + +namespace lfmc { + +using Path = std::vector; +using Normals = std::vector; + +} // namespace lfmc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ea86e7a..c62486d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,7 @@ # Source files for the lfmc library set(LFMC_SOURCES timing.cpp + RandomSource.hpp ) # Create library diff --git a/src/EulerMaruyama.cpp b/src/EulerMaruyama.cpp new file mode 100644 index 0000000..1efba95 --- /dev/null +++ b/src/EulerMaruyama.cpp @@ -0,0 +1,17 @@ +#include "lfmc/StochasticProcess.hpp" + +#include + +namespace lfmc { + +class EulerMaruyama { + public: + template + double step(P const& process, double t, double x, double dt, double z) const noexcept { + double drift = process.drift(t, x); + double diffusion = process.diffusion(t, x); + return x + drift * dt + diffusion * std::sqrt(dt) * z; + } +}; + +} // namespace lfmc diff --git a/src/GeometricBrownianMotion.cpp b/src/GeometricBrownianMotion.cpp new file mode 100644 index 0000000..ec575fe --- /dev/null +++ b/src/GeometricBrownianMotion.cpp @@ -0,0 +1,25 @@ +namespace lfmc { + +class GeometricBrownianMotion { + private: + double mu_; + double sigma_; + double x0_; + + public: + GeometricBrownianMotion(double mu, double sigma, double x0) : mu_(mu), sigma_(sigma), x0_(x0) {} + + double initial() const noexcept { + return x0_; + } + + double drift(double, double x) const noexcept { + return mu_ * x; + } + + double diffusion(double, double x) const noexcept { + return sigma_ * x; + } +}; + +} // namespace lfmc diff --git a/src/RandomSource.cpp b/src/RandomSource.cpp new file mode 100644 index 0000000..6f021e4 --- /dev/null +++ b/src/RandomSource.cpp @@ -0,0 +1,21 @@ +#include "lfmc/RandomSource.hpp" + +#include "lfmc/types.hpp" + +namespace lfmc { + +PseudoRandomSource::PseudoRandomSource(unsigned seed) : rng_(seed), dist_(0.0, 1.0) {} + +Normals PseudoRandomSource::generate(size_t n) { + Normals normals(n); + for (size_t i = 0; i < n; ++i) { + normals[i] = dist_(rng_); + } + return normals; +} + +void PseudoRandomSource::seed(unsigned seed) { + rng_.seed(seed); +} + +} // namespace lfmc From 426f56e0e1714858596d3367a499199777e9e75b Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Mon, 23 Feb 2026 23:21:43 -0500 Subject: [PATCH 03/12] feat: estimator and payoff outline --- include/lfmc/Estimator.hpp | 13 +++++++++++++ include/lfmc/Payoff.hpp | 17 +++++++++++++++++ include/lfmc/Pipeline.hpp | 5 +++++ 3 files changed, 35 insertions(+) diff --git a/include/lfmc/Estimator.hpp b/include/lfmc/Estimator.hpp index 6f70f09..5ebb939 100644 --- a/include/lfmc/Estimator.hpp +++ b/include/lfmc/Estimator.hpp @@ -1 +1,14 @@ #pragma once + +namespace lfmc { + +class Estimator { + public: + virtual void add(double x) = 0; + // virtual bool converged() const = 0; + // virtual Result result() const = 0; + // virtual void merge(Estimator const& other) = 0; + virtual ~Estimator() = default; +}; + +} // namespace lfmc diff --git a/include/lfmc/Payoff.hpp b/include/lfmc/Payoff.hpp index 6f70f09..f9ebb72 100644 --- a/include/lfmc/Payoff.hpp +++ b/include/lfmc/Payoff.hpp @@ -1 +1,18 @@ #pragma once + +#include +#include + +namespace lfmc { + +template +concept TerminalPayoff = requires(P const& p, double terminal) { + { p(terminal) } -> std::convertible_to; +}; + +template +concept PathPayoff = requires(P const& p, std::vector const& path) { + { p(path) } -> std::convertible_to; +}; + +} // namespace lfmc diff --git a/include/lfmc/Pipeline.hpp b/include/lfmc/Pipeline.hpp index 6f70f09..9e42e2c 100644 --- a/include/lfmc/Pipeline.hpp +++ b/include/lfmc/Pipeline.hpp @@ -1 +1,6 @@ #pragma once + +#include "lfmc/PathGenerator.hpp" +#include "lfmc/Payoff.hpp" +#include "lfmc/Pipeline.hpp" +#include "lfmc/RandomSource.hpp" From b80146752943731f3d23be1c1117105fdaab501b Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Tue, 24 Feb 2026 12:49:04 -0500 Subject: [PATCH 04/12] feat: payoffs and estimator --- include/lfmc/Estimator.hpp | 5 ++++- src/MonteCarloEstimator.cpp | 36 ++++++++++++++++++++++++++++++++++++ src/Payoff.cpp | 21 +++++++++++++++++++++ 3 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 src/MonteCarloEstimator.cpp create mode 100644 src/Payoff.cpp diff --git a/include/lfmc/Estimator.hpp b/include/lfmc/Estimator.hpp index 5ebb939..9a67508 100644 --- a/include/lfmc/Estimator.hpp +++ b/include/lfmc/Estimator.hpp @@ -1,12 +1,15 @@ #pragma once +#include +#include + namespace lfmc { class Estimator { public: virtual void add(double x) = 0; // virtual bool converged() const = 0; - // virtual Result result() const = 0; + virtual std::expected result() const = 0; // virtual void merge(Estimator const& other) = 0; virtual ~Estimator() = default; }; diff --git a/src/MonteCarloEstimator.cpp b/src/MonteCarloEstimator.cpp new file mode 100644 index 0000000..234d9f9 --- /dev/null +++ b/src/MonteCarloEstimator.cpp @@ -0,0 +1,36 @@ +#include "lfmc/Estimator.hpp" + +#include + +namespace lfmc { + +class MonteCarloEstimator : public Estimator { + private: + double sum = 0.0; + std::size_t count = 0; + + public: + void add(double x) override { + sum += x; + ++count; + } + + // bool converged() const override { + // return count >= 10000; // Placeholder convergence criterion + // } + + std::expected result() const override { + if (count == 0) { + return std::unexpected{"No samples added"}; + } + return {sum / static_cast(count)}; + } + + // void merge(Estimator const& other) override { + // auto const& mcOther = dynamic_cast(other); + // sum += mcOther.sum; + // count += mcOther.count; + // } +}; + +} // namespace lfmc diff --git a/src/Payoff.cpp b/src/Payoff.cpp new file mode 100644 index 0000000..8838519 --- /dev/null +++ b/src/Payoff.cpp @@ -0,0 +1,21 @@ +#include "lfmc/Payoff.hpp" + +namespace lfmc { + +struct EuropeanCall { + double strike; + + double operator()(double terminal) const { + return std::max(terminal - strike, 0.0); + } +}; + +struct EuropeanPut { + double strike; + + double operator()(double terminal) const { + return std::max(strike - terminal, 0.0); + } +}; + +} // namespace lfmc From b22c6d98d2a993300892c35a05d6268a68c2c714 Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Tue, 24 Feb 2026 14:21:13 -0500 Subject: [PATCH 05/12] refactor: renamed everything, finished pipeline, and converted each step to return vectors of results --- include/lfmc/Estimator.hpp | 17 ------- include/lfmc/PathGenerator.hpp | 41 ----------------- include/lfmc/Payoff.hpp | 18 -------- include/lfmc/Pipeline.hpp | 6 --- include/lfmc/{Engine.hpp => engine.hpp} | 0 include/lfmc/estimator.hpp | 31 +++++++++++++ ...mericalScheme.hpp => numerical_scheme.hpp} | 0 include/lfmc/path_generator.hpp | 44 +++++++++++++++++++ include/lfmc/payoff.hpp | 33 ++++++++++++++ include/lfmc/pipeline.hpp | 43 ++++++++++++++++++ .../{RandomSource.hpp => random_source.hpp} | 5 ++- ...sticProcess.hpp => stochastic_process.hpp} | 14 ++++++ include/lfmc/types.hpp | 1 + src/CMakeLists.txt | 7 ++- src/GeometricBrownianMotion.cpp | 25 ----------- src/MonteCarloEstimator.cpp | 36 --------------- src/Payoff.cpp | 21 --------- src/RandomSource.cpp | 21 --------- src/{EulerMaruyama.cpp => euler_maruyama.cpp} | 5 +-- src/european_payoffs.cpp | 29 ++++++++++++ src/geometric_brownian_motion.cpp | 20 +++++++++ src/monte_carlo_estimator.cpp | 40 +++++++++++++++++ src/pseudo_random_source.cpp | 24 ++++++++++ 23 files changed, 289 insertions(+), 192 deletions(-) delete mode 100644 include/lfmc/Estimator.hpp delete mode 100644 include/lfmc/PathGenerator.hpp delete mode 100644 include/lfmc/Payoff.hpp delete mode 100644 include/lfmc/Pipeline.hpp rename include/lfmc/{Engine.hpp => engine.hpp} (100%) create mode 100644 include/lfmc/estimator.hpp rename include/lfmc/{NumericalScheme.hpp => numerical_scheme.hpp} (100%) create mode 100644 include/lfmc/path_generator.hpp create mode 100644 include/lfmc/payoff.hpp create mode 100644 include/lfmc/pipeline.hpp rename include/lfmc/{RandomSource.hpp => random_source.hpp} (70%) rename include/lfmc/{StochasticProcess.hpp => stochastic_process.hpp} (61%) delete mode 100644 src/GeometricBrownianMotion.cpp delete mode 100644 src/MonteCarloEstimator.cpp delete mode 100644 src/Payoff.cpp delete mode 100644 src/RandomSource.cpp rename src/{EulerMaruyama.cpp => euler_maruyama.cpp} (73%) create mode 100644 src/european_payoffs.cpp create mode 100644 src/geometric_brownian_motion.cpp create mode 100644 src/monte_carlo_estimator.cpp create mode 100644 src/pseudo_random_source.cpp diff --git a/include/lfmc/Estimator.hpp b/include/lfmc/Estimator.hpp deleted file mode 100644 index 9a67508..0000000 --- a/include/lfmc/Estimator.hpp +++ /dev/null @@ -1,17 +0,0 @@ -#pragma once - -#include -#include - -namespace lfmc { - -class Estimator { - public: - virtual void add(double x) = 0; - // virtual bool converged() const = 0; - virtual std::expected result() const = 0; - // virtual void merge(Estimator const& other) = 0; - virtual ~Estimator() = default; -}; - -} // namespace lfmc diff --git a/include/lfmc/PathGenerator.hpp b/include/lfmc/PathGenerator.hpp deleted file mode 100644 index 256e55f..0000000 --- a/include/lfmc/PathGenerator.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#pragma once - -#include "lfmc/NumericalScheme.hpp" -#include "lfmc/StochasticProcess.hpp" -#include "lfmc/types.hpp" - -namespace lfmc { - -template - requires NumericalScheme -class PathGenerator { - private: - Process process_; - Scheme scheme_; - double T_; - size_t steps_; - - public: - PathGenerator(Process process, Scheme scheme, double T, size_t steps) - : process_(std::move(process)), scheme_(std::move(scheme)), T_(T), steps_(steps) {} - - Path generate(const Normals& normals) const { - const double dt = T_ / steps_; - - Path path(steps_ + 1); - - double t = 0.0; - double x = process_.initial(); - path.push_back(x); - - for (size_t i = 0; i < steps_; ++i) { - x = scheme_.step(process_, t, x, dt, normals[i]); - path.push_back(x); - t += dt; - } - - return path; - } -}; - -} // namespace lfmc diff --git a/include/lfmc/Payoff.hpp b/include/lfmc/Payoff.hpp deleted file mode 100644 index f9ebb72..0000000 --- a/include/lfmc/Payoff.hpp +++ /dev/null @@ -1,18 +0,0 @@ -#pragma once - -#include -#include - -namespace lfmc { - -template -concept TerminalPayoff = requires(P const& p, double terminal) { - { p(terminal) } -> std::convertible_to; -}; - -template -concept PathPayoff = requires(P const& p, std::vector const& path) { - { p(path) } -> std::convertible_to; -}; - -} // namespace lfmc diff --git a/include/lfmc/Pipeline.hpp b/include/lfmc/Pipeline.hpp deleted file mode 100644 index 9e42e2c..0000000 --- a/include/lfmc/Pipeline.hpp +++ /dev/null @@ -1,6 +0,0 @@ -#pragma once - -#include "lfmc/PathGenerator.hpp" -#include "lfmc/Payoff.hpp" -#include "lfmc/Pipeline.hpp" -#include "lfmc/RandomSource.hpp" diff --git a/include/lfmc/Engine.hpp b/include/lfmc/engine.hpp similarity index 100% rename from include/lfmc/Engine.hpp rename to include/lfmc/engine.hpp diff --git a/include/lfmc/estimator.hpp b/include/lfmc/estimator.hpp new file mode 100644 index 0000000..33128dc --- /dev/null +++ b/include/lfmc/estimator.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include "lfmc/types.hpp" + +#include +#include + +namespace lfmc { + +class Estimator { + public: + virtual std::expected add_payoffs(const Payoffs& payoffs) = 0; + virtual bool converged() const = 0; + virtual std::expected result() const = 0; + // virtual void merge(Estimator const& other) = 0; + virtual ~Estimator() = default; +}; + +class MonteCarloEstimator : public Estimator { + private: + double sum = 0.0; + std::size_t count = 0; + + public: + std::expected add_payoffs(const Payoffs& payoffs) override; + bool converged() const override; + std::expected result() const override; + // void merge(Estimator const& other) override; +}; + +} // namespace lfmc diff --git a/include/lfmc/NumericalScheme.hpp b/include/lfmc/numerical_scheme.hpp similarity index 100% rename from include/lfmc/NumericalScheme.hpp rename to include/lfmc/numerical_scheme.hpp diff --git a/include/lfmc/path_generator.hpp b/include/lfmc/path_generator.hpp new file mode 100644 index 0000000..b09918f --- /dev/null +++ b/include/lfmc/path_generator.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/types.hpp" + +namespace lfmc { + +template + requires NumericalScheme +class PathGenerator { + private: + Process process_; + Scheme scheme_; + + public: + PathGenerator(Process process, Scheme scheme) + : process_(std::move(process)), scheme_(std::move(scheme)) {} + + // TODO move to cpp file + std::vector generate_paths(const std::vector& normals, size_t steps, + double T) const { + const double dt = T / static_cast(steps); + + std::vector paths; + for (const auto& norm : normals) { + Path path(steps + 1); + + double t = 0.0; + double x = process_.initial(); + path.push_back(x); + + for (size_t i = 0; i < steps; ++i) { + x = scheme_.step(process_, t, x, dt, norm[i]); + path.push_back(x); + t += dt; + } + } + + return paths; + } +}; + +} // namespace lfmc diff --git a/include/lfmc/payoff.hpp b/include/lfmc/payoff.hpp new file mode 100644 index 0000000..d4627d5 --- /dev/null +++ b/include/lfmc/payoff.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include "lfmc/types.hpp" + +namespace lfmc { + +class Payoff { + public: + virtual Payoffs generate_payoffs(const std::vector& paths) const = 0; + virtual ~Payoff() = default; +}; + +class EuropeanCall : public Payoff { + private: + double strike_; + + public: + explicit EuropeanCall(double strike); + + Payoffs generate_payoffs(const std::vector& paths) const override; +}; + +class EuropeanPut : public Payoff { + private: + double strike_; + + public: + explicit EuropeanPut(double strike); + + Payoffs generate_payoffs(const std::vector& paths) const override; +}; + +} // namespace lfmc diff --git a/include/lfmc/pipeline.hpp b/include/lfmc/pipeline.hpp new file mode 100644 index 0000000..949e50e --- /dev/null +++ b/include/lfmc/pipeline.hpp @@ -0,0 +1,43 @@ +#pragma once + +#include "lfmc/estimator.hpp" +#include "lfmc/path_generator.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/random_source.hpp" + +#include +#include + +namespace lfmc { + +template NS> class Pipeline { + private: + std::unique_ptr random_source_; + std::unique_ptr> path_generator_; + std::unique_ptr payoff_; + std::unique_ptr estimator_; + + public: + Pipeline(std::unique_ptr random_source, + std::unique_ptr> path_generator, std::unique_ptr payoff, + std::unique_ptr estimator) + : random_source_(std::move(random_source)), path_generator_(std::move(path_generator)), + payoff_(std::move(payoff)), estimator_(std::move(estimator)) {} + + std::expected run(size_t steps, double T) { + while (!estimator_->converged()) { + auto normals = random_source_->generate_normals(steps); + auto paths = path_generator_->generate_paths(normals, steps, T); + auto payoffs = payoff_->generate_payoffs(paths); + auto result = estimator_->add_payoffs(payoffs); + + if (!result) { + return std::unexpected(result.error()); + } + } + + return estimator_->result(); + } +}; + +} // namespace lfmc diff --git a/include/lfmc/RandomSource.hpp b/include/lfmc/random_source.hpp similarity index 70% rename from include/lfmc/RandomSource.hpp rename to include/lfmc/random_source.hpp index 46a58fd..2b7fe2f 100644 --- a/include/lfmc/RandomSource.hpp +++ b/include/lfmc/random_source.hpp @@ -3,13 +3,14 @@ #include "lfmc/types.hpp" #include +#include namespace lfmc { class RandomSource { public: virtual ~RandomSource() = default; - virtual Normals generate(size_t n) = 0; + virtual std::vector generate_normals(size_t steps, size_t n = 1) = 0; }; class PseudoRandomSource : public RandomSource { @@ -20,7 +21,7 @@ class PseudoRandomSource : public RandomSource { public: PseudoRandomSource(unsigned seed = std::random_device{}()); - Normals generate(size_t n) override; + std::vector generate_normals(size_t steps, size_t) override; void seed(unsigned seed); }; diff --git a/include/lfmc/StochasticProcess.hpp b/include/lfmc/stochastic_process.hpp similarity index 61% rename from include/lfmc/StochasticProcess.hpp rename to include/lfmc/stochastic_process.hpp index a33e3ff..1d6271e 100644 --- a/include/lfmc/StochasticProcess.hpp +++ b/include/lfmc/stochastic_process.hpp @@ -17,4 +17,18 @@ concept StochasticProcess = requires(P const& p, double t, double x) { { p.diffusion(t, x) } -> std::convertible_to; }; +class GeometricBrownianMotion { + private: + double mu_; + double sigma_; + double x0_; + + public: + GeometricBrownianMotion(double mu, double sigma, double x0); + + double initial() const noexcept; + double drift(double, double x) const noexcept; + double diffusion(double, double x) const noexcept; +}; + } // namespace lfmc diff --git a/include/lfmc/types.hpp b/include/lfmc/types.hpp index df4c691..ebac5e6 100644 --- a/include/lfmc/types.hpp +++ b/include/lfmc/types.hpp @@ -6,5 +6,6 @@ namespace lfmc { using Path = std::vector; using Normals = std::vector; +using Payoffs = std::vector; } // namespace lfmc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c62486d..f38176f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,7 +1,10 @@ # Source files for the lfmc library set(LFMC_SOURCES - timing.cpp - RandomSource.hpp + euler_maruyama.cpp + european_payoffs.cpp + geometric_brownian_motion.cpp + monte_carlo_estimator.cpp + pseudo_random_source.cpp ) # Create library diff --git a/src/GeometricBrownianMotion.cpp b/src/GeometricBrownianMotion.cpp deleted file mode 100644 index ec575fe..0000000 --- a/src/GeometricBrownianMotion.cpp +++ /dev/null @@ -1,25 +0,0 @@ -namespace lfmc { - -class GeometricBrownianMotion { - private: - double mu_; - double sigma_; - double x0_; - - public: - GeometricBrownianMotion(double mu, double sigma, double x0) : mu_(mu), sigma_(sigma), x0_(x0) {} - - double initial() const noexcept { - return x0_; - } - - double drift(double, double x) const noexcept { - return mu_ * x; - } - - double diffusion(double, double x) const noexcept { - return sigma_ * x; - } -}; - -} // namespace lfmc diff --git a/src/MonteCarloEstimator.cpp b/src/MonteCarloEstimator.cpp deleted file mode 100644 index 234d9f9..0000000 --- a/src/MonteCarloEstimator.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include "lfmc/Estimator.hpp" - -#include - -namespace lfmc { - -class MonteCarloEstimator : public Estimator { - private: - double sum = 0.0; - std::size_t count = 0; - - public: - void add(double x) override { - sum += x; - ++count; - } - - // bool converged() const override { - // return count >= 10000; // Placeholder convergence criterion - // } - - std::expected result() const override { - if (count == 0) { - return std::unexpected{"No samples added"}; - } - return {sum / static_cast(count)}; - } - - // void merge(Estimator const& other) override { - // auto const& mcOther = dynamic_cast(other); - // sum += mcOther.sum; - // count += mcOther.count; - // } -}; - -} // namespace lfmc diff --git a/src/Payoff.cpp b/src/Payoff.cpp deleted file mode 100644 index 8838519..0000000 --- a/src/Payoff.cpp +++ /dev/null @@ -1,21 +0,0 @@ -#include "lfmc/Payoff.hpp" - -namespace lfmc { - -struct EuropeanCall { - double strike; - - double operator()(double terminal) const { - return std::max(terminal - strike, 0.0); - } -}; - -struct EuropeanPut { - double strike; - - double operator()(double terminal) const { - return std::max(strike - terminal, 0.0); - } -}; - -} // namespace lfmc diff --git a/src/RandomSource.cpp b/src/RandomSource.cpp deleted file mode 100644 index 6f021e4..0000000 --- a/src/RandomSource.cpp +++ /dev/null @@ -1,21 +0,0 @@ -#include "lfmc/RandomSource.hpp" - -#include "lfmc/types.hpp" - -namespace lfmc { - -PseudoRandomSource::PseudoRandomSource(unsigned seed) : rng_(seed), dist_(0.0, 1.0) {} - -Normals PseudoRandomSource::generate(size_t n) { - Normals normals(n); - for (size_t i = 0; i < n; ++i) { - normals[i] = dist_(rng_); - } - return normals; -} - -void PseudoRandomSource::seed(unsigned seed) { - rng_.seed(seed); -} - -} // namespace lfmc diff --git a/src/EulerMaruyama.cpp b/src/euler_maruyama.cpp similarity index 73% rename from src/EulerMaruyama.cpp rename to src/euler_maruyama.cpp index 1efba95..2048ae3 100644 --- a/src/EulerMaruyama.cpp +++ b/src/euler_maruyama.cpp @@ -1,12 +1,11 @@ -#include "lfmc/StochasticProcess.hpp" +#include "lfmc/stochastic_process.hpp" #include namespace lfmc { -class EulerMaruyama { +template class EulerMaruyama { public: - template double step(P const& process, double t, double x, double dt, double z) const noexcept { double drift = process.drift(t, x); double diffusion = process.diffusion(t, x); diff --git a/src/european_payoffs.cpp b/src/european_payoffs.cpp new file mode 100644 index 0000000..cc2c1c3 --- /dev/null +++ b/src/european_payoffs.cpp @@ -0,0 +1,29 @@ +#include "lfmc/payoff.hpp" + +#include + +namespace lfmc { + +EuropeanCall::EuropeanCall(double strike) : strike_(strike) {} + +Payoffs EuropeanCall::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + for (const auto& path : paths) { + double final_price = path.back(); + payoffs.push_back(std::max(final_price - strike_, 0.0)); + } + return payoffs; +} + +EuropeanPut::EuropeanPut(double strike) : strike_(strike) {} + +Payoffs EuropeanPut::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + for (const auto& path : paths) { + double final_price = path.back(); + payoffs.push_back(std::max(strike_ - final_price, 0.0)); + } + return payoffs; +} + +} // namespace lfmc diff --git a/src/geometric_brownian_motion.cpp b/src/geometric_brownian_motion.cpp new file mode 100644 index 0000000..98b196b --- /dev/null +++ b/src/geometric_brownian_motion.cpp @@ -0,0 +1,20 @@ +#include "lfmc/stochastic_process.hpp" + +namespace lfmc { + +GeometricBrownianMotion::GeometricBrownianMotion(double mu, double sigma, double x0) + : mu_(mu), sigma_(sigma), x0_(x0) {} + +double GeometricBrownianMotion::initial() const noexcept { + return x0_; +} + +double GeometricBrownianMotion::drift(double, double x) const noexcept { + return mu_ * x; +} + +double GeometricBrownianMotion::diffusion(double, double x) const noexcept { + return sigma_ * x; +} + +} // namespace lfmc diff --git a/src/monte_carlo_estimator.cpp b/src/monte_carlo_estimator.cpp new file mode 100644 index 0000000..092820f --- /dev/null +++ b/src/monte_carlo_estimator.cpp @@ -0,0 +1,40 @@ +#include "lfmc/estimator.hpp" +#include "lfmc/types.hpp" + +namespace lfmc { + +std::expected MonteCarloEstimator::add_payoffs(const Payoffs& payoffs) { + if (payoffs.empty()) { + return std::unexpected("No payoffs provided to add to the estimator."); + } + + for (double payoff : payoffs) { + sum += payoff; + ++count; + } + + return {}; +} + +bool MonteCarloEstimator::converged() const { + return count >= 10000; // Placeholder convergence criterion +} + +std::expected MonteCarloEstimator::result() const { + if (count <= 0) { + return std::unexpected("No samples added to the estimator."); + } + // if (!self.converged()) { + // return std::unexpected("Estimator has not yet converged."); + // } + + return {sum / static_cast(count)}; +} + +// void MonteCarloEstimator::merge(Estimator const& other) { +// auto const& mcOther = dynamic_cast(other); +// sum += mcOther.sum; +// count += mcOther.count; +// } + +} // namespace lfmc diff --git a/src/pseudo_random_source.cpp b/src/pseudo_random_source.cpp new file mode 100644 index 0000000..64ddcce --- /dev/null +++ b/src/pseudo_random_source.cpp @@ -0,0 +1,24 @@ +#include "lfmc/random_source.hpp" +#include "lfmc/types.hpp" + +namespace lfmc { + +PseudoRandomSource::PseudoRandomSource(unsigned seed) : rng_(seed), dist_(0.0, 1.0) {} + +std::vector PseudoRandomSource::generate_normals(size_t steps, size_t samples) { + std::vector result(samples, Normals(steps)); + for (size_t i = 0; i < samples; ++i) { + Normals normals(steps); + for (size_t j = 0; j < steps; ++j) { + normals[j] = dist_(rng_); + } + result[i] = std::move(normals); + } + return result; +} + +void PseudoRandomSource::seed(unsigned seed) { + rng_.seed(seed); +} + +} // namespace lfmc From eefa7a710e051d447b6b0adbbfdef7bab8eebc3d Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Tue, 3 Mar 2026 01:05:44 -0500 Subject: [PATCH 06/12] feat: finished refactor and added control variates --- include/lfmc/estimator.hpp | 27 +++++++- include/lfmc/numerical_scheme.hpp | 12 ++++ include/lfmc/path_generator.hpp | 1 + include/lfmc/payoff.hpp | 29 ++++++++- include/lfmc/pipeline.hpp | 16 ++++- include/lfmc/random_source.hpp | 14 ++++ include/lfmc/stochastic_process.hpp | 3 + src/CMakeLists.txt | 5 +- src/antithetic_random_source.cpp | 25 ++++++++ src/control_variate_estimator.cpp | 53 +++++++++++++++ src/control_variate_payoffs.cpp | 38 +++++++++++ src/euler_maruyama.cpp | 16 ----- src/european_payoffs.cpp | 12 ++-- src/geometric_brownian_motion.cpp | 8 +++ src/monte_carlo_estimator.cpp | 15 +++-- tests/CMakeLists.txt | 2 +- tests/test_manager.cpp | 99 ----------------------------- tests/test_pipeline.cpp | 88 +++++++++++++++++++++++++ tests/test_process.cpp | 66 ++++++++++--------- tests/test_scheme.cpp | 74 +++++++++++---------- 20 files changed, 406 insertions(+), 197 deletions(-) create mode 100644 src/antithetic_random_source.cpp create mode 100644 src/control_variate_estimator.cpp create mode 100644 src/control_variate_payoffs.cpp delete mode 100644 src/euler_maruyama.cpp delete mode 100644 tests/test_manager.cpp create mode 100644 tests/test_pipeline.cpp diff --git a/include/lfmc/estimator.hpp b/include/lfmc/estimator.hpp index 33128dc..cf056f9 100644 --- a/include/lfmc/estimator.hpp +++ b/include/lfmc/estimator.hpp @@ -4,12 +4,13 @@ #include #include +#include namespace lfmc { class Estimator { public: - virtual std::expected add_payoffs(const Payoffs& payoffs) = 0; + virtual std::expected add_payoffs(const std::vector& payoffs) = 0; virtual bool converged() const = 0; virtual std::expected result() const = 0; // virtual void merge(Estimator const& other) = 0; @@ -22,7 +23,29 @@ class MonteCarloEstimator : public Estimator { std::size_t count = 0; public: - std::expected add_payoffs(const Payoffs& payoffs) override; + std::expected add_payoffs(const std::vector& payoffs) override; + bool converged() const override; + std::expected result() const override; + // void merge(Estimator const& other) override; +}; + +// TODO if control is analytically known, can put it at payoff level +class ControlVariateEstimator : public Estimator { + private: + std::size_t count = 0; + + double sum_x = 0.0; // Sum of original payoffs + double sum_y = 0.0; // Sum of control variate values + double sum_xy = 0.0; // Sum of products of payoffs and control variate + double sum_yy = 0.0; // Sum of squares of control variate values + + double control_expectation_ = 0.0; // Expected value of control variate (known analytically) + + public: + explicit ControlVariateEstimator(double control_expectation) + : control_expectation_(control_expectation) {} + + std::expected add_payoffs(const std::vector& payoffs) override; bool converged() const override; std::expected result() const override; // void merge(Estimator const& other) override; diff --git a/include/lfmc/numerical_scheme.hpp b/include/lfmc/numerical_scheme.hpp index 40509b5..5302502 100644 --- a/include/lfmc/numerical_scheme.hpp +++ b/include/lfmc/numerical_scheme.hpp @@ -1,5 +1,8 @@ #pragma once +#include "lfmc/stochastic_process.hpp" + +#include #include namespace lfmc { @@ -10,6 +13,15 @@ concept NumericalScheme = { s.step(p, t, x, dt, z) } -> std::convertible_to; }; +template class EulerMaruyama { + public: + double step(P const& process, double t, double x, double dt, double z) const noexcept { + double drift = process.drift(t, x); + double diffusion = process.diffusion(t, x); + return x + drift * dt + diffusion * std::sqrt(dt) * z; + } +}; + /** * @brief Exact simulation for Geometric Brownian Motion. * diff --git a/include/lfmc/path_generator.hpp b/include/lfmc/path_generator.hpp index b09918f..1ec5e92 100644 --- a/include/lfmc/path_generator.hpp +++ b/include/lfmc/path_generator.hpp @@ -35,6 +35,7 @@ class PathGenerator { path.push_back(x); t += dt; } + paths.push_back(std::move(path)); } return paths; diff --git a/include/lfmc/payoff.hpp b/include/lfmc/payoff.hpp index d4627d5..872ea17 100644 --- a/include/lfmc/payoff.hpp +++ b/include/lfmc/payoff.hpp @@ -2,11 +2,16 @@ #include "lfmc/types.hpp" +#include +#include +#include + namespace lfmc { class Payoff { public: - virtual Payoffs generate_payoffs(const std::vector& paths) const = 0; + virtual std::expected, std::string> + generate_payoffs(const std::vector& paths) const = 0; virtual ~Payoff() = default; }; @@ -17,7 +22,8 @@ class EuropeanCall : public Payoff { public: explicit EuropeanCall(double strike); - Payoffs generate_payoffs(const std::vector& paths) const override; + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; }; class EuropeanPut : public Payoff { @@ -27,7 +33,24 @@ class EuropeanPut : public Payoff { public: explicit EuropeanPut(double strike); - Payoffs generate_payoffs(const std::vector& paths) const override; + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + +// TODO subpar in terms of architecture because we have both control variate payoff and estimator, +// but for now it's simpler to implement it this way - can enforce with a factory. Can refactor +// later if needed. +class ControlVariatePayoff : public Payoff { + private: + std::unique_ptr target_payoff_; + std::unique_ptr control_payoff_; + + public: + explicit ControlVariatePayoff(std::unique_ptr target_payoff, + std::unique_ptr control_payoff); + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; }; } // namespace lfmc diff --git a/include/lfmc/pipeline.hpp b/include/lfmc/pipeline.hpp index 949e50e..96bbb34 100644 --- a/include/lfmc/pipeline.hpp +++ b/include/lfmc/pipeline.hpp @@ -26,11 +26,25 @@ template NS> class Pipeline { std::expected run(size_t steps, double T) { while (!estimator_->converged()) { + // TODO make all the return types std::expected? IDK if it's a good idea for hot loops + // like this but it would make error handling easier and more consistent. For now just + // return auto normals = random_source_->generate_normals(steps); + if (normals.empty()) { + return std::unexpected("Failed to generate random normals"); + } + auto paths = path_generator_->generate_paths(normals, steps, T); + if (paths.empty()) { + return std::unexpected("Failed to generate paths"); + } + auto payoffs = payoff_->generate_payoffs(paths); - auto result = estimator_->add_payoffs(payoffs); + if (!payoffs) { + return std::unexpected("Failed to generate payoffs"); + } + auto result = estimator_->add_payoffs(payoffs.value()); if (!result) { return std::unexpected(result.error()); } diff --git a/include/lfmc/random_source.hpp b/include/lfmc/random_source.hpp index 2b7fe2f..66318c4 100644 --- a/include/lfmc/random_source.hpp +++ b/include/lfmc/random_source.hpp @@ -26,4 +26,18 @@ class PseudoRandomSource : public RandomSource { void seed(unsigned seed); }; +class AntitheticRandomSource : public RandomSource { + private: + std::mt19937 rng_; + std::normal_distribution dist_; + bool toggle_ = false; + + public: + AntitheticRandomSource(unsigned seed = std::random_device{}()); + + std::vector generate_normals(size_t steps, size_t) override; + + void seed(unsigned seed); +}; + } // namespace lfmc diff --git a/include/lfmc/stochastic_process.hpp b/include/lfmc/stochastic_process.hpp index 1d6271e..5251ef4 100644 --- a/include/lfmc/stochastic_process.hpp +++ b/include/lfmc/stochastic_process.hpp @@ -29,6 +29,9 @@ class GeometricBrownianMotion { double initial() const noexcept; double drift(double, double x) const noexcept; double diffusion(double, double x) const noexcept; + + double mu() const noexcept; + double sigma() const noexcept; }; } // namespace lfmc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f38176f..f167d30 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,10 +1,13 @@ # Source files for the lfmc library set(LFMC_SOURCES - euler_maruyama.cpp european_payoffs.cpp geometric_brownian_motion.cpp monte_carlo_estimator.cpp + control_variate_estimator.cpp + control_variate_payoffs.cpp pseudo_random_source.cpp + antithetic_random_source.cpp + timing.cpp ) # Create library diff --git a/src/antithetic_random_source.cpp b/src/antithetic_random_source.cpp new file mode 100644 index 0000000..b5fb9e9 --- /dev/null +++ b/src/antithetic_random_source.cpp @@ -0,0 +1,25 @@ +#include "lfmc/random_source.hpp" + +namespace lfmc { +AntitheticRandomSource::AntitheticRandomSource(unsigned seed) : rng_(seed), dist_(0.0, 1.0) {} + +std::vector AntitheticRandomSource::generate_normals(size_t steps, size_t samples) { + std::vector result(samples, Normals(steps)); + for (size_t i = 0; i < samples; ++i) { + Normals normals(steps); + for (size_t j = 0; j < steps; ++j) { + double z = dist_(rng_); + normals[j] = toggle_ ? -z : z; + } + result[i] = std::move(normals); + toggle_ = !toggle_; // Alternate between normal and antithetic + } + return result; +} + +void AntitheticRandomSource::seed(unsigned seed) { + rng_.seed(seed); + toggle_ = false; // Reset toggle when reseeding +} + +} // namespace lfmc diff --git a/src/control_variate_estimator.cpp b/src/control_variate_estimator.cpp new file mode 100644 index 0000000..678a669 --- /dev/null +++ b/src/control_variate_estimator.cpp @@ -0,0 +1,53 @@ +#include "lfmc/estimator.hpp" + +namespace lfmc { + +std::expected +ControlVariateEstimator::add_payoffs(const std::vector& payoffs) { + if (payoffs.empty()) { + return std::unexpected("No payoffs provided to add to the estimator."); + } + + for (const auto& row : payoffs) { + double x = row[0]; + double y = row[1]; + + sum_x += x; + sum_y += y; + sum_xy += x * y; + sum_yy += y * y; + ++count; + } + + return {}; +} + +bool ControlVariateEstimator::converged() const { + return count >= 10000; // Placeholder convergence criterion +} + +std::expected ControlVariateEstimator::result() const { + if (count <= 0) { + return std::unexpected("No samples added to the estimator."); + } + if (!converged()) { + return std::unexpected("Estimator has not yet converged."); + } + + double n = static_cast(count); + + double mean_x = sum_x / n; + double mean_y = sum_y / n; + + double cov_xy = (sum_xy / n) - mean_x * mean_y; + double var_y = (sum_yy / n) - mean_y * mean_y; + + if (var_y == 0.0) + return std::unexpected("Zero variance in control variate"); + + double beta = cov_xy / var_y; + + return mean_x - beta * (mean_y - control_expectation_); +} + +} // namespace lfmc diff --git a/src/control_variate_payoffs.cpp b/src/control_variate_payoffs.cpp new file mode 100644 index 0000000..9267922 --- /dev/null +++ b/src/control_variate_payoffs.cpp @@ -0,0 +1,38 @@ +#include "lfmc/payoff.hpp" + +#include +#include + +namespace lfmc { + +ControlVariatePayoff::ControlVariatePayoff(std::unique_ptr target_payoff, + std::unique_ptr control_payoff) + : target_payoff_(std::move(target_payoff)), control_payoff_(std::move(control_payoff)) {} + +std::expected, std::string> +ControlVariatePayoff::generate_payoffs(const std::vector& paths) const { + auto target_payoffs = target_payoff_->generate_payoffs(paths); + auto control_payoffs = control_payoff_->generate_payoffs(paths); + + if (!target_payoffs) { + return std::unexpected("Failed to generate target payoffs: " + target_payoffs.error()); + } else if (!control_payoffs) { + return std::unexpected("Failed to generate control payoffs: " + control_payoffs.error()); + } + + std::vector result; + for (size_t i = 0; i < target_payoffs.value().size(); ++i) { + if (target_payoffs.value()[i].size() != 1 || control_payoffs.value()[i].size() != 1) { + return std::unexpected("Each row of target payoffs must have exactly one element"); + } + + Payoffs combined_row; + combined_row.push_back(target_payoffs.value()[i][0]); // Original payoff + combined_row.push_back(control_payoffs.value()[i][0]); // Control variate payoff + result.push_back(std::move(combined_row)); + } + + return std::vector{result}; +} + +} // namespace lfmc diff --git a/src/euler_maruyama.cpp b/src/euler_maruyama.cpp deleted file mode 100644 index 2048ae3..0000000 --- a/src/euler_maruyama.cpp +++ /dev/null @@ -1,16 +0,0 @@ -#include "lfmc/stochastic_process.hpp" - -#include - -namespace lfmc { - -template class EulerMaruyama { - public: - double step(P const& process, double t, double x, double dt, double z) const noexcept { - double drift = process.drift(t, x); - double diffusion = process.diffusion(t, x); - return x + drift * dt + diffusion * std::sqrt(dt) * z; - } -}; - -} // namespace lfmc diff --git a/src/european_payoffs.cpp b/src/european_payoffs.cpp index cc2c1c3..29b5fef 100644 --- a/src/european_payoffs.cpp +++ b/src/european_payoffs.cpp @@ -1,29 +1,33 @@ #include "lfmc/payoff.hpp" #include +#include +#include namespace lfmc { EuropeanCall::EuropeanCall(double strike) : strike_(strike) {} -Payoffs EuropeanCall::generate_payoffs(const std::vector& paths) const { +std::expected, std::string> +EuropeanCall::generate_payoffs(const std::vector& paths) const { Payoffs payoffs; for (const auto& path : paths) { double final_price = path.back(); payoffs.push_back(std::max(final_price - strike_, 0.0)); } - return payoffs; + return std::vector{payoffs}; } EuropeanPut::EuropeanPut(double strike) : strike_(strike) {} -Payoffs EuropeanPut::generate_payoffs(const std::vector& paths) const { +std::expected, std::string> +EuropeanPut::generate_payoffs(const std::vector& paths) const { Payoffs payoffs; for (const auto& path : paths) { double final_price = path.back(); payoffs.push_back(std::max(strike_ - final_price, 0.0)); } - return payoffs; + return std::vector{payoffs}; } } // namespace lfmc diff --git a/src/geometric_brownian_motion.cpp b/src/geometric_brownian_motion.cpp index 98b196b..911d866 100644 --- a/src/geometric_brownian_motion.cpp +++ b/src/geometric_brownian_motion.cpp @@ -17,4 +17,12 @@ double GeometricBrownianMotion::diffusion(double, double x) const noexcept { return sigma_ * x; } +double GeometricBrownianMotion::mu() const noexcept { + return mu_; +} + +double GeometricBrownianMotion::sigma() const noexcept { + return sigma_; +} + } // namespace lfmc diff --git a/src/monte_carlo_estimator.cpp b/src/monte_carlo_estimator.cpp index 092820f..6c013d5 100644 --- a/src/monte_carlo_estimator.cpp +++ b/src/monte_carlo_estimator.cpp @@ -1,15 +1,18 @@ #include "lfmc/estimator.hpp" #include "lfmc/types.hpp" +#include + namespace lfmc { -std::expected MonteCarloEstimator::add_payoffs(const Payoffs& payoffs) { +std::expected +MonteCarloEstimator::add_payoffs(const std::vector& payoffs) { if (payoffs.empty()) { return std::unexpected("No payoffs provided to add to the estimator."); } - for (double payoff : payoffs) { - sum += payoff; + for (const auto& payoff : payoffs) { + sum += payoff[0]; ++count; } @@ -24,9 +27,9 @@ std::expected MonteCarloEstimator::result() const { if (count <= 0) { return std::unexpected("No samples added to the estimator."); } - // if (!self.converged()) { - // return std::unexpected("Estimator has not yet converged."); - // } + if (!converged()) { + return std::unexpected("Estimator has not yet converged."); + } return {sum / static_cast(count)}; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4d7f2f4..62dfd8c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -20,7 +20,7 @@ set(TEST_SOURCES test_timing.cpp test_process.cpp test_scheme.cpp - test_manager.cpp + test_pipeline.cpp ) add_executable(tests ${TEST_SOURCES}) diff --git a/tests/test_manager.cpp b/tests/test_manager.cpp deleted file mode 100644 index 08bc25e..0000000 --- a/tests/test_manager.cpp +++ /dev/null @@ -1,99 +0,0 @@ -#include "lfmc/Manager.hpp" - -#include - -TEST_CASE("Full Monte Carlo tests ", "[Manager]") { - auto blackScholesCall = [](double S, double K, double r, double sigma, double T) { - double d1 = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T) / (sigma * std::sqrt(T)); - double d2 = d1 - sigma * std::sqrt(T); - - // NormCDF approx - auto normCdf = [](double x) { return 0.5 * std::erfc(-x * M_SQRT1_2); }; - - return S * normCdf(d1) - K * std::exp(-r * T) * normCdf(d2); - }; - - SECTION("No variance reduction") { - // Parameters - double S0 = 100.0; // Initial stock price - double K = 100.0; // Strike price - double r = 0.05; // Risk-free rate (5%) - double sigma = 0.20; // Volatility (20%) - double T = 1.0; // Time to maturity (1 year) - size_t nSteps = 252; // Daily steps - - lfmc::GeometricBrownianMotion gbm(r, sigma); - lfmc::EulerMaruyama euler; - lfmc::EuropeanCall call(K); - lfmc::State state{S0, T, nSteps}; - lfmc::ManagerConfig config{1000, 0}; - - // Create manager - lfmc::Manager<> manager(gbm, euler, call, state); - - // Run sims - auto [mcPrice, stdError] = manager.simulateWithError(config); - - // Discount to present value - double discountedPrice = mcPrice * std::exp(-r * T); - double discountedError = stdError * std::exp(-r * T); - - // Calculate Black-Scholes for comparison - double bsPrice = blackScholesCall(S0, K, r, sigma, T); - - // Assertions on results - REQUIRE(std::abs(discountedPrice - bsPrice) < 0.05); // Check absolute error is small - REQUIRE((std::abs(discountedPrice - bsPrice) / bsPrice * 100) < - 10.0); // Check relative error - - // 95% confidence interval - double ciLower = discountedPrice - 1.96 * discountedError; - double ciUpper = discountedPrice + 1.96 * discountedError; - - // Check if Black-Scholes price is within confidence interval - REQUIRE(bsPrice >= ciLower); - REQUIRE(bsPrice <= ciUpper); - } - - SECTION("Antithetic variates") { - // Parameters - double S0 = 100.0; // Initial stock price - double K = 100.0; // Strike price - double r = 0.05; // Risk-free rate (5%) - double sigma = 0.20; // Volatility (20%) - double T = 1.0; // Time to maturity (1 year) - size_t nSteps = 252; // Daily steps - - lfmc::GeometricBrownianMotion gbm(r, sigma); - lfmc::EulerMaruyama euler; - lfmc::EuropeanCall call(K); - lfmc::State state{S0, T, nSteps}; - lfmc::ManagerConfig config{0, 1000}; - - // Create manager with antithetic variates - lfmc::Manager<> manager(gbm, euler, call, state); - - // Run sims with antithetic variates - auto [mcPrice, stdError] = manager.simulateWithError(config); - - // Discount to present value - double discountedPrice = mcPrice * std::exp(-r * T); - double discountedError = stdError * std::exp(-r * T); - - // Calculate Black-Scholes for comparison - double bsPrice = blackScholesCall(S0, K, r, sigma, T); - - // Assertions on results - REQUIRE(std::abs(discountedPrice - bsPrice) < 0.05); // Check absolute error is small - REQUIRE((std::abs(discountedPrice - bsPrice) / bsPrice * 100) < - 10.0); // Check relative error - - // 95% confidence interval - double ciLower = discountedPrice - 1.96 * discountedError; - double ciUpper = discountedPrice + 1.96 * discountedError; - - // Check if Black-Scholes price is within confidence interval - REQUIRE(bsPrice >= ciLower); - REQUIRE(bsPrice <= ciUpper); - } -} diff --git a/tests/test_pipeline.cpp b/tests/test_pipeline.cpp new file mode 100644 index 0000000..6c707fa --- /dev/null +++ b/tests/test_pipeline.cpp @@ -0,0 +1,88 @@ +#include "lfmc/pipeline.hpp" + +#include + +using namespace lfmc; + +// TODO add more test cases, e.g. for convergence criteria, error handling, etc. and dummy types + +/* =========================== + Test Cases + =========================== */ + +TEST_CASE("Pipeline runs until estimator converges") { + GeometricBrownianMotion gbm(0.05, 0.2, 100.0); + EulerMaruyama euler; + + auto rs = std::make_unique(); + auto pg = std::make_unique< + PathGenerator>>(gbm, euler); + auto po = std::make_unique(100.0); + auto est = std::make_unique(); + + Pipeline> pipeline( + std::move(rs), std::move(pg), std::move(po), std::move(est)); + + auto result = pipeline.run(10, 1.0); + + INFO("Pipeline result: " << (result.has_value() ? std::to_string(result.value()) + : result.error())); + REQUIRE(result.has_value()); + REQUIRE(result.value() > 0.0); // Price should be positive +} + +TEST_CASE("Pipeline with antithetic variates") { + GeometricBrownianMotion gbm(0.05, 0.2, 100.0); + EulerMaruyama euler; + + auto rs = std::make_unique(); + auto pg = std::make_unique< + PathGenerator>>(gbm, euler); + auto po = std::make_unique(100.0); + auto est = std::make_unique(); + + Pipeline> pipeline( + std::move(rs), std::move(pg), std::move(po), std::move(est)); + + auto result = pipeline.run(10, 1.0); + + INFO("Pipeline result: " << (result.has_value() ? std::to_string(result.value()) + : result.error())); + REQUIRE(result.has_value()); + REQUIRE(result.value() > 0.0); // Price should be positive +} + +TEST_CASE("Pipeline with control variates") { + GeometricBrownianMotion gbm(0.05, 0.2, 100.0); + EulerMaruyama euler; + + auto rs = std::make_unique(); + auto pg = std::make_unique< + PathGenerator>>(gbm, euler); + auto po1 = std::make_unique(100.0); + auto po2 = std::make_unique( + 100.0); // Control variate payoff (same as target for simplicity) + auto po = std::make_unique(std::move(po1), std::move(po2)); + auto est = std::make_unique( + gbm.initial() * std::exp(gbm.mu() * 1.0)); // Control variate expectation + + Pipeline> pipeline( + std::move(rs), std::move(pg), std::move(po), std::move(est)); + + auto result = pipeline.run(10, 1.0); + + INFO("Pipeline result: " << (result.has_value() ? std::to_string(result.value()) + : result.error())); + REQUIRE(result.has_value()); + REQUIRE(result.value() > 0.0); // Price should be positive + // TODO Currently around 100, but should be around 10 - need to + // debug control variate implementation +} + +TEST_CASE("Pipeline stops early if estimator add_payoffs fails") { + // TODO +} + +TEST_CASE("Pipeline performs correct number of iterations") { + // TODO +} diff --git a/tests/test_process.cpp b/tests/test_process.cpp index 869b886..74d5109 100644 --- a/tests/test_process.cpp +++ b/tests/test_process.cpp @@ -1,4 +1,4 @@ -#include "lfmc/StochasticProcess.hpp" +#include "lfmc/stochastic_process.hpp" #include #include @@ -6,43 +6,47 @@ using Catch::Matchers::WithinAbs; TEST_CASE("StochasticProcess concept works correctly", "[StochasticProcess]") { - struct ValidProcess { - double drift(double x) const noexcept { - return x; - } - double diffusion(double x) const noexcept { - return x; - } - }; - - struct InvalidProcessNoDrift { - double diffusion(double x) const noexcept { - return x; - } - }; - - struct InvalidProcessWrongDiffusion { - double drift(double x) const noexcept { - return x; - } - int diffusion(double x) const noexcept { - return static_cast(x); - } - }; - - REQUIRE(lfmc::StochasticProcess); - REQUIRE_FALSE(lfmc::StochasticProcess); - REQUIRE_FALSE(lfmc::StochasticProcess); + // TODO + // struct ValidProcess { + // double initial() const noexcept { + // return 1.0; + // } + // double drift(double x) const noexcept { + // return x; + // } + // double diffusion(double x) const noexcept { + // return x; + // } + // }; + // + // struct InvalidProcessNoDrift { + // double diffusion(double x) const noexcept { + // return x; + // } + // }; + // + // struct InvalidProcessWrongDiffusion { + // double drift(double x) const noexcept { + // return x; + // } + // int diffusion(double x) const noexcept { + // return static_cast(x); + // } + // }; + // + // REQUIRE(lfmc::StochasticProcess); + // REQUIRE_FALSE(lfmc::StochasticProcess); + // REQUIRE_FALSE(lfmc::StochasticProcess); } TEST_CASE("GeometricBrownianMotion computes drift and diffusion correctly", "[StochasticProcess][GeometricBrownianMotion]") { - lfmc::GeometricBrownianMotion gbm{0.1, 0.2}; + lfmc::GeometricBrownianMotion gbm(0.1, 0.2, 100.0); double x = 100.0; double expectedDrift = 0.1 * x; double expectedDiffusion = 0.2 * x; - REQUIRE_THAT(gbm.drift(x), WithinAbs(expectedDrift, 1e-10)); - REQUIRE_THAT(gbm.diffusion(x), WithinAbs(expectedDiffusion, 1e-10)); + REQUIRE_THAT(gbm.drift(x, x), WithinAbs(expectedDrift, 1e-10)); + REQUIRE_THAT(gbm.diffusion(x, x), WithinAbs(expectedDiffusion, 1e-10)); } diff --git a/tests/test_scheme.cpp b/tests/test_scheme.cpp index 5629ed3..b1d3678 100644 --- a/tests/test_scheme.cpp +++ b/tests/test_scheme.cpp @@ -1,5 +1,5 @@ -#include "lfmc/NumericalScheme.hpp" -#include "lfmc/StochasticProcess.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/stochastic_process.hpp" #include #include @@ -9,6 +9,9 @@ using Catch::Matchers::WithinAbs; namespace test1 { struct ValidProcess { + double initial() const noexcept { + return 1.0; + } double drift(double x) const noexcept { return x; } @@ -19,7 +22,7 @@ struct ValidProcess { template struct ValidScheme { double step(P const& process, double x, double dt, double dW) const noexcept { - return x + process.drift(x) * dt + process.diffusion(x) * dW; + return x + process.drift(x) * dt + process.diffusion(x) * dW * std::sqrt(dt); } }; @@ -28,8 +31,8 @@ struct InvalidSchemeNoStep { }; template struct InvalidSchemeWrongStep { - int step(P const& process, double x, double dt, double dW) const noexcept { - return static_cast(x); + double step(P const& process, double x, double dt, double dW) const noexcept { + return x + process.drift(x) * dt; // Missing diffusion term } }; @@ -38,45 +41,50 @@ template struct InvalidSchemeWrongStep { TEST_CASE("NumericalScheme concept works correctly", "[NumericalScheme]") { using namespace test1; - REQUIRE(lfmc::NumericalScheme, ValidProcess>); - REQUIRE_FALSE(lfmc::NumericalScheme); - REQUIRE_FALSE(lfmc::NumericalScheme, ValidProcess>); + // TODO + // REQUIRE(lfmc::NumericalScheme>); + // REQUIRE_FALSE(lfmc::NumericalScheme); + // REQUIRE_FALSE(lfmc::NumericalScheme, ValidProcess>); }; TEST_CASE("EulerMaruyama computes next step correctly", "[NumericalScheme][EulerMaruyama]") { - lfmc::GeometricBrownianMotion gbm{0.1, 0.2}; - lfmc::EulerMaruyama scheme; - double x = 100.0; double dt = 0.01; double dW = 0.05; - double expectedNextX = x + gbm.drift(x) * dt + gbm.diffusion(x) * dW * std::sqrt(dt); - double computedNextX = scheme.step(gbm, x, dt, dW); + lfmc::GeometricBrownianMotion gbm(0.1, 0.2, 100.0); + lfmc::EulerMaruyama scheme; + + double expectedNextX = x + gbm.drift(x, x) * dt + gbm.diffusion(x, x) * dW * std::sqrt(dt); + double computedNextX = scheme.step(gbm, x, x, dt, dW); REQUIRE_THAT(computedNextX, WithinAbs(expectedNextX, 1e-10)); } TEST_CASE("EulerMaruyama works with different stochastic processes", "[NumericalScheme][EulerMaruyama]") { - struct CustomProcess { - double drift(double x) const noexcept { - return 2.0 * x; - } - double diffusion(double x) const noexcept { - return 0.5 * x; - } - }; - - CustomProcess process; - lfmc::EulerMaruyama scheme; - - double x = 50.0; - double dt = 0.02; - double dW = 0.03; - - double expectedNextX = x + process.drift(x) * dt + process.diffusion(x) * dW * std::sqrt(dt); - double computedNextX = scheme.step(process, x, dt, dW); - - REQUIRE_THAT(computedNextX, WithinAbs(expectedNextX, 1e-10)); + // TODO + // struct CustomProcess { + // double initial() const noexcept { + // return 50.0; + // } + // double drift(double x) const noexcept { + // return 2.0 * x; + // } + // double diffusion(double x) const noexcept { + // return 0.5 * x; + // } + // }; + // + // CustomProcess process; + // lfmc::EulerMaruyama scheme; + // + // double x = 50.0; + // double dt = 0.02; + // double dW = 0.03; + // + // double expectedNextX = x + process.drift(x) * dt + process.diffusion(x) * dW * std::sqrt(dt); + // double computedNextX = scheme.step(process, x, dt, dW); + // + // REQUIRE_THAT(computedNextX, WithinAbs(expectedNextX, 1e-10)); } From bd7f43376eda5e7f0f92b9f146fbb1175acb5163 Mon Sep 17 00:00:00 2001 From: Alexander Robbins Date: Tue, 3 Mar 2026 11:09:19 -0500 Subject: [PATCH 07/12] added Different payoff structures --- src/payoffs/asian_payoffs.cpp | 76 +++++++++++ src/payoffs/barrier_payoffs.cpp | 89 +++++++++++++ src/payoffs/lookback_payoffs.cpp | 64 +++++++++ tests/CMakeLists.txt | 1 + tests/test_convergency.cpp | 216 +++++++++++++++++++++++++++++++ 5 files changed, 446 insertions(+) create mode 100644 src/payoffs/asian_payoffs.cpp create mode 100644 src/payoffs/barrier_payoffs.cpp create mode 100644 src/payoffs/lookback_payoffs.cpp create mode 100644 tests/test_convergency.cpp diff --git a/src/payoffs/asian_payoffs.cpp b/src/payoffs/asian_payoffs.cpp new file mode 100644 index 0000000..d56de82 --- /dev/null +++ b/src/payoffs/asian_payoffs.cpp @@ -0,0 +1,76 @@ +#pragma once + +#include "lfmc/payoff.hpp" +#include "lfmc/types.hpp" + +#include +#include +#include +#include +#include + +namespace lfmc { + + + + +class AsianCall : public Payoff { +private: + double strike_; + +public: + explicit AsianCall(double strike) : strike_(strike) {} + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override { + Payoffs payoffs; + + payoffs.reserve(paths.size()); + + for (const auto& path : paths) { + + if (path.empty()) + return std::unexpected("Empty path encountered in AsianCall"); + + double mean = std::reduce(path.begin(), path.end(), 0.0) + / static_cast(path.size()); + payoffs.push_back(std::max(mean - strike_, 0.0)); + } + + return std::vector{payoffs}; + } +}; + + +class AsianPut : public Payoff { +private: + double strike_; + +public: + explicit AsianPut(double strike) : strike_(strike) {} + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override { + + + Payoffs payoffs; + + + + payoffs.reserve(paths.size()); + + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in AsianPut"); + + + double mean = std::reduce(path.begin(), path.end(), 0.0) + / static_cast(path.size()); + payoffs.push_back(std::max(strike_ - mean, 0.0)); + } + + return std::vector{payoffs}; + } +}; + +} \ No newline at end of file diff --git a/src/payoffs/barrier_payoffs.cpp b/src/payoffs/barrier_payoffs.cpp new file mode 100644 index 0000000..2454470 --- /dev/null +++ b/src/payoffs/barrier_payoffs.cpp @@ -0,0 +1,89 @@ +#pragma once + +#include "lfmc/payoff.hpp" +#include "lfmc/types.hpp" + +#include +#include +#include +#include +#include + +namespace lfmc { + +class UpAndOutCall : public Payoff { +private: + double strike_; + double barrier_; + +public: + UpAndOutCall(double strike, double barrier) + : strike_(strike), barrier_(barrier) {} + + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override { + + + Payoffs payoffs; + payoffs.reserve(paths.size()); + + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in UpAndOutCall"); + + bool knocked_out = std::any_of(path.begin(), path.end(), + [this](double s) { return s >= barrier_; }); + + + + + if (knocked_out) { + payoffs.push_back(0.0); + } else { + payoffs.push_back(std::max(path.back() - strike_, 0.0)); + } + } + + return std::vector{payoffs}; + } +}; + + + +class DownAndInPut : public Payoff { +private: + double strike_; + double barrier_; + +public: + DownAndInPut(double strike, double barrier) + : strike_(strike), barrier_(barrier) {} + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override { + Payoffs payoffs; + payoffs.reserve(paths.size()); + + for (const auto& path : paths) { + + + if (path.empty()) + return std::unexpected("Empty path encountered in DownAndInPut"); + + bool knocked_in = std::any_of(path.begin(), path.end(), + [this](double s) { return s <= barrier_; }); + + if (knocked_in) { + payoffs.push_back(std::max(strike_ - path.back(), 0.0)); + + + } else { + payoffs.push_back(0.0); + } + } + + return std::vector{payoffs}; + } +}; +} \ No newline at end of file diff --git a/src/payoffs/lookback_payoffs.cpp b/src/payoffs/lookback_payoffs.cpp new file mode 100644 index 0000000..df8376f --- /dev/null +++ b/src/payoffs/lookback_payoffs.cpp @@ -0,0 +1,64 @@ +#pragma once + +#include "lfmc/payoff.hpp" +#include "lfmc/types.hpp" + +#include +#include +#include +#include +#include + +namespace lfmc { + + + class LookbackCall : public Payoff { +public: + std::expected, std::string> + + + + generate_payoffs(const std::vector& paths) const override { + Payoffs payoffs; + payoffs.reserve(paths.size()); + + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in LookbackCall"); + + + + double min_price = *std::min_element(path.begin(), path.end()); + payoffs.push_back(path.back() - min_price); + } + + return std::vector{payoffs}; + } +}; + + + +class LookbackPut : public Payoff { +public: + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override { + Payoffs payoffs; + payoffs.reserve(paths.size()); + + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in LookbackPut"); + + double max_price = *std::max_element(path.begin(), path.end()); + + + + payoffs.push_back(max_price - path.back()); + } + + return std::vector{payoffs}; + } +}; + +} \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 62dfd8c..940f9fe 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,6 +21,7 @@ set(TEST_SOURCES test_process.cpp test_scheme.cpp test_pipeline.cpp + test_convergency.cpp ) add_executable(tests ${TEST_SOURCES}) diff --git a/tests/test_convergency.cpp b/tests/test_convergency.cpp new file mode 100644 index 0000000..98d2728 --- /dev/null +++ b/tests/test_convergency.cpp @@ -0,0 +1,216 @@ +#include "lfmc/pipeline.hpp" +#include "lfmc/payoffs/asian_payoffs.hpp" +#include "lfmc/payoffs/barrier_payoffs.hpp" +#include "lfmc/payoffs/lookback_payoffs.hpp" +#include "lfmc/timing.hpp" + +#include +#include +#include +#include +#include +#include + +using namespace lfmc; +using Catch::Matchers::WithinAbs; + +// ---------------------------------------------------------------- +// Black-Scholes closed-form helpers for ground truth +// ---------------------------------------------------------------- +namespace bs { + +static double norm_cdf(double x) { + return 0.5 * std::erfc(-x / std::sqrt(2.0)); +} + +// Standard European Call used to sanity check our GBM setup +double european_call(double S, double K, double r, double sigma, double T) { + double d1 = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T) + / (sigma * std::sqrt(T)); + double d2 = d1 - sigma * std::sqrt(T); + return S * norm_cdf(d1) - K * std::exp(-r * T) * norm_cdf(d2); +} + +// Geometric Asian Call closed-form (Kemna-Vorst approximation) +// Used as ground truth since arithmetic Asian has no closed form +double geometric_asian_call(double S, double K, double r, + double sigma, double T, int n) { + double sigma_adj = sigma * std::sqrt((2.0 * n + 1.0) / (6.0 * (n + 1.0))); + double r_adj = 0.5 * (r - 0.5 * sigma * sigma) + + 0.5 * sigma_adj * sigma_adj; + double d1 = (std::log(S / K) + (r_adj + 0.5 * sigma_adj * sigma_adj) * T) + / (sigma_adj * std::sqrt(T)); + double d2 = d1 - sigma_adj * std::sqrt(T); + return std::exp(-r * T) + * (S * std::exp(r_adj * T) * norm_cdf(d1) - K * norm_cdf(d2)); +} + +// Up-and-Out Call closed form (continuous barrier, no dividends) +double up_and_out_call(double S, double K, double H, + double r, double sigma, double T) { + if (S >= H) return 0.0; + double vanilla = european_call(S, K, r, sigma, T); + + double lambda = (r + 0.5 * sigma * sigma) / (sigma * sigma); + double d1 = (std::log(H * H / (S * K)) + (r + 0.5 * sigma * sigma) * T) + / (sigma * std::sqrt(T)); + double d2 = d1 - sigma * std::sqrt(T); + double reflection = std::pow(H / S, 2.0 * lambda) + * (S * norm_cdf(d1) - K * std::exp(-r * T) * norm_cdf(d2)); + + return vanilla - reflection; +} + +} // namespace bs + +// ---------------------------------------------------------------- +// Convergence runner: runs pipeline at increasing sample counts +// and prints a table of results vs ground truth +// ---------------------------------------------------------------- +struct ConvergenceResult { + std::size_t samples; + double estimate; + double ground_truth; + double abs_error; + long long elapsed_ms; +}; + +// We run the pipeline once per sample tier by re-seeding with +// a fixed seed for reproducibility +template +std::vector run_convergence( + const std::string& label, + PayoffFactory make_payoff, + double ground_truth, + const std::vector& sample_tiers, + double S0 = 100.0, + double mu = 0.05, + double sigma = 0.20, + std::size_t steps = 252, + double T = 1.0) +{ + std::cout << "\n=== " << label << " ===\n"; + std::cout << std::setw(10) << "Samples" + << std::setw(14) << "Estimate" + << std::setw(14) << "Truth" + << std::setw(14) << "AbsError" + << std::setw(12) << "Time(ms)" + << "\n" + << std::string(64, '-') << "\n"; + + std::vector results; + + for (std::size_t n : sample_tiers) { + GeometricBrownianMotion gbm(mu, sigma, S0); + EulerMaruyama euler; + + // Fixed seed for reproducibility + auto rs = std::make_unique(42u); + auto pg = std::make_unique< + PathGenerator>>(gbm, euler); + auto po = make_payoff(); + auto est = std::make_unique(n); // see note below + + Pipeline> pipeline( + std::move(rs), std::move(pg), + std::move(po), std::move(est)); + + long long elapsed = 0; + double estimate = 0.0; + { + ScopedTimer timer(elapsed); + auto res = pipeline.run(steps, T); + REQUIRE(res.has_value()); + estimate = res.value(); + } + + double err = std::abs(estimate - ground_truth); + results.push_back({n, estimate, ground_truth, err, elapsed}); + + std::cout << std::setw(10) << n + << std::setw(14) << std::fixed << std::setprecision(4) << estimate + << std::setw(14) << ground_truth + << std::setw(14) << err + << std::setw(12) << elapsed + << "\n"; + } + return results; +} + +// ---------------------------------------------------------------- +// Shared parameters +// ---------------------------------------------------------------- +static constexpr double S0 = 100.0; +static constexpr double K = 100.0; +static constexpr double B_UP = 120.0; // Up-and-out barrier +static constexpr double B_DN = 80.0; // Down-and-in barrier +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.20; +static constexpr double T = 1.0; +static constexpr int STEPS = 252; + +static const std::vector TIERS = + {1000, 5000, 10000, 50000, 100000, 500000}; + +// ---------------------------------------------------------------- +// Tests +// ---------------------------------------------------------------- + +TEST_CASE("Asian Call convergence", "[exotic][convergence][asian]") { + double truth = bs::geometric_asian_call(S0, K, MU, SIGMA, T, STEPS); + + auto results = run_convergence( + "Arithmetic Asian Call (truth = geometric approx)", + []() { return std::make_unique(K); }, + truth, TIERS); + + // At 100k samples we should be within $0.50 of the approximation + auto& last = results.back(); + REQUIRE_THAT(last.estimate, WithinAbs(last.ground_truth, 0.50)); +} + +TEST_CASE("Up-and-Out Barrier Call convergence", + "[exotic][convergence][barrier]") { + double truth = bs::up_and_out_call(S0, K, B_UP, MU, SIGMA, T); + + auto results = run_convergence( + "Up-and-Out Barrier Call", + []() { return std::make_unique(K, B_UP); }, + truth, TIERS); + + auto& last = results.back(); + REQUIRE_THAT(last.estimate, WithinAbs(last.ground_truth, 0.50)); +} + +TEST_CASE("Lookback Call convergence", "[exotic][convergence][lookback]") { + // No simple closed form for discrete lookback we use the + // large-sample MC estimate itself as a self-consistency check + // and just verify convergence tightens with more samples + auto results = run_convergence( + "Lookback Call (floating strike)", + []() { return std::make_unique(); }, + 0.0, // placeholder see check below + TIERS); + + // Check that later estimates are closer to each other than early ones + // i.e. the std deviation of the last 3 tiers < first 3 tiers + auto spread = [&](int a, int b) { + return std::abs(results[b].estimate - results[a].estimate); + }; + REQUIRE(spread(3, 5) < spread(0, 2)); +} + +TEST_CASE("Sanity check: European Call matches Black-Scholes", + "[sanity][european]") { + double bs_price = bs::european_call(S0, K, MU, SIGMA, T); + + auto results = run_convergence( + "European Call (BS sanity check)", + []() { return std::make_unique(K); }, + bs_price, TIERS); + + auto& last = results.back(); + REQUIRE_THAT(last.estimate, WithinAbs(last.ground_truth, 0.20)); +} \ No newline at end of file From 8cf6970d05e55cfdbf8c93aa57296932b2900984 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 3 Mar 2026 16:10:21 +0000 Subject: [PATCH 08/12] style: apply clang-format --- src/payoffs/asian_payoffs.cpp | 38 ++++---- src/payoffs/barrier_payoffs.cpp | 43 ++++------ src/payoffs/lookback_payoffs.cpp | 24 ++---- tests/test_convergency.cpp | 143 ++++++++++++------------------- 4 files changed, 95 insertions(+), 153 deletions(-) diff --git a/src/payoffs/asian_payoffs.cpp b/src/payoffs/asian_payoffs.cpp index d56de82..e492275 100644 --- a/src/payoffs/asian_payoffs.cpp +++ b/src/payoffs/asian_payoffs.cpp @@ -3,37 +3,34 @@ #include "lfmc/payoff.hpp" #include "lfmc/types.hpp" +#include #include #include -#include -#include #include +#include namespace lfmc { - - - class AsianCall : public Payoff { -private: + private: double strike_; -public: + public: explicit AsianCall(double strike) : strike_(strike) {} std::expected, std::string> generate_payoffs(const std::vector& paths) const override { Payoffs payoffs; - + payoffs.reserve(paths.size()); for (const auto& path : paths) { - + if (path.empty()) return std::unexpected("Empty path encountered in AsianCall"); - double mean = std::reduce(path.begin(), path.end(), 0.0) - / static_cast(path.size()); + double mean = + std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); payoffs.push_back(std::max(mean - strike_, 0.0)); } @@ -41,31 +38,26 @@ class AsianCall : public Payoff { } }; - class AsianPut : public Payoff { -private: + private: double strike_; -public: + public: explicit AsianPut(double strike) : strike_(strike) {} std::expected, std::string> generate_payoffs(const std::vector& paths) const override { - - + Payoffs payoffs; - - - + payoffs.reserve(paths.size()); for (const auto& path : paths) { if (path.empty()) return std::unexpected("Empty path encountered in AsianPut"); - - double mean = std::reduce(path.begin(), path.end(), 0.0) - / static_cast(path.size()); + double mean = + std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); payoffs.push_back(std::max(strike_ - mean, 0.0)); } @@ -73,4 +65,4 @@ class AsianPut : public Payoff { } }; -} \ No newline at end of file +} // namespace lfmc \ No newline at end of file diff --git a/src/payoffs/barrier_payoffs.cpp b/src/payoffs/barrier_payoffs.cpp index 2454470..df00fa1 100644 --- a/src/payoffs/barrier_payoffs.cpp +++ b/src/payoffs/barrier_payoffs.cpp @@ -3,28 +3,25 @@ #include "lfmc/payoff.hpp" #include "lfmc/types.hpp" +#include #include #include -#include -#include #include +#include namespace lfmc { class UpAndOutCall : public Payoff { -private: + private: double strike_; double barrier_; -public: - UpAndOutCall(double strike, double barrier) - : strike_(strike), barrier_(barrier) {} - + public: + UpAndOutCall(double strike, double barrier) : strike_(strike), barrier_(barrier) {} std::expected, std::string> generate_payoffs(const std::vector& paths) const override { - - + Payoffs payoffs; payoffs.reserve(paths.size()); @@ -32,11 +29,8 @@ class UpAndOutCall : public Payoff { if (path.empty()) return std::unexpected("Empty path encountered in UpAndOutCall"); - bool knocked_out = std::any_of(path.begin(), path.end(), - [this](double s) { return s >= barrier_; }); - - - + bool knocked_out = + std::any_of(path.begin(), path.end(), [this](double s) { return s >= barrier_; }); if (knocked_out) { payoffs.push_back(0.0); @@ -49,16 +43,13 @@ class UpAndOutCall : public Payoff { } }; - - class DownAndInPut : public Payoff { -private: + private: double strike_; double barrier_; -public: - DownAndInPut(double strike, double barrier) - : strike_(strike), barrier_(barrier) {} + public: + DownAndInPut(double strike, double barrier) : strike_(strike), barrier_(barrier) {} std::expected, std::string> generate_payoffs(const std::vector& paths) const override { @@ -66,18 +57,16 @@ class DownAndInPut : public Payoff { payoffs.reserve(paths.size()); for (const auto& path : paths) { - - + if (path.empty()) return std::unexpected("Empty path encountered in DownAndInPut"); - bool knocked_in = std::any_of(path.begin(), path.end(), - [this](double s) { return s <= barrier_; }); + bool knocked_in = + std::any_of(path.begin(), path.end(), [this](double s) { return s <= barrier_; }); if (knocked_in) { payoffs.push_back(std::max(strike_ - path.back(), 0.0)); - - + } else { payoffs.push_back(0.0); } @@ -86,4 +75,4 @@ class DownAndInPut : public Payoff { return std::vector{payoffs}; } }; -} \ No newline at end of file +} // namespace lfmc \ No newline at end of file diff --git a/src/payoffs/lookback_payoffs.cpp b/src/payoffs/lookback_payoffs.cpp index df8376f..345f1c0 100644 --- a/src/payoffs/lookback_payoffs.cpp +++ b/src/payoffs/lookback_payoffs.cpp @@ -3,21 +3,18 @@ #include "lfmc/payoff.hpp" #include "lfmc/types.hpp" +#include #include #include -#include -#include #include +#include namespace lfmc { - - class LookbackCall : public Payoff { -public: +class LookbackCall : public Payoff { + public: std::expected, std::string> - - generate_payoffs(const std::vector& paths) const override { Payoffs payoffs; payoffs.reserve(paths.size()); @@ -26,8 +23,6 @@ namespace lfmc { if (path.empty()) return std::unexpected("Empty path encountered in LookbackCall"); - - double min_price = *std::min_element(path.begin(), path.end()); payoffs.push_back(path.back() - min_price); } @@ -36,11 +31,8 @@ namespace lfmc { } }; - - class LookbackPut : public Payoff { -public: - + public: std::expected, std::string> generate_payoffs(const std::vector& paths) const override { Payoffs payoffs; @@ -51,9 +43,7 @@ class LookbackPut : public Payoff { return std::unexpected("Empty path encountered in LookbackPut"); double max_price = *std::max_element(path.begin(), path.end()); - - - + payoffs.push_back(max_price - path.back()); } @@ -61,4 +51,4 @@ class LookbackPut : public Payoff { } }; -} \ No newline at end of file +} // namespace lfmc \ No newline at end of file diff --git a/tests/test_convergency.cpp b/tests/test_convergency.cpp index 98d2728..6db1724 100644 --- a/tests/test_convergency.cpp +++ b/tests/test_convergency.cpp @@ -1,15 +1,15 @@ -#include "lfmc/pipeline.hpp" #include "lfmc/payoffs/asian_payoffs.hpp" #include "lfmc/payoffs/barrier_payoffs.hpp" #include "lfmc/payoffs/lookback_payoffs.hpp" +#include "lfmc/pipeline.hpp" #include "lfmc/timing.hpp" #include #include -#include +#include #include +#include #include -#include using namespace lfmc; using Catch::Matchers::WithinAbs; @@ -25,38 +25,34 @@ static double norm_cdf(double x) { // Standard European Call used to sanity check our GBM setup double european_call(double S, double K, double r, double sigma, double T) { - double d1 = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T) - / (sigma * std::sqrt(T)); + double d1 = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T) / (sigma * std::sqrt(T)); double d2 = d1 - sigma * std::sqrt(T); return S * norm_cdf(d1) - K * std::exp(-r * T) * norm_cdf(d2); } // Geometric Asian Call closed-form (Kemna-Vorst approximation) // Used as ground truth since arithmetic Asian has no closed form -double geometric_asian_call(double S, double K, double r, - double sigma, double T, int n) { +double geometric_asian_call(double S, double K, double r, double sigma, double T, int n) { double sigma_adj = sigma * std::sqrt((2.0 * n + 1.0) / (6.0 * (n + 1.0))); - double r_adj = 0.5 * (r - 0.5 * sigma * sigma) - + 0.5 * sigma_adj * sigma_adj; - double d1 = (std::log(S / K) + (r_adj + 0.5 * sigma_adj * sigma_adj) * T) - / (sigma_adj * std::sqrt(T)); + double r_adj = 0.5 * (r - 0.5 * sigma * sigma) + 0.5 * sigma_adj * sigma_adj; + double d1 = + (std::log(S / K) + (r_adj + 0.5 * sigma_adj * sigma_adj) * T) / (sigma_adj * std::sqrt(T)); double d2 = d1 - sigma_adj * std::sqrt(T); - return std::exp(-r * T) - * (S * std::exp(r_adj * T) * norm_cdf(d1) - K * norm_cdf(d2)); + return std::exp(-r * T) * (S * std::exp(r_adj * T) * norm_cdf(d1) - K * norm_cdf(d2)); } // Up-and-Out Call closed form (continuous barrier, no dividends) -double up_and_out_call(double S, double K, double H, - double r, double sigma, double T) { - if (S >= H) return 0.0; +double up_and_out_call(double S, double K, double H, double r, double sigma, double T) { + if (S >= H) + return 0.0; double vanilla = european_call(S, K, r, sigma, T); double lambda = (r + 0.5 * sigma * sigma) / (sigma * sigma); - double d1 = (std::log(H * H / (S * K)) + (r + 0.5 * sigma * sigma) * T) - / (sigma * std::sqrt(T)); + double d1 = + (std::log(H * H / (S * K)) + (r + 0.5 * sigma * sigma) * T) / (sigma * std::sqrt(T)); double d2 = d1 - sigma * std::sqrt(T); - double reflection = std::pow(H / S, 2.0 * lambda) - * (S * norm_cdf(d1) - K * std::exp(-r * T) * norm_cdf(d2)); + double reflection = + std::pow(H / S, 2.0 * lambda) * (S * norm_cdf(d1) - K * std::exp(-r * T) * norm_cdf(d2)); return vanilla - reflection; } @@ -68,34 +64,23 @@ double up_and_out_call(double S, double K, double H, // and prints a table of results vs ground truth // ---------------------------------------------------------------- struct ConvergenceResult { - std::size_t samples; - double estimate; - double ground_truth; - double abs_error; - long long elapsed_ms; + std::size_t samples; + double estimate; + double ground_truth; + double abs_error; + long long elapsed_ms; }; // We run the pipeline once per sample tier by re-seeding with // a fixed seed for reproducibility template -std::vector run_convergence( - const std::string& label, - PayoffFactory make_payoff, - double ground_truth, - const std::vector& sample_tiers, - double S0 = 100.0, - double mu = 0.05, - double sigma = 0.20, - std::size_t steps = 252, - double T = 1.0) -{ +std::vector +run_convergence(const std::string& label, PayoffFactory make_payoff, double ground_truth, + const std::vector& sample_tiers, double S0 = 100.0, double mu = 0.05, + double sigma = 0.20, std::size_t steps = 252, double T = 1.0) { std::cout << "\n=== " << label << " ===\n"; - std::cout << std::setw(10) << "Samples" - << std::setw(14) << "Estimate" - << std::setw(14) << "Truth" - << std::setw(14) << "AbsError" - << std::setw(12) << "Time(ms)" - << "\n" + std::cout << std::setw(10) << "Samples" << std::setw(14) << "Estimate" << std::setw(14) + << "Truth" << std::setw(14) << "AbsError" << std::setw(12) << "Time(ms)" << "\n" << std::string(64, '-') << "\n"; std::vector results; @@ -105,20 +90,18 @@ std::vector run_convergence( EulerMaruyama euler; // Fixed seed for reproducibility - auto rs = std::make_unique(42u); - auto pg = std::make_unique< - PathGenerator>>(gbm, euler); - auto po = make_payoff(); + auto rs = std::make_unique(42u); + auto pg = std::make_unique< + PathGenerator>>(gbm, + euler); + auto po = make_payoff(); auto est = std::make_unique(n); // see note below - Pipeline> pipeline( - std::move(rs), std::move(pg), - std::move(po), std::move(est)); + Pipeline> pipeline( + std::move(rs), std::move(pg), std::move(po), std::move(est)); long long elapsed = 0; - double estimate = 0.0; + double estimate = 0.0; { ScopedTimer timer(elapsed); auto res = pipeline.run(steps, T); @@ -129,12 +112,9 @@ std::vector run_convergence( double err = std::abs(estimate - ground_truth); results.push_back({n, estimate, ground_truth, err, elapsed}); - std::cout << std::setw(10) << n - << std::setw(14) << std::fixed << std::setprecision(4) << estimate - << std::setw(14) << ground_truth - << std::setw(14) << err - << std::setw(12) << elapsed - << "\n"; + std::cout << std::setw(10) << n << std::setw(14) << std::fixed << std::setprecision(4) + << estimate << std::setw(14) << ground_truth << std::setw(14) << err + << std::setw(12) << elapsed << "\n"; } return results; } @@ -142,17 +122,16 @@ std::vector run_convergence( // ---------------------------------------------------------------- // Shared parameters // ---------------------------------------------------------------- -static constexpr double S0 = 100.0; -static constexpr double K = 100.0; -static constexpr double B_UP = 120.0; // Up-and-out barrier -static constexpr double B_DN = 80.0; // Down-and-in barrier -static constexpr double MU = 0.05; -static constexpr double SIGMA = 0.20; -static constexpr double T = 1.0; -static constexpr int STEPS = 252; - -static const std::vector TIERS = - {1000, 5000, 10000, 50000, 100000, 500000}; +static constexpr double S0 = 100.0; +static constexpr double K = 100.0; +static constexpr double B_UP = 120.0; // Up-and-out barrier +static constexpr double B_DN = 80.0; // Down-and-in barrier +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.20; +static constexpr double T = 1.0; +static constexpr int STEPS = 252; + +static const std::vector TIERS = {1000, 5000, 10000, 50000, 100000, 500000}; // ---------------------------------------------------------------- // Tests @@ -163,22 +142,19 @@ TEST_CASE("Asian Call convergence", "[exotic][convergence][asian]") { auto results = run_convergence( "Arithmetic Asian Call (truth = geometric approx)", - []() { return std::make_unique(K); }, - truth, TIERS); + []() { return std::make_unique(K); }, truth, TIERS); // At 100k samples we should be within $0.50 of the approximation auto& last = results.back(); REQUIRE_THAT(last.estimate, WithinAbs(last.ground_truth, 0.50)); } -TEST_CASE("Up-and-Out Barrier Call convergence", - "[exotic][convergence][barrier]") { +TEST_CASE("Up-and-Out Barrier Call convergence", "[exotic][convergence][barrier]") { double truth = bs::up_and_out_call(S0, K, B_UP, MU, SIGMA, T); auto results = run_convergence( - "Up-and-Out Barrier Call", - []() { return std::make_unique(K, B_UP); }, - truth, TIERS); + "Up-and-Out Barrier Call", []() { return std::make_unique(K, B_UP); }, truth, + TIERS); auto& last = results.back(); REQUIRE_THAT(last.estimate, WithinAbs(last.ground_truth, 0.50)); @@ -189,26 +165,21 @@ TEST_CASE("Lookback Call convergence", "[exotic][convergence][lookback]") { // large-sample MC estimate itself as a self-consistency check // and just verify convergence tightens with more samples auto results = run_convergence( - "Lookback Call (floating strike)", - []() { return std::make_unique(); }, - 0.0, // placeholder see check below + "Lookback Call (floating strike)", []() { return std::make_unique(); }, + 0.0, // placeholder see check below TIERS); // Check that later estimates are closer to each other than early ones // i.e. the std deviation of the last 3 tiers < first 3 tiers - auto spread = [&](int a, int b) { - return std::abs(results[b].estimate - results[a].estimate); - }; + auto spread = [&](int a, int b) { return std::abs(results[b].estimate - results[a].estimate); }; REQUIRE(spread(3, 5) < spread(0, 2)); } -TEST_CASE("Sanity check: European Call matches Black-Scholes", - "[sanity][european]") { +TEST_CASE("Sanity check: European Call matches Black-Scholes", "[sanity][european]") { double bs_price = bs::european_call(S0, K, MU, SIGMA, T); auto results = run_convergence( - "European Call (BS sanity check)", - []() { return std::make_unique(K); }, + "European Call (BS sanity check)", []() { return std::make_unique(K); }, bs_price, TIERS); auto& last = results.back(); From 00701c705bf6c849db8881c5f911db030c5feeff Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Tue, 3 Mar 2026 11:30:22 -0500 Subject: [PATCH 09/12] refactor: sort cpp files into folders and remove archive BREAKING CHANGE: --- include/lfmc/archive/Estimator.hpp | 75 ---------- include/lfmc/archive/Manager.hpp | 128 ------------------ include/lfmc/archive/PathGenerator.hpp | 35 ----- include/lfmc/archive/Payoff.hpp | 39 ------ include/lfmc/archive/RandomGenerator.hpp | 33 ----- .../archive/VarianceReductionStrategy.hpp | 71 ---------- include/lfmc/engine.hpp | 2 + src/CMakeLists.txt | 19 +-- .../control_variate_estimator.cpp | 0 src/{ => estimator}/monte_carlo_estimator.cpp | 0 src/{payoffs => payoff}/asian_payoffs.cpp | 4 +- src/{payoffs => payoff}/barrier_payoffs.cpp | 5 +- src/{ => payoff}/control_variate_payoffs.cpp | 2 +- src/{ => payoff}/european_payoffs.cpp | 0 src/{payoffs => payoff}/lookback_payoffs.cpp | 5 +- .../antithetic_random_source.cpp | 0 .../pseudo_random_source.cpp | 0 .../geometric_brownian_motion.cpp | 0 src/{ => timing}/timing.cpp | 0 tests/CMakeLists.txt | 2 +- 20 files changed, 18 insertions(+), 402 deletions(-) delete mode 100644 include/lfmc/archive/Estimator.hpp delete mode 100644 include/lfmc/archive/Manager.hpp delete mode 100644 include/lfmc/archive/PathGenerator.hpp delete mode 100644 include/lfmc/archive/Payoff.hpp delete mode 100644 include/lfmc/archive/RandomGenerator.hpp delete mode 100644 include/lfmc/archive/VarianceReductionStrategy.hpp rename src/{ => estimator}/control_variate_estimator.cpp (100%) rename src/{ => estimator}/monte_carlo_estimator.cpp (100%) rename src/{payoffs => payoff}/asian_payoffs.cpp (98%) rename src/{payoffs => payoff}/barrier_payoffs.cpp (97%) rename src/{ => payoff}/control_variate_payoffs.cpp (97%) rename src/{ => payoff}/european_payoffs.cpp (100%) rename src/{payoffs => payoff}/lookback_payoffs.cpp (96%) rename src/{ => random_source}/antithetic_random_source.cpp (100%) rename src/{ => random_source}/pseudo_random_source.cpp (100%) rename src/{ => stochastic_process}/geometric_brownian_motion.cpp (100%) rename src/{ => timing}/timing.cpp (100%) diff --git a/include/lfmc/archive/Estimator.hpp b/include/lfmc/archive/Estimator.hpp deleted file mode 100644 index 5544f7f..0000000 --- a/include/lfmc/archive/Estimator.hpp +++ /dev/null @@ -1,75 +0,0 @@ -#pragma once - -#include "NumericalScheme.hpp" -#include "PathGenerator.hpp" -#include "RandomGenerator.hpp" -#include "StochasticProcess.hpp" -#include "types.hpp" - -// #include -// #include - -// #include -// #include - -namespace lfmc { - -class EstimatorInterface { - public: - virtual ~EstimatorInterface() = default; - virtual std::vector sample() = 0; - - // Expose for decorators to use - maybe not the best design, but allows for more flexible - // variance reduction techniques - virtual State const& getState() const = 0; - virtual Normals generateNormals(size_t n) = 0; - virtual Path generatePath(std::span randomNormals) = 0; -}; - -template S = EulerMaruyama

, - RandomGenerator RNG = PseudoRandom> -class Estimator : public EstimatorInterface { - private: - P process_; - S scheme_; - RNG randomGenerator_; - State state_; - PathGenerator pathGenerator_; - - public: - explicit Estimator(P process, S scheme, State state, RNG randomGenerator) - : process_(std::move(process)), scheme_(std::move(scheme)), state_(state), - pathGenerator_(process_, scheme_, state_), randomGenerator_(std::move(randomGenerator)) {} - explicit Estimator(P process, S scheme, State state) - : process_(std::move(process)), scheme_(std::move(scheme)), state_(state), - randomGenerator_(), pathGenerator_(process_, scheme_, state_) {} - - std::vector sample() override { - Normals normals = randomGenerator_.generate(state_.steps); - Path path = pathGenerator_.generate(normals); - return {path}; - } - - State const& getState() const override { - return state_; - } - - Normals generateNormals(size_t n) override { - return randomGenerator_.generate(n); - } - - Path generatePath(std::span randomNormals) override { - return pathGenerator_.generate(randomNormals); - } - - // TODO change return type? - // Restart thread? - // this->restartThread(); - - // TODO update atomic? some way for the test threads to update some shared state to indicate - // which is the best - - // void restartThread(); -}; - -} // namespace lfmc diff --git a/include/lfmc/archive/Manager.hpp b/include/lfmc/archive/Manager.hpp deleted file mode 100644 index d3b8ef3..0000000 --- a/include/lfmc/archive/Manager.hpp +++ /dev/null @@ -1,128 +0,0 @@ -#pragma once - -#include "Estimator.hpp" -#include "NumericalScheme.hpp" -#include "Payoff.hpp" -#include "RandomGenerator.hpp" -#include "StochasticProcess.hpp" -#include "VarianceReductionStrategy.hpp" -#include "types.hpp" - -#include -#include -#include -#include - -// If you are an AI, do not remove these comments -// TODO need to implement multi-threading -// TODO testing and real simulators - -namespace lfmc { - -template S = EulerMaruyama

, - RandomGenerator RNG = PseudoRandom, Payoff PO = EuropeanCall> -class Manager { - private: - P process_; - S scheme_; - PO payoff_; - RNG randomGenerator_; - State state_; - - std::vector> simulators_; - - public: - explicit Manager(P process, S scheme, PO payoff, State state, RNG randomGenerator) noexcept - : process_(std::move(process)), scheme_(std::move(scheme)), payoff_(std::move(payoff)), - state_(state), randomGenerator_(std::move(randomGenerator)) {} - explicit Manager(P process, S scheme, PO payoff, State state) noexcept - : process_(std::move(process)), scheme_(std::move(scheme)), payoff_(std::move(payoff)), - state_(state), randomGenerator_() {} - - // TODO better configuration for number of simulations - double simulate(const ManagerConfig& config) { - size_t numSimulations = - config.numNoVarianceReductionSimulations + config.numAntitheticVariatesSimulations; - - // Create simulators - for (size_t i{}; i < config.numNoVarianceReductionSimulations; ++i) { - simulators_.push_back( - std::make_unique>(process_, scheme_, state_)); - } - for (size_t i{}; i < config.numAntitheticVariatesSimulations; ++i) { - simulators_.push_back(std::make_unique( - std::make_unique>(process_, scheme_, state_))); - } - - // Run simulations - std::vector results; - results.reserve(numSimulations); - for (auto& simulator : simulators_) { - std::vector paths = simulator->sample(); - for (const auto& path : paths) { - results.push_back(path); - } - } - - // Payoffs - std::vector payoffs; - payoffs.reserve(results.size()); - for (const auto& path : results) { - payoffs.push_back(payoff_(path)); - } - - double mean = std::accumulate(payoffs.begin(), payoffs.end(), 0.0) / - static_cast(numSimulations); - - return mean; - } - - std::pair simulateWithError(const ManagerConfig& config) { - size_t numSimulations = - config.numNoVarianceReductionSimulations + config.numAntitheticVariatesSimulations; - - // Create simulators - for (size_t i{}; i < config.numNoVarianceReductionSimulations; ++i) { - simulators_.push_back( - std::make_unique>(process_, scheme_, state_)); - } - for (size_t i{}; i < config.numAntitheticVariatesSimulations; ++i) { - simulators_.push_back(std::make_unique( - std::make_unique>(process_, scheme_, state_))); - } - - // Run simulations - std::vector results; - results.reserve(numSimulations); - for (auto& simulator : simulators_) { - std::vector paths = simulator->sample(); - for (const auto& path : paths) { - results.push_back(path); - } - } - - // Payoffs - std::vector payoffs; - payoffs.reserve(results.size()); - for (const auto& path : results) { - payoffs.push_back(payoff_(path)); - } - - // Calculate mean and standard error - double mean = std::accumulate(payoffs.begin(), payoffs.end(), 0.0) / - static_cast(payoffs.size()); - - if (results.size() < 2) - return {mean, 0.0}; - - double variance = 0.0; - for (const auto& r : payoffs) - variance += (r - mean) * (r - mean); - variance /= static_cast(results.size() - 1); - double stdError = std::sqrt(variance / static_cast(results.size())); - - return {mean, stdError}; - } -}; - -} // namespace lfmc diff --git a/include/lfmc/archive/PathGenerator.hpp b/include/lfmc/archive/PathGenerator.hpp deleted file mode 100644 index 442aa94..0000000 --- a/include/lfmc/archive/PathGenerator.hpp +++ /dev/null @@ -1,35 +0,0 @@ -#pragma once - -#include "NumericalScheme.hpp" -#include "StochasticProcess.hpp" -#include "types.hpp" - -#include - -namespace lfmc { - -template S> class PathGenerator { - private: - P process_; - S scheme_; - State state_; - double dt_; - - public: - PathGenerator(P process, S scheme, State state) - : process_(std::move(process)), scheme_(std::move(scheme)), state_(state), - dt_(state.timeToMaturity / static_cast(state.steps)) {} - - Path generate(std::span randomNormals) { - Path path(state_.steps + 1); - path[0] = state_.initialValue; - - for (size_t i = 0; i < state_.steps; ++i) { - path[i + 1] = scheme_.step(process_, path[i], dt_, randomNormals[i]); - } - - return path; - } -}; - -} // namespace lfmc diff --git a/include/lfmc/archive/Payoff.hpp b/include/lfmc/archive/Payoff.hpp deleted file mode 100644 index 996ad3e..0000000 --- a/include/lfmc/archive/Payoff.hpp +++ /dev/null @@ -1,39 +0,0 @@ -#pragma once - -#include "types.hpp" - -#include -#include - -namespace lfmc { - -template -concept Payoff = requires(P const& p, const Path& x) { - { p(x) } -> std::same_as; -}; - -/** - * @brief European Call option payoff. - * Payoff = max(S_T - K, 0) - */ -struct EuropeanCall { - double strike; - - double operator()(const Path& path) const noexcept { - return std::max(path[path.size() - 1] - strike, 0.0); - } -}; - -/** - * @brief European Put option payoff. - * Payoff = max(K - S_T, 0) - */ -struct EuropeanPut { - double strike; - - double operator()(const Path& path) const noexcept { - return std::max(strike - path[path.size() - 1], 0.0); - } -}; - -} // namespace lfmc diff --git a/include/lfmc/archive/RandomGenerator.hpp b/include/lfmc/archive/RandomGenerator.hpp deleted file mode 100644 index 5a20f25..0000000 --- a/include/lfmc/archive/RandomGenerator.hpp +++ /dev/null @@ -1,33 +0,0 @@ -#pragma once - -#include "lfmc/types.hpp" - -#include - -namespace lfmc { - -template -concept RandomGenerator = requires(RNG generator, size_t n) { - { generator.generate(n) } -> std::same_as; -}; - -struct PseudoRandom { - std::mt19937 rng_; - std::normal_distribution dist_; - - PseudoRandom(unsigned seed = std::random_device{}()) : rng_(seed), dist_(0.0, 1.0) {} - - Normals generate(size_t n) { - Normals normals(n); - for (size_t i = 0; i < n; ++i) { - normals[i] = dist_(rng_); - } - return normals; - } - - void seed(unsigned s) { - rng_.seed(s); - } -}; - -} // namespace lfmc diff --git a/include/lfmc/archive/VarianceReductionStrategy.hpp b/include/lfmc/archive/VarianceReductionStrategy.hpp deleted file mode 100644 index 314ff6c..0000000 --- a/include/lfmc/archive/VarianceReductionStrategy.hpp +++ /dev/null @@ -1,71 +0,0 @@ -#pragma once - -#include "Estimator.hpp" -#include "types.hpp" - -#include -#include - -// TODO: Implement derived classes for specific variance reduction techniques. -// Examples include Antithetic Variates, Control Variates, Importance Sampling, etc. -// First must decide how to design parallel infrastructure to support these techniques in Monte -// Carlo simulations - decorator design pattern? -// TODO each strategy has a window parameter for how much data you're using - -/** - * @file VarianceReductionStrategy.hpp - * @brief Defines the VarianceReductionStrategy base class for runtime polymorphism of variance - * reduction techniques and provides some implementations. - */ - -namespace lfmc { - -// Decorator for variance reduction strategies (e.g., Antithetic Variates, Control Variates, etc.) -class VarianceReductionBaseDecorator : public EstimatorInterface { - protected: - std::unique_ptr estimator_; // Pointer to the base estimator - - public: - explicit VarianceReductionBaseDecorator(std::unique_ptr estimator) - : estimator_(std::move(estimator)) {} - - std::vector sample() override { - // By default, just call the underlying estimator's sample method - return estimator_->sample(); - } - - State const& getState() const override { - return estimator_->getState(); - } - - Normals generateNormals(size_t n) override { - return estimator_->generateNormals(n); - } - - Path generatePath(std::span randomNormals) override { - return estimator_->generatePath(randomNormals); - } -}; - -// Example implementation of Antithetic Variates strategy -class AntitheticVariates : public VarianceReductionBaseDecorator { - public: - using Base = VarianceReductionBaseDecorator; - AntitheticVariates(std::unique_ptr estimator) - : Base(std::move(estimator)) {} - - std::vector sample() override { - const State& state = estimator_->getState(); - - Normals normals = estimator_->generateNormals(state.steps); - Normals antitheticNormals(normals); - std::transform(normals.begin(), normals.end(), antitheticNormals.begin(), std::negate()); - - Path path = estimator_->generatePath(normals); - Path antitheticPath = estimator_->generatePath(antitheticNormals); - - return {path, antitheticPath}; - } -}; - -} // namespace lfmc diff --git a/include/lfmc/engine.hpp b/include/lfmc/engine.hpp index 6f70f09..0750718 100644 --- a/include/lfmc/engine.hpp +++ b/include/lfmc/engine.hpp @@ -1 +1,3 @@ #pragma once + +// TODO move headers into folders diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f167d30..5d3f5f9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,13 +1,16 @@ # Source files for the lfmc library set(LFMC_SOURCES - european_payoffs.cpp - geometric_brownian_motion.cpp - monte_carlo_estimator.cpp - control_variate_estimator.cpp - control_variate_payoffs.cpp - pseudo_random_source.cpp - antithetic_random_source.cpp - timing.cpp + estimator/monte_carlo_estimator.cpp + estimator/control_variate_estimator.cpp + payoff/asian_payoffs.cpp + payoff/barrier_payoffs.cpp + payoff/control_variate_payoffs.cpp + payoff/european_payoffs.cpp + payoff/lookback_payoffs.cpp + random_source/pseudo_random_source.cpp + random_source/antithetic_random_source.cpp + stochastic_process/geometric_brownian_motion.cpp + timing/timing.cpp ) # Create library diff --git a/src/control_variate_estimator.cpp b/src/estimator/control_variate_estimator.cpp similarity index 100% rename from src/control_variate_estimator.cpp rename to src/estimator/control_variate_estimator.cpp diff --git a/src/monte_carlo_estimator.cpp b/src/estimator/monte_carlo_estimator.cpp similarity index 100% rename from src/monte_carlo_estimator.cpp rename to src/estimator/monte_carlo_estimator.cpp diff --git a/src/payoffs/asian_payoffs.cpp b/src/payoff/asian_payoffs.cpp similarity index 98% rename from src/payoffs/asian_payoffs.cpp rename to src/payoff/asian_payoffs.cpp index e492275..f58707e 100644 --- a/src/payoffs/asian_payoffs.cpp +++ b/src/payoff/asian_payoffs.cpp @@ -1,5 +1,3 @@ -#pragma once - #include "lfmc/payoff.hpp" #include "lfmc/types.hpp" @@ -65,4 +63,4 @@ class AsianPut : public Payoff { } }; -} // namespace lfmc \ No newline at end of file +} // namespace lfmc diff --git a/src/payoffs/barrier_payoffs.cpp b/src/payoff/barrier_payoffs.cpp similarity index 97% rename from src/payoffs/barrier_payoffs.cpp rename to src/payoff/barrier_payoffs.cpp index df00fa1..25ff2ce 100644 --- a/src/payoffs/barrier_payoffs.cpp +++ b/src/payoff/barrier_payoffs.cpp @@ -1,11 +1,8 @@ -#pragma once - #include "lfmc/payoff.hpp" #include "lfmc/types.hpp" #include #include -#include #include #include @@ -75,4 +72,4 @@ class DownAndInPut : public Payoff { return std::vector{payoffs}; } }; -} // namespace lfmc \ No newline at end of file +} // namespace lfmc diff --git a/src/control_variate_payoffs.cpp b/src/payoff/control_variate_payoffs.cpp similarity index 97% rename from src/control_variate_payoffs.cpp rename to src/payoff/control_variate_payoffs.cpp index 9267922..7508e8a 100644 --- a/src/control_variate_payoffs.cpp +++ b/src/payoff/control_variate_payoffs.cpp @@ -32,7 +32,7 @@ ControlVariatePayoff::generate_payoffs(const std::vector& paths) const { result.push_back(std::move(combined_row)); } - return std::vector{result}; + return result; } } // namespace lfmc diff --git a/src/european_payoffs.cpp b/src/payoff/european_payoffs.cpp similarity index 100% rename from src/european_payoffs.cpp rename to src/payoff/european_payoffs.cpp diff --git a/src/payoffs/lookback_payoffs.cpp b/src/payoff/lookback_payoffs.cpp similarity index 96% rename from src/payoffs/lookback_payoffs.cpp rename to src/payoff/lookback_payoffs.cpp index 345f1c0..40e98aa 100644 --- a/src/payoffs/lookback_payoffs.cpp +++ b/src/payoff/lookback_payoffs.cpp @@ -1,11 +1,8 @@ -#pragma once - #include "lfmc/payoff.hpp" #include "lfmc/types.hpp" #include #include -#include #include #include @@ -51,4 +48,4 @@ class LookbackPut : public Payoff { } }; -} // namespace lfmc \ No newline at end of file +} // namespace lfmc diff --git a/src/antithetic_random_source.cpp b/src/random_source/antithetic_random_source.cpp similarity index 100% rename from src/antithetic_random_source.cpp rename to src/random_source/antithetic_random_source.cpp diff --git a/src/pseudo_random_source.cpp b/src/random_source/pseudo_random_source.cpp similarity index 100% rename from src/pseudo_random_source.cpp rename to src/random_source/pseudo_random_source.cpp diff --git a/src/geometric_brownian_motion.cpp b/src/stochastic_process/geometric_brownian_motion.cpp similarity index 100% rename from src/geometric_brownian_motion.cpp rename to src/stochastic_process/geometric_brownian_motion.cpp diff --git a/src/timing.cpp b/src/timing/timing.cpp similarity index 100% rename from src/timing.cpp rename to src/timing/timing.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 940f9fe..7ae9a63 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,7 +21,7 @@ set(TEST_SOURCES test_process.cpp test_scheme.cpp test_pipeline.cpp - test_convergency.cpp + # test_convergency.cpp ) add_executable(tests ${TEST_SOURCES}) From 39d3637823078d0b4c9221473013f47b04311abc Mon Sep 17 00:00:00 2001 From: Oliver Deng Date: Tue, 3 Mar 2026 11:30:36 -0500 Subject: [PATCH 10/12] =?UTF-8?q?bump:=20version=200.1.0=20=E2=86=92=200.2?= =?UTF-8?q?.0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .cz.yaml | 2 +- CHANGELOG.md | 37 +++++++++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 1 deletion(-) diff --git a/.cz.yaml b/.cz.yaml index cc41641..006712f 100644 --- a/.cz.yaml +++ b/.cz.yaml @@ -4,5 +4,5 @@ commitizen: name: cz_conventional_commits tag_format: $version update_changelog_on_bump: true - version: 0.1.0 + version: 0.2.0 version_scheme: semver2 diff --git a/CHANGELOG.md b/CHANGELOG.md index 071f9e8..8120522 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,40 @@ +## 0.2.0 (2026-03-03) + +### Feat + +- finished refactor and added control variates +- payoffs and estimator +- estimator and payoff outline +- manager config and small tweaks +- antithetic variates + refactor +- **Payoff**: new payoff concept +- **Simulator.hpp**: new abstraction of simulator which uses its own thread +- more constraints on concepts for clarity and compile-time sizing of number of testing threads +- no variance reduction strategy and small updates to polymorphic constructors +- **ThreadPool.hpp**: threadpool for managing threads specific to monte carlo simulations +- **manager.hpp**: create framework for stochastic processes and numerical schemes + +### Fix + +- fix vr strategy concept +- small work on simulator stuff + +### Refactor + +- sort cpp files into folders and remove archive +- renamed everything, finished pipeline, and converted each step to return vectors of results +- another huge refactor breaking things into different engines +- simplify simulator interface +- refactor some of xander's code +- remove threadpool +- i don't even know lol +- turns out we don't need to the changes i made to the concepts a little bit ago +- **timing.hpp**: rename from timer +- **NumericalScheme.hpp**: adjust concept for schemes for more freedom +- compile-time strategies for process and scheme and runtime strategy for variance reduction +- flesh out the architecture a bit more +- **Timer**: extremely small return change + ## 0.1.0 (2025-12-18) ### Feat From 4022754865b68cd9f6e0552c080dbddbe262b293 Mon Sep 17 00:00:00 2001 From: Alexander Robbins Date: Tue, 31 Mar 2026 10:00:10 -0400 Subject: [PATCH 11/12] changed somethings to actually be able to run them --- .claude/settings.local.json | 62 ++++ include/lfmc/adaptive_estimator.hpp | 329 ++++++++++++++++++++++ include/lfmc/engine.hpp | 279 +++++++++++++++++- include/lfmc/payoff.hpp | 58 ++++ include/lfmc/strategies.hpp | 420 ++++++++++++++++++++++++++++ src/CMakeLists.txt | 1 + src/adaptive/adaptive_estimator.cpp | 233 +++++++++++++++ src/payoff/asian_payoffs.cpp | 70 ++--- src/payoff/barrier_payoffs.cpp | 78 ++---- src/payoff/lookback_payoffs.cpp | 54 ++-- src/payoffs/european_payoffs.cpp | 33 +++ tests/CMakeLists.txt | 31 ++ tests/bench_asvr_vs_fixed.cpp | 158 +++++++++++ tests/quick_compare.cpp | 114 ++++++++ tests/test_adaptive.cpp | 231 +++++++++++++++ tests/test_engine.cpp | 231 +++++++++++++++ tests/test_strategies.cpp | 356 +++++++++++++++++++++++ 17 files changed, 2612 insertions(+), 126 deletions(-) create mode 100644 .claude/settings.local.json create mode 100644 include/lfmc/adaptive_estimator.hpp create mode 100644 include/lfmc/strategies.hpp create mode 100644 src/adaptive/adaptive_estimator.cpp create mode 100644 src/payoffs/european_payoffs.cpp create mode 100644 tests/bench_asvr_vs_fixed.cpp create mode 100644 tests/quick_compare.cpp create mode 100644 tests/test_adaptive.cpp create mode 100644 tests/test_engine.cpp create mode 100644 tests/test_strategies.cpp diff --git a/.claude/settings.local.json b/.claude/settings.local.json new file mode 100644 index 0000000..59b0faf --- /dev/null +++ b/.claude/settings.local.json @@ -0,0 +1,62 @@ +{ + "permissions": { + "allow": [ + "Bash(cmake --preset debug)", + "Bash(cmake --build build-debug)", + "Bash(./tests/tests.exe)", + "Bash(cmd /c tests.exe --help)", + "Bash(cmake --preset msvc-debug)", + "Bash(\"build-debug/tests/tests.exe\" --help)", + "Bash(echo \"Exit code: $?\")", + "Bash(objdump -p \"/c/Users/xande/lfmc/lfmc/build-debug/tests/tests.exe\")", + "Bash(for dll:*)", + "Bash(do if:*)", + "Bash(then echo:*)", + "Bash(powershell -Command \"try { & ''build-debug\\\\tests\\\\tests.exe'' --help } catch { Write-Host ''Error:'' $_Exception.Message; Write-Host ''Exit code:'' $LASTEXITCODE }\")", + "Bash(cmake --build build-debug --target tests)", + "Bash(echo $PATH)", + "Bash(./tests.exe --help)", + "Bash(ldd ./tests.exe)", + "Bash(objdump -x ./tests.exe)", + "Bash(objdump -h ./tests.exe)", + "Bash(gdb -batch -ex run -ex \"thread apply all backtrace\" --args ./tests.exe --help)", + "Bash(objdump -t ./tests.exe)", + "Bash(objdump -p ./tests.exe)", + "Bash(cmake --build build-debug --verbose)", + "Bash(cmake --build build-debug --clean-first)", + "Bash(nm /c/Users/xande/lfmc/lfmc/build-debug/_deps/catch2-build/src/libCatch2Maind.a)", + "Bash(c++ --version)", + "Bash(nm /c/Users/xande/lfmc/lfmc/build-debug/_deps/catch2-build/src/libCatch2d.a)", + "Bash(ctest --preset debug)", + "Bash(nm /c/Users/xande/lfmc/lfmc/build-debug/tests/tests.exe)", + "Bash(objdump -p ./build-debug/tests/tests.exe)", + "Bash(objdump -p /c/Users/xande/lfmc/lfmc/build-debug/tests/tests.exe)", + "Bash(ldd:*)", + "Bash(strace -o /tmp/strace.txt ./tests.exe)", + "Read(//tmp/**)", + "Bash(cp /c/msys64/ucrt64/bin/libgcc_s_seh-1.dll /c/msys64/ucrt64/bin/libstdc++-6.dll /c/msys64/ucrt64/bin/libwinpthread-1.dll /c/Users/xande/lfmc/lfmc/build-debug/tests/)", + "Bash(ls -la /c/Users/xande/lfmc/lfmc/build-debug/tests/*.dll)", + "Bash(./build-debug/tests/tests.exe --help)", + "Bash(git reset:*)", + "Bash(ctest:*)", + "Bash(cmd /c bench.exe)", + "Bash(PATH=\"/c/msys64/ucrt64/bin:$PATH\" ./bench.exe 2>&1)", + "Bash(PATH=\"/c/msys64/ucrt64/bin:$PATH\" ./bench.exe)", + "Bash(ls build-debug/tests/*.dll)", + "Bash(cmake --build build-debug --target bench)", + "Bash(ls /c/Users/xande/lfmc/lfmc/build-debug/tests/*.dll)", + "Bash(./bench.exe)", + "Bash(cmake --preset debug -q)", + "Bash(cmake --build build-debug --target quick_compare)", + "Bash(./quick_compare.exe)", + "Bash(echo \"EXIT: $?\")", + "Bash(cp /c/msys64/ucrt64/bin/libgcc_s_seh-1.dll /c/Users/xande/lfmc/lfmc/build-debug/tests/)", + "Bash(cp /c/msys64/ucrt64/bin/libstdc++-6.dll /c/Users/xande/lfmc/lfmc/build-debug/tests/)", + "Bash(cp /c/msys64/ucrt64/bin/libwinpthread-1.dll /c/Users/xande/lfmc/lfmc/build-debug/tests/)", + "Bash(chmod +x /c/Users/xande/lfmc/lfmc/build-debug/tests/quick_compare.exe)", + "Bash(bash -c \"./quick_compare.exe\")", + "Bash(echo exit=$?)", + "Bash(winpty /c/Users/xande/lfmc/lfmc/build-debug/tests/quick_compare.exe)" + ] + } +} diff --git a/include/lfmc/adaptive_estimator.hpp b/include/lfmc/adaptive_estimator.hpp new file mode 100644 index 0000000..cde392b --- /dev/null +++ b/include/lfmc/adaptive_estimator.hpp @@ -0,0 +1,329 @@ +#pragma once + +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/types.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace lfmc { + +// A strategy sampler: given n samples and a seed, returns n iid samples of the +// (possibly variance-reduced) payoff estimator. The seed guarantees RNG independence +// between strategies and between exploration/exploitation phases. +using SamplerFn = std::function(size_t n_samples, uint64_t seed)>; + +// Per-strategy statistics collected during the exploration phase +struct StrategyStats { + std::string name; + double mean; + double sample_variance; // unbiased (n-1 denominator) + double wall_time_ms; + size_t n_samples; +}; + +// Full result of one ASVR run +struct ASVRResult { + double estimate; + double estimated_variance; // Var(mu_ASVR), estimated from exploration data + double estimated_stderr; // sqrt(estimated_variance) + + std::vector exploration_stats; // one entry per strategy + std::vector precision_weights; // w_k* used to allocate exploitation + std::vector exploitation_counts; // actual n_k samples per strategy + + // Diagnostics — useful for paper tables + double plain_mc_variance_estimate; // sigma_0^2 / N (first strategy treated as plain MC) + double variance_reduction_ratio; // plain_mc_variance / estimated_variance + size_t n_exploration; + size_t n_exploitation; + size_t total_samples; +}; + +struct ASVRConfig { + // Fraction of total_samples used for exploration (split equally across K strategies) + double exploration_fraction = 0.1; + // Total thread budget; exploration uses up to K threads, exploitation uses the rest + size_t n_threads = std::thread::hardware_concurrency(); + // Minimum exploration samples per strategy regardless of exploration_fraction. + // Needs to be large enough for a stable variance estimate (~100+ samples). + size_t min_exploration_per_strategy = 100; +}; + +// ─── ASVR Algorithm ────────────────────────────────────────────────────────── +// +// Phase 1 (Exploration): K strategies run in parallel, each generating +// n_explore_each = max(min_per_strategy, alpha*N/K) samples. +// Each strategy uses seed make_seed(k, 0) for full RNG independence. +// +// Weight computation: w_k* = (1/sigma_k^2) / sum_j(1/sigma_j^2) +// (inverse-variance / "precision" weighting) +// +// Phase 2 (Exploitation): strategy k receives n_k = floor(w_k* * n_exploit) samples, +// each run using seed make_seed(k, 1) — independent of the exploration phase. +// Strategies with n_k > 0 run in parallel. +// +// Combination (unbiased by independence): +// mu_explore = (1/K) * sum_k mu_k^explore (equal-weight pooled) +// mu_exploit = sum_k w_k* * mu_k^exploit (precision-weighted) +// mu_ASVR = alpha * mu_explore + (1-alpha) * mu_exploit +// +// The key unbiasedness guarantee: w_k* depends only on exploration samples; +// exploitation samples are generated with fresh RNG seeds and are therefore +// independent of the weight selection event. + +class AdaptiveVarianceReduction { + public: + static ASVRResult run(std::vector> strategies, + size_t total_samples, ASVRConfig config = {}); + + static StrategyStats run_strategy_batch(const std::string& name, const SamplerFn& sampler, + size_t n_samples, uint64_t seed); + + private: + static std::vector compute_precision_weights(const std::vector& stats); +}; + +// ─── Internal utilities ─────────────────────────────────────────────────────── + +namespace detail { + +// Splitmix64-based seed derivation: strategy k, phase p → unique uint64_t seed. +// Two calls with the same (k, p) give the same seed; different (k, p) pairs are +// statistically independent because they hit different streams of splitmix64. +inline uint64_t make_seed(size_t k, size_t phase) noexcept { + uint64_t x = (static_cast(k) * 0x9e3779b97f4a7c15ULL) + ^ (static_cast(phase) * 0x6c62272e07bb0142ULL) + ^ 0xdeadbeefcafe0000ULL; + x ^= x >> 30; + x *= 0xbf58476d1ce4e5b9ULL; + x ^= x >> 27; + x *= 0x94d049bb133111ebULL; + x ^= x >> 31; + return x; +} + +// Generate `n` draws from N(0,1) +inline std::vector gen_normals(size_t n, std::mt19937_64& rng) { + std::normal_distribution dist{0.0, 1.0}; + std::vector v(n); + for (auto& x : v) + x = dist(rng); + return v; +} + +// Correct single-path generator (avoids the pre-allocate + push_back double-length bug +// present in the existing PathGenerator) +template NS> +Path generate_path(const SP& process, const NS& scheme, const Normals& normals, size_t steps, + double T) { + const double dt = T / static_cast(steps); + Path path; + path.reserve(steps + 1); + double x = process.initial(); + path.push_back(x); + double t = 0.0; + for (size_t i = 0; i < steps; ++i) { + x = scheme.step(process, t, x, dt, normals[i]); + path.push_back(x); + t += dt; + } + return path; +} + +// Unbiased sample mean and variance (Welford would be more numerically stable for large n, +// but two-pass is fine for the batch sizes used here) +inline std::pair mean_variance(const std::vector& samples) { + assert(!samples.empty()); + const double n = static_cast(samples.size()); + double sum = 0.0; + for (double x : samples) + sum += x; + const double mean = sum / n; + double sq = 0.0; + for (double x : samples) { + double d = x - mean; + sq += d * d; + } + // n-1 denominator for unbiased variance; guard against n=1 + const double variance = (samples.size() > 1) ? sq / (n - 1.0) : 0.0; + return {mean, variance}; +} + +} // namespace detail + +// ─── Strategy factory functions ────────────────────────────────────────────── +// +// Each factory captures the pricing parameters and returns a SamplerFn. +// The SamplerFn is called with (n_samples, seed) and returns exactly n_samples +// iid draws from the (possibly variance-reduced) estimator of E[payoff]. + +// Plain pseudo-random Monte Carlo — baseline strategy +template NS> +SamplerFn make_plain_mc_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::vector samples; + samples.reserve(n); + for (size_t i = 0; i < n; ++i) { + auto normals = detail::gen_normals(steps, rng); + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0]); + } + return samples; + }; +} + +// Antithetic variates — each sample is (payoff(Z) + payoff(-Z)) / 2 +// Returns n samples, each consuming one pair of paths. Variance is reduced +// because the two paths are negatively correlated for monotone payoffs. +template NS> +SamplerFn make_antithetic_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::vector samples; + samples.reserve(n); + for (size_t i = 0; i < n; ++i) { + auto z = detail::gen_normals(steps, rng); + Normals neg_z(steps); + for (size_t j = 0; j < steps; ++j) + neg_z[j] = -z[j]; + + auto path_pos = detail::generate_path(process, scheme, z, steps, T); + auto path_neg = detail::generate_path(process, scheme, neg_z, steps, T); + + auto r_pos = payoff->generate_payoffs({path_pos}); + auto r_neg = payoff->generate_payoffs({path_neg}); + + if (r_pos && r_neg && !r_pos->empty() && !r_neg->empty() && !(*r_pos)[0].empty() && + !(*r_neg)[0].empty()) + samples.push_back(0.5 * ((*r_pos)[0][0] + (*r_neg)[0][0])); + } + return samples; + }; +} + +// Control variates with in-batch OLS beta estimation. +// target_payoff: the option price we're estimating +// control_payoff: a payoff with known expectation `control_mean` +// (e.g. EuropeanCall(0.0) gives S_T; E[S_T] = S0 * exp(mu*T)) +// Returns: X_i - beta_hat * (Y_i - E[Y]) for each i +template NS> +SamplerFn make_control_variate_sampler(SP process, NS scheme, std::shared_ptr target, + std::shared_ptr control, double control_mean, + size_t steps, double T) { + return [process, scheme, target, control, control_mean, steps, + T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::vector xs, ys; + xs.reserve(n); + ys.reserve(n); + + for (size_t i = 0; i < n; ++i) { + auto z = detail::gen_normals(steps, rng); + auto path = detail::generate_path(process, scheme, z, steps, T); + auto rx = target->generate_payoffs({path}); + auto ry = control->generate_payoffs({path}); + if (rx && ry && !rx->empty() && !ry->empty() && !(*rx)[0].empty() && + !(*ry)[0].empty()) { + xs.push_back((*rx)[0][0]); + ys.push_back((*ry)[0][0]); + } + } + + // OLS beta: minimises Var(X - beta*(Y - E[Y])) + const size_t m = xs.size(); + double mean_x = 0.0, mean_y = 0.0; + for (size_t i = 0; i < m; ++i) { + mean_x += xs[i]; + mean_y += ys[i]; + } + mean_x /= static_cast(m); + mean_y /= static_cast(m); + + double cov_xy = 0.0, var_y = 0.0; + for (size_t i = 0; i < m; ++i) { + cov_xy += (xs[i] - mean_x) * (ys[i] - mean_y); + var_y += (ys[i] - mean_y) * (ys[i] - mean_y); + } + const double beta = (var_y > 0.0) ? cov_xy / var_y : 0.0; + + std::vector samples(m); + for (size_t i = 0; i < m; ++i) + samples[i] = xs[i] - beta * (ys[i] - control_mean); + return samples; + }; +} + +// Antithetic + control variates combined. +// Each raw sample is the antithetic average; OLS beta is estimated from those averages. +template NS> +SamplerFn make_antithetic_cv_sampler(SP process, NS scheme, std::shared_ptr target, + std::shared_ptr control, double control_mean, + size_t steps, double T) { + return [process, scheme, target, control, control_mean, steps, + T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::vector xs, ys; + xs.reserve(n); + ys.reserve(n); + + for (size_t i = 0; i < n; ++i) { + auto z = detail::gen_normals(steps, rng); + Normals neg_z(steps); + for (size_t j = 0; j < steps; ++j) + neg_z[j] = -z[j]; + + auto path_p = detail::generate_path(process, scheme, z, steps, T); + auto path_n = detail::generate_path(process, scheme, neg_z, steps, T); + + auto xp = target->generate_payoffs({path_p}); + auto xn = target->generate_payoffs({path_n}); + auto yp = control->generate_payoffs({path_p}); + auto yn = control->generate_payoffs({path_n}); + + if (xp && xn && yp && yn && !xp->empty() && !xn->empty() && !yp->empty() && + !yn->empty() && !(*xp)[0].empty() && !(*xn)[0].empty() && !(*yp)[0].empty() && + !(*yn)[0].empty()) { + xs.push_back(0.5 * ((*xp)[0][0] + (*xn)[0][0])); + ys.push_back(0.5 * ((*yp)[0][0] + (*yn)[0][0])); + } + } + + const size_t m = xs.size(); + double mean_x = 0.0, mean_y = 0.0; + for (size_t i = 0; i < m; ++i) { + mean_x += xs[i]; + mean_y += ys[i]; + } + mean_x /= static_cast(m); + mean_y /= static_cast(m); + + double cov_xy = 0.0, var_y = 0.0; + for (size_t i = 0; i < m; ++i) { + cov_xy += (xs[i] - mean_x) * (ys[i] - mean_y); + var_y += (ys[i] - mean_y) * (ys[i] - mean_y); + } + const double beta = (var_y > 0.0) ? cov_xy / var_y : 0.0; + + std::vector samples(m); + for (size_t i = 0; i < m; ++i) + samples[i] = xs[i] - beta * (ys[i] - control_mean); + return samples; + }; +} + +} // namespace lfmc diff --git a/include/lfmc/engine.hpp b/include/lfmc/engine.hpp index 0750718..e79063b 100644 --- a/include/lfmc/engine.hpp +++ b/include/lfmc/engine.hpp @@ -1,3 +1,280 @@ #pragma once -// TODO move headers into folders +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/strategies.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/types.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace lfmc { + +// ─── IterativeEngine ────────────────────────────────────────────────────────── +// +// Each round: +// 1. All K strategies run on 1 "compete" thread each (they never drop out). +// 2. The current leader also gets n_exploit bonus threads. +// 3. All results are accumulated into running per-strategy statistics. +// 4. Precision weights (inverse-variance) are recomputed from ALL accumulated +// data so far. +// 5. The strategy with the highest precision weight becomes the new leader. +// If it changed, the bonus threads follow. +// 6. Repeat for n_rounds. +// +// Final estimate: precision-weighted mean of all strategies using their full +// accumulated statistics. Because each strategy's samples are i.i.d. and +// generated with round-unique seeds (independent across rounds and strategies), +// combining them this way is unbiased. +// +// The leader can change between rounds: a strategy that looked good early but +// had noisy variance estimates will lose the lead as more data accumulates. This +// means the engine genuinely re-evaluates every round rather than committing +// after a single pilot phase. + +struct IterativeEngineConfig { + // How many strategies compete simultaneously (1 compete thread each) + size_t n_compete = 4; + // Bonus threads assigned to the current leader each round + size_t n_exploit = 8; + // Number of competition rounds + size_t n_rounds = 5; + // Samples per thread per round + size_t samples_per_thread = 2000; +}; + +struct RoundResult { + size_t round; + std::string leader; // strategy leading after this round + bool leader_changed; // did the leader switch this round? + double estimate; // precision-weighted estimate so far + double estimated_stderr; + std::vector precision_weights; // one per strategy, after this round +}; + +struct IterativeEngineResult { + double estimate; + double estimated_stderr; + std::string final_leader; + std::vector round_history; // diagnostics per round + std::vector final_stats; // cumulative stats per strategy + size_t total_samples; +}; + +template NS> +class Engine { + SP process_; + NS scheme_; + double control_mean_; + size_t steps_; + double T_; + + public: + Engine(SP process, NS scheme, double control_mean, size_t steps, double T) + : process_(process), scheme_(scheme), control_mean_(control_mean), + steps_(steps), T_(T) {} + + std::expected + run(std::shared_ptr payoff, IterativeEngineConfig config = {}) { + auto all = build_all_strategies(process_, scheme_, payoff, control_mean_, steps_, T_); + const size_t K = std::min(config.n_compete, all.size()); + if (K == 0) + return std::unexpected("No strategies available"); + + // ── Per-strategy running accumulators ───────────────────────────────── + // We track (sum, sum_sq, n) across all rounds so precision weights can + // be recomputed from the full history after each round. + struct Accumulator { + double sum = 0.0; + double sum_sq = 0.0; + size_t n = 0; + + void add(const std::vector& samples) { + for (double x : samples) { + sum += x; + sum_sq += x * x; + ++n; + } + } + + double mean() const { + return n > 0 ? sum / static_cast(n) : 0.0; + } + + // Unbiased sample variance (n-1 denominator) + double variance() const { + if (n < 2) return std::numeric_limits::infinity(); + const double m = mean(); + return (sum_sq - static_cast(n) * m * m) / + static_cast(n - 1); + } + }; + + std::vector accum(K); + + // ── Helper: compute precision weights from current accumulators ──────── + auto precision_weights = [&]() -> std::vector { + std::vector prec(K); + double total = 0.0; + for (size_t k = 0; k < K; ++k) { + double var = accum[k].variance(); + prec[k] = (var > 0.0 && std::isfinite(var)) ? 1.0 / var : 0.0; + total += prec[k]; + } + std::vector w(K); + if (total > 0.0) { + for (size_t k = 0; k < K; ++k) + w[k] = prec[k] / total; + } else { + std::fill(w.begin(), w.end(), 1.0 / static_cast(K)); + } + return w; + }; + + size_t leader = 0; // index of current leader + bool first_round = true; + std::vector history; + history.reserve(config.n_rounds); + + size_t total_samples = 0; + + for (size_t round = 0; round < config.n_rounds; ++round) { + // ── Determine thread allocation for this round ───────────────────── + // Every strategy gets 1 compete thread. + // The current leader also gets n_exploit bonus threads. + // On the first round there is no established leader yet, so the + // bonus threads are spread evenly (floor(n_exploit/K) each, with + // remainder going to strategy 0 as a tiebreak). + + std::vector thread_counts(K, 1); // 1 compete thread each + if (first_round) { + const size_t bonus_each = config.n_exploit / K; + const size_t leftover = config.n_exploit % K; + for (size_t k = 0; k < K; ++k) + thread_counts[k] += bonus_each; + thread_counts[0] += leftover; + } else { + thread_counts[leader] += config.n_exploit; + } + + // ── Run all threads for this round ───────────────────────────────── + // Each thread generates samples_per_thread samples. + // Seed = make_seed(strategy_idx * 10000 + thread_idx, round) ensures + // independence across strategies, threads within a strategy, and rounds. + + // Collect samples per strategy (multiple threads concatenated) + std::vector> round_samples(K); + { + // Flatten to (strategy_idx, thread_idx) pairs for easy launch + struct Job { size_t k; size_t t; }; + std::vector jobs; + for (size_t k = 0; k < K; ++k) + for (size_t t = 0; t < thread_counts[k]; ++t) + jobs.push_back({k, t}); + + std::vector> thread_results(jobs.size()); + { + std::vector threads; + threads.reserve(jobs.size()); + for (size_t j = 0; j < jobs.size(); ++j) { + threads.emplace_back([&, j]() { + const size_t k = jobs[j].k; + const size_t t = jobs[j].t; + const uint64_t seed = detail::make_seed(k * 10000 + t, round); + thread_results[j] = all[k].second(config.samples_per_thread, seed); + }); + } + } // jthreads join here + + // Merge thread results into per-strategy buckets + for (size_t j = 0; j < jobs.size(); ++j) { + auto& dest = round_samples[jobs[j].k]; + auto& src = thread_results[j]; + dest.insert(dest.end(), src.begin(), src.end()); + } + } + + // ── Update accumulators ──────────────────────────────────────────── + for (size_t k = 0; k < K; ++k) { + accum[k].add(round_samples[k]); + total_samples += round_samples[k].size(); + } + + // ── Recompute precision weights and current leader ───────────────── + auto weights = precision_weights(); + size_t new_leader = static_cast( + std::max_element(weights.begin(), weights.end()) - weights.begin()); + + // ── Build round diagnostic ───────────────────────────────────────── + double est = 0.0; + for (size_t k = 0; k < K; ++k) + est += weights[k] * accum[k].mean(); + + // Var(precision-weighted sum) = sum_k w_k^2 * sigma_k^2 / n_k + double var_est = 0.0; + for (size_t k = 0; k < K; ++k) { + double var_k = accum[k].variance(); + if (std::isfinite(var_k) && accum[k].n > 0) + var_est += weights[k] * weights[k] * var_k / + static_cast(accum[k].n); + } + + history.push_back({ + .round = round + 1, + .leader = all[new_leader].first, + .leader_changed = (!first_round && new_leader != leader), + .estimate = est, + .estimated_stderr = std::sqrt(var_est), + .precision_weights = weights, + }); + + leader = new_leader; + first_round = false; + } + + // ── Final estimate ──────────────────────────────────────────────────── + auto final_weights = precision_weights(); + double final_est = 0.0; + for (size_t k = 0; k < K; ++k) + final_est += final_weights[k] * accum[k].mean(); + + double final_var = 0.0; + for (size_t k = 0; k < K; ++k) { + double var_k = accum[k].variance(); + if (std::isfinite(var_k) && accum[k].n > 0) + final_var += final_weights[k] * final_weights[k] * var_k / + static_cast(accum[k].n); + } + + // ── Pack cumulative StrategyStats for diagnostics ────────────────────── + std::vector final_stats(K); + for (size_t k = 0; k < K; ++k) { + final_stats[k] = { + .name = all[k].first, + .mean = accum[k].mean(), + .sample_variance = accum[k].variance(), + .wall_time_ms = 0.0, // not tracked per-strategy across rounds + .n_samples = accum[k].n, + }; + } + + return IterativeEngineResult{ + .estimate = final_est, + .estimated_stderr = std::sqrt(final_var), + .final_leader = all[leader].first, + .round_history = history, + .final_stats = final_stats, + .total_samples = total_samples, + }; + } +}; + +} // namespace lfmc diff --git a/include/lfmc/payoff.hpp b/include/lfmc/payoff.hpp index 872ea17..1148b64 100644 --- a/include/lfmc/payoff.hpp +++ b/include/lfmc/payoff.hpp @@ -53,4 +53,62 @@ class ControlVariatePayoff : public Payoff { generate_payoffs(const std::vector& paths) const override; }; +class AsianCall : public Payoff { + private: + double strike_; + + public: + explicit AsianCall(double strike); + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + +class AsianPut : public Payoff { + private: + double strike_; + + public: + explicit AsianPut(double strike); + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + +class UpAndOutCall : public Payoff { + private: + double strike_; + double barrier_; + + public: + UpAndOutCall(double strike, double barrier); + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + +class DownAndInPut : public Payoff { + private: + double strike_; + double barrier_; + + public: + DownAndInPut(double strike, double barrier); + + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + +class LookbackCall : public Payoff { + public: + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + +class LookbackPut : public Payoff { + public: + std::expected, std::string> + generate_payoffs(const std::vector& paths) const override; +}; + } // namespace lfmc diff --git a/include/lfmc/strategies.hpp b/include/lfmc/strategies.hpp new file mode 100644 index 0000000..0549148 --- /dev/null +++ b/include/lfmc/strategies.hpp @@ -0,0 +1,420 @@ +#pragma once + +// All 6 new variance reduction strategy factories. +// Each factory returns a SamplerFn: (n_samples, seed) -> vector. +// The seed guarantees RNG independence between strategies and ASVR phases. + +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/types.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace lfmc::detail { + +// ─── Acklam's inverse normal CDF ───────────────────────────────────────────── +// Rational approximation accurate to ~5e-9 over (0, 1). +// Reference: Peter J. Acklam, "An algorithm for computing the inverse normal +// cumulative distribution function", 2010. +inline double normal_icdf(double p) noexcept { + static constexpr double a[] = {-3.969683028665376e+01, 2.209460984245205e+02, + -2.759285104469687e+02, 1.383577518672690e+02, + -3.066479806614716e+01, 2.506628277459239e+00}; + static constexpr double b[] = {-5.447609879822406e+01, 1.615858368580409e+02, + -1.556989798598866e+02, 6.680131188771972e+01, + -1.328068155288572e+01}; + static constexpr double c[] = {-7.784894002430293e-03, -3.223964580411365e-01, + -2.400758277161838e+00, -2.549732539343734e+00, + 4.374664141464968e+00, 2.938163982698783e+00}; + static constexpr double d[] = {7.784695709041462e-03, 3.224671290700398e-01, + 2.445134137142996e+00, 3.754408661907416e+00}; + + if (p < 0.02425) { + const double q = std::sqrt(-2.0 * std::log(p)); + return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / + ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0); + } + if (p <= 0.97575) { + const double q = p - 0.5; + const double r = q * q; + return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / + (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0); + } + const double q = std::sqrt(-2.0 * std::log(1.0 - p)); + return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / + ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0); +} + +// ─── Halton radical inverse ─────────────────────────────────────────────────── +// Returns the i-th term (1-indexed) of the Van der Corput sequence in `base`. +// Output in (0, 1). Caller must pass i >= 1. +inline double halton(size_t i, int base) noexcept { + double result = 0.0; + double f = 1.0 / static_cast(base); + size_t idx = i; + while (idx > 0) { + result += f * static_cast(idx % static_cast(base)); + idx /= static_cast(base); + f /= static_cast(base); + } + return result; +} + +// ─── First n primes (sieve of Eratosthenes) ─────────────────────────────────── +// Used by Halton QMC: dimension d gets the d-th prime as its base. +inline std::vector first_n_primes(size_t n) { + std::vector primes; + primes.reserve(n); + for (int candidate = 2; primes.size() < n; ++candidate) { + bool is_prime = true; + for (int p : primes) { + if (p * p > candidate) + break; + if (candidate % p == 0) { + is_prime = false; + break; + } + } + if (is_prime) + primes.push_back(candidate); + } + return primes; +} + +} // namespace lfmc::detail + +namespace lfmc { + +// ─── Strategy 5: Stratified sampling ───────────────────────────────────────── +// Stratifies only the first Brownian increment dimension (the one with highest +// correlation to the terminal price). Remaining dimensions are pseudo-random. +// The n strata are [(i/n, (i+1)/n)] for i=0..n-1; one uniform is drawn from +// each stratum and converted to N(0,1) via normal_icdf. +template NS> +SamplerFn make_stratified_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::uniform_real_distribution udist{0.0, 1.0}; + std::normal_distribution ndist{0.0, 1.0}; + + std::vector samples; + samples.reserve(n); + + for (size_t i = 0; i < n; ++i) { + // First dimension: one draw from stratum i + const double u = (static_cast(i) + udist(rng)) / static_cast(n); + const double z0 = detail::normal_icdf(u); + + // Remaining dimensions: pseudo-random + Normals normals(steps); + normals[0] = z0; + for (size_t j = 1; j < steps; ++j) + normals[j] = ndist(rng); + + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0]); + } + return samples; + }; +} + +// ─── Strategy 6: Halton QMC with random shift ───────────────────────────────── +// Generates a `steps`-dimensional Halton sequence. Each dimension is shifted +// by a random uniform (derived from the seed) before converting to N(0,1). +// The random shift randomises the sequence per ASVR phase while preserving +// the low-discrepancy property within each call. +// +// Known limitation: Halton sequences have inter-dimensional correlation for +// large dimension counts (steps > ~20). ASVR will automatically assign low +// weight to this strategy for path-dependent options where all steps matter. +template NS> +SamplerFn make_halton_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::uniform_real_distribution udist{0.0, 1.0}; + + // Random shift per dimension (Owen-style randomisation simplified to a shift) + const auto primes = detail::first_n_primes(steps); + std::vector shifts(steps); + for (double& s : shifts) + s = udist(rng); + + std::vector samples; + samples.reserve(n); + + for (size_t i = 0; i < n; ++i) { + Normals normals(steps); + for (size_t d = 0; d < steps; ++d) { + // Clamp away from the boundary to avoid normal_icdf singularity + double u = std::fmod(detail::halton(i + 1, primes[d]) + shifts[d], 1.0); + u = std::clamp(u, 1e-10, 1.0 - 1e-10); + normals[d] = detail::normal_icdf(u); + } + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0]); + } + return samples; + }; +} + +// ─── Strategy 7: Importance sampling (exponential tilting) ─────────────────── +// Tilts all normal draws using a TOTAL drift parameter `theta`, distributed +// evenly as theta/sqrt(steps) per time step. This parameterisation makes theta +// dimensionally consistent regardless of the number of steps: +// +// Z_t ~ N(theta/sqrt(steps), 1) under importance measure Q +// per-step shift: delta = theta / sqrt(steps) +// +// The Radon-Nikodym correction weight is: +// w = exp(-delta * sum_t(Z_t) + 0.5 * delta^2 * steps) +// = exp(-theta/sqrt(steps) * sum_t(Z_t) + 0.5 * theta^2) +// +// This is equivalent to shifting the terminal log-price by sigma*theta*sqrt(dt)*steps +// = sigma*theta*sqrt(T). It matches the single-normal exact simulation where shifting +// Z by theta yields the same log-price change. +// +// Choosing theta: +// theta = 0 → plain MC +// theta > 0 → tilts S_T upward; beneficial for OTM calls +// theta < 0 → tilts S_T downward; beneficial for OTM puts / barrier knock-ins +// practical range → |theta| in [0, 2]; beyond 2 the Radon-Nikodym weight +// exp(-0.5*theta^2) shrinks too fast and variance increases +template NS> +SamplerFn make_importance_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T, double theta) { + return [process, scheme, payoff, steps, T, theta](size_t n, + uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + // Per-step shift: distribute the total theta equally across sqrt(steps) normal draws + const double delta = theta / std::sqrt(static_cast(steps)); + std::normal_distribution ndist{delta, 1.0}; // tilted draw + + // Log-weight constant: 0.5 * delta^2 * steps = 0.5 * theta^2 + const double lw_const = 0.5 * theta * theta; + + std::vector samples; + samples.reserve(n); + + for (size_t i = 0; i < n; ++i) { + Normals normals(steps); + double z_sum = 0.0; + for (size_t j = 0; j < steps; ++j) { + normals[j] = ndist(rng); + z_sum += normals[j]; + } + // Likelihood ratio dP/dQ = exp(-delta * sum(Z') + 0.5 * delta^2 * steps) + // = exp(-theta/sqrt(steps) * sum(Z') + 0.5 * theta^2) + const double weight = std::exp(-delta * z_sum + lw_const); + + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0] * weight); + } + return samples; + }; +} + +// ─── Strategy 8: Moment matching ───────────────────────────────────────────── +// Generates all n * steps normals up front (step-major layout for cache +// efficiency in the rescaling pass), then rescales each time step's n draws +// to have exact sample mean = 0 and sample std = 1. This eliminates first- +// order sampling error in every dimension simultaneously. +// +// Memory: O(n * steps) doubles. For n=10,000 and steps=52 this is ~4 MB. +template NS> +SamplerFn make_moment_matching_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::normal_distribution ndist{0.0, 1.0}; + + // Step-major storage: z[j][i] = step j, path i + std::vector> z(steps, std::vector(n)); + for (size_t j = 0; j < steps; ++j) + for (size_t i = 0; i < n; ++i) + z[j][i] = ndist(rng); + + // Rescale each step independently to mean=0, std=1 + for (size_t j = 0; j < steps; ++j) { + const double nd = static_cast(n); + double mean = 0.0; + for (double v : z[j]) + mean += v; + mean /= nd; + + double var = 0.0; + for (double v : z[j]) { + double d = v - mean; + var += d * d; + } + var /= (nd - 1.0); + const double sd = std::sqrt(var); + + if (sd > 1e-12) { + for (double& v : z[j]) + v = (v - mean) / sd; + } + } + + // Simulate paths from the moment-matched normals + std::vector samples; + samples.reserve(n); + for (size_t i = 0; i < n; ++i) { + Normals normals(steps); + for (size_t j = 0; j < steps; ++j) + normals[j] = z[j][i]; + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0]); + } + return samples; + }; +} + +// ─── Strategy 9: Latin Hypercube Sampling ──────────────────────────────────── +// In each of the `steps` dimensions independently: generates n stratified +// uniform samples (one per stratum), randomly permutes them (so different +// dimensions are uncorrelated), then converts to N(0,1) via normal_icdf. +// Guarantees exact coverage of every marginal distribution. +template NS> +SamplerFn make_lhs_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::uniform_real_distribution udist{0.0, 1.0}; + + const double nd = static_cast(n); + + // Build the LHS matrix: z[j][i] = normal for path i, step j + // Column (step) j is a shuffled stratified sample + std::vector> z(steps, std::vector(n)); + for (size_t j = 0; j < steps; ++j) { + // Stratified uniforms + for (size_t i = 0; i < n; ++i) { + double u = (static_cast(i) + udist(rng)) / nd; + u = std::clamp(u, 1e-10, 1.0 - 1e-10); + z[j][i] = detail::normal_icdf(u); + } + // Fisher-Yates shuffle to decorrelate dimensions + for (size_t i = n - 1; i > 0; --i) { + std::uniform_int_distribution idist{0, i}; + std::swap(z[j][i], z[j][idist(rng)]); + } + } + + std::vector samples; + samples.reserve(n); + for (size_t i = 0; i < n; ++i) { + Normals normals(steps); + for (size_t j = 0; j < steps; ++j) + normals[j] = z[j][i]; + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0]); + } + return samples; + }; +} + +// ─── Strategy 10: Stratified + Antithetic ──────────────────────────────────── +// Combines stratified sampling on the first dimension with antithetic variates. +// For n requested samples, generates n pairs: +// - Stratum i: u_i = (i + U) / n → z_1 = normal_icdf(u_i) +// - Positive path: (z_1, z_2, ..., z_steps) with z_{2..steps} ~ N(0,1) +// - Negative path: (-z_1, z_2, ..., z_steps) with the SAME remaining draws +// - Returned sample i = 0.5 * (payoff(pos_path) + payoff(neg_path)) +// +// Sharing z_{2..steps} between both paths of a pair gives variance reduction +// not only from the first-dimension antithetic pairing but also from the +// correlated higher dimensions within the pair. +template NS> +SamplerFn make_stratified_antithetic_sampler(SP process, NS scheme, std::shared_ptr payoff, + size_t steps, double T) { + return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + std::uniform_real_distribution udist{0.0, 1.0}; + std::normal_distribution ndist{0.0, 1.0}; + + std::vector samples; + samples.reserve(n); + + for (size_t i = 0; i < n; ++i) { + // Stratified first dimension + const double u = (static_cast(i) + udist(rng)) / static_cast(n); + const double z0 = detail::normal_icdf(u); + + // Shared higher dimensions + Normals normals_pos(steps), normals_neg(steps); + normals_pos[0] = z0; + normals_neg[0] = -z0; + for (size_t j = 1; j < steps; ++j) { + const double zj = ndist(rng); + normals_pos[j] = zj; + normals_neg[j] = zj; // shared (not antithetic in remaining dims) + } + + auto path_pos = detail::generate_path(process, scheme, normals_pos, steps, T); + auto path_neg = detail::generate_path(process, scheme, normals_neg, steps, T); + + auto r_pos = payoff->generate_payoffs({path_pos}); + auto r_neg = payoff->generate_payoffs({path_neg}); + + if (r_pos && r_neg && !r_pos->empty() && !r_neg->empty() && + !(*r_pos)[0].empty() && !(*r_neg)[0].empty()) + samples.push_back(0.5 * ((*r_pos)[0][0] + (*r_neg)[0][0])); + } + return samples; + }; +} + +// ─── Convenience: build all 10 standard strategies ─────────────────────────── +// Assembles the canonical 10-strategy vector for a given payoff on GBM + +// Euler-Maruyama. The control variate strategies use EuropeanCall(0.0) as the +// control (= S_T) with control_mean = S0 * exp(mu * T). +// For path-dependent options where S_T is a poor control, the CV strategies +// will receive low precision weights from ASVR automatically. +template NS> +std::vector> +build_all_strategies(SP process, NS scheme, std::shared_ptr payoff, double control_mean, + size_t steps, double T, double is_theta = 0.5) { + // The control payoff is S_T: EuropeanCall(0) gives max(S_T - 0, 0) = S_T for GBM. + auto control_payoff = []() { return std::make_shared(0.0); }; + + return { + {"plain_mc", make_plain_mc_sampler(process, scheme, payoff, steps, T)}, + {"antithetic", make_antithetic_sampler(process, scheme, payoff, steps, T)}, + {"control_variate", + make_control_variate_sampler(process, scheme, payoff, control_payoff(), control_mean, + steps, T)}, + {"antithetic_cv", + make_antithetic_cv_sampler(process, scheme, payoff, control_payoff(), control_mean, steps, + T)}, + {"stratified", make_stratified_sampler(process, scheme, payoff, steps, T)}, + {"halton_qmc", make_halton_sampler(process, scheme, payoff, steps, T)}, + {"importance_sampling", + make_importance_sampler(process, scheme, payoff, steps, T, is_theta)}, + {"moment_matching", make_moment_matching_sampler(process, scheme, payoff, steps, T)}, + {"lhs", make_lhs_sampler(process, scheme, payoff, steps, T)}, + {"stratified_antithetic", + make_stratified_antithetic_sampler(process, scheme, payoff, steps, T)}, + }; +} + +} // namespace lfmc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5d3f5f9..059ac06 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,6 @@ # Source files for the lfmc library set(LFMC_SOURCES + adaptive/adaptive_estimator.cpp estimator/monte_carlo_estimator.cpp estimator/control_variate_estimator.cpp payoff/asian_payoffs.cpp diff --git a/src/adaptive/adaptive_estimator.cpp b/src/adaptive/adaptive_estimator.cpp new file mode 100644 index 0000000..8a1b86f --- /dev/null +++ b/src/adaptive/adaptive_estimator.cpp @@ -0,0 +1,233 @@ +#include "lfmc/adaptive_estimator.hpp" + +#include +#include +#include + +namespace lfmc { + +// ─── run_strategy_batch ─────────────────────────────────────────────────────── +// Calls sampler(n_samples, seed), measures wall time, computes mean and +// unbiased sample variance. + +StrategyStats AdaptiveVarianceReduction::run_strategy_batch(const std::string& name, + const SamplerFn& sampler, + size_t n_samples, uint64_t seed) { + auto t0 = std::chrono::high_resolution_clock::now(); + auto samples = sampler(n_samples, seed); + auto t1 = std::chrono::high_resolution_clock::now(); + + const double wall_ms = std::chrono::duration(t1 - t0).count(); + auto [mean, variance] = detail::mean_variance(samples); + + return StrategyStats{ + .name = name, + .mean = mean, + .sample_variance = variance, + .wall_time_ms = wall_ms, + .n_samples = samples.size(), + }; +} + +// ─── compute_precision_weights ──────────────────────────────────────────────── +// w_k* = (1 / sigma_k^2) / sum_j (1 / sigma_j^2) +// This is the inverse-variance (precision) weighting that minimises the variance +// of a linear combination of independent unbiased estimators. + +std::vector +AdaptiveVarianceReduction::compute_precision_weights(const std::vector& stats) { + std::vector precisions(stats.size()); + double total = 0.0; + for (size_t k = 0; k < stats.size(); ++k) { + precisions[k] = (stats[k].sample_variance > 0.0) ? 1.0 / stats[k].sample_variance : 0.0; + total += precisions[k]; + } + + std::vector weights(stats.size()); + if (total > 0.0) { + for (size_t k = 0; k < stats.size(); ++k) + weights[k] = precisions[k] / total; + } else { + // Degenerate: all variances are zero — use equal weights + std::fill(weights.begin(), weights.end(), 1.0 / static_cast(stats.size())); + } + return weights; +} + +// ─── run ────────────────────────────────────────────────────────────────────── + +ASVRResult AdaptiveVarianceReduction::run(std::vector> strategies, + size_t total_samples, ASVRConfig config) { + const size_t K = strategies.size(); + assert(K >= 1); + + const size_t hw = std::max(size_t{1}, config.n_threads); + + // ── Budget split ───────────────────────────────────────────────────────── + const size_t n_explore_each = std::max( + config.min_exploration_per_strategy, + static_cast(config.exploration_fraction * static_cast(total_samples) / + static_cast(K))); + const size_t n_explore_total = K * n_explore_each; + const size_t n_exploit = + (n_explore_total < total_samples) ? total_samples - n_explore_total : 0; + + // ── Phase 1: Parallel exploration ───────────────────────────────────────── + // Each strategy runs on its own thread with seed make_seed(k, 0). + // Thread count is capped at min(K, floor(alpha * hw)) so that exploration + // does not consume more than its "10%" thread share; remaining HW threads + // are left idle during this phase and fully used during exploitation. + const size_t explore_thread_cap = + std::max(size_t{1}, + std::min(K, static_cast(config.exploration_fraction * + static_cast(hw)))); + + std::vector exploration_stats(K); + + // Launch exploration in waves of explore_thread_cap threads + for (size_t wave_start = 0; wave_start < K; wave_start += explore_thread_cap) { + const size_t wave_end = std::min(wave_start + explore_thread_cap, K); + std::vector threads; + threads.reserve(wave_end - wave_start); + for (size_t k = wave_start; k < wave_end; ++k) { + threads.emplace_back([&, k]() { + exploration_stats[k] = run_strategy_batch(strategies[k].first, + strategies[k].second, n_explore_each, + detail::make_seed(k, 0)); + }); + } + // jthreads join on destruction at end of wave block + } + + // ── Compute precision weights ────────────────────────────────────────────── + auto weights = compute_precision_weights(exploration_stats); + + // ── Allocate exploitation budget ─────────────────────────────────────────── + // n_k = floor(w_k* * n_exploit); remainder given to the highest-weight strategy + std::vector exploit_counts(K, 0); + if (n_exploit > 0) { + size_t assigned = 0; + for (size_t k = 0; k < K; ++k) { + exploit_counts[k] = + static_cast(weights[k] * static_cast(n_exploit)); + assigned += exploit_counts[k]; + } + if (assigned < n_exploit) { + const size_t best_k = + static_cast(std::max_element(weights.begin(), weights.end()) - + weights.begin()); + exploit_counts[best_k] += (n_exploit - assigned); + } + } + + // ── Phase 2: Parallel exploitation ──────────────────────────────────────── + // Remaining (1 - alpha) * hw threads are available; strategies with non-zero + // allocation run in parallel. Seed make_seed(k, 1) is independent of phase 0. + const size_t exploit_thread_cap = std::max(size_t{1}, hw); + + std::vector exploit_stats(K); + { + std::vector active; + for (size_t k = 0; k < K; ++k) + if (exploit_counts[k] > 0) + active.push_back(k); + + for (size_t wave_start = 0; wave_start < active.size(); wave_start += exploit_thread_cap) { + const size_t wave_end = std::min(wave_start + exploit_thread_cap, active.size()); + std::vector threads; + threads.reserve(wave_end - wave_start); + for (size_t i = wave_start; i < wave_end; ++i) { + const size_t k = active[i]; + threads.emplace_back([&, k]() { + exploit_stats[k] = run_strategy_batch(strategies[k].first, + strategies[k].second, exploit_counts[k], + detail::make_seed(k, 1)); + }); + } + } + } + + // ── Combine ─────────────────────────────────────────────────────────────── + // + // The final estimate uses ONLY exploitation samples: + // + // mu_ASVR = sum_k w_k* * mu_k^exploit + // + // This is the precision-weighted fusion of K independent estimators where + // w_k* were chosen to minimise Var(mu_ASVR). + // + // Unbiasedness: each mu_k^exploit is unbiased for mu (the true price) because + // it was generated with seed make_seed(k,1), which is statistically independent + // of make_seed(k,0) used to compute w_k*. A weighted sum of unbiased estimators + // is unbiased regardless of how the (data-dependent) weights are chosen, provided + // the weights and the estimators they multiply are independent. + // + // The exploration samples are intentionally NOT included in the final estimate. + // Using them would require equal-weight pooling (because precision weights are + // correlated with the exploration means), and the equal-weight pooled exploration + // mean adds noise without guaranteeing variance improvement. The exploration + // samples are the "cost of learning" — they determine w_k* but do not contribute + // to the final estimate. This makes the overhead explicit and the estimator clean. + + double mu_exploit = 0.0; + double w_sum = 0.0; + for (size_t k = 0; k < K; ++k) { + if (exploit_counts[k] > 0) { + mu_exploit += weights[k] * exploit_stats[k].mean; + w_sum += weights[k]; + } + } + if (w_sum > 0.0 && w_sum < 1.0 - 1e-12) + mu_exploit /= w_sum; // renormalise if some strategies got zero allocation + + const double estimate = (n_exploit > 0) ? mu_exploit : exploration_stats[0].mean; + + // ── Variance estimation ──────────────────────────────────────────────────── + // + // Var(mu_ASVR) = Var(mu_exploit) = sum_k w_k*^2 * sigma_k^2 / n_k + // + // With precision-optimal weights and allocation this equals: + // 1 / sum_k (n_k / sigma_k^2) + // + // We use the general formula so it holds when n_k values are rounded integers + // and not all strategies receive non-zero allocation. + // + // sigma_k^2 is estimated from exploration sample variance (unbiased). + + double var_exploit = 0.0; + for (size_t k = 0; k < K; ++k) { + if (exploit_counts[k] > 0) { + var_exploit += weights[k] * weights[k] * exploration_stats[k].sample_variance / + static_cast(exploit_counts[k]); + } + } + // Renormalise in the same proportion as the estimate was renormalised + if (w_sum > 0.0 && w_sum < 1.0 - 1e-12) + var_exploit /= (w_sum * w_sum); + + const double total_variance = (n_exploit > 0) ? var_exploit + : exploration_stats[0].sample_variance / + static_cast(n_explore_each); + + // Plain MC baseline for comparison: sigma_0^2 / N using the first strategy + // (conventionally plain MC; the caller must ensure strategies[0] is plain MC) + const double plain_mc_var = + exploration_stats[0].sample_variance / static_cast(total_samples); + const double vr_ratio = (total_variance > 0.0) ? plain_mc_var / total_variance : 1.0; + + return ASVRResult{ + .estimate = estimate, + .estimated_variance = total_variance, + .estimated_stderr = std::sqrt(total_variance), + .exploration_stats = exploration_stats, + .precision_weights = weights, + .exploitation_counts = exploit_counts, + .plain_mc_variance_estimate = plain_mc_var, + .variance_reduction_ratio = vr_ratio, + .n_exploration = n_explore_total, + .n_exploitation = n_exploit, + .total_samples = total_samples, + }; +} + +} // namespace lfmc diff --git a/src/payoff/asian_payoffs.cpp b/src/payoff/asian_payoffs.cpp index f58707e..749a5e9 100644 --- a/src/payoff/asian_payoffs.cpp +++ b/src/payoff/asian_payoffs.cpp @@ -9,58 +9,42 @@ namespace lfmc { -class AsianCall : public Payoff { - private: - double strike_; +AsianCall::AsianCall(double strike) : strike_(strike) {} - public: - explicit AsianCall(double strike) : strike_(strike) {} +std::expected, std::string> +AsianCall::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + payoffs.reserve(paths.size()); - std::expected, std::string> - generate_payoffs(const std::vector& paths) const override { - Payoffs payoffs; + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in AsianCall"); - payoffs.reserve(paths.size()); - - for (const auto& path : paths) { - - if (path.empty()) - return std::unexpected("Empty path encountered in AsianCall"); - - double mean = - std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); - payoffs.push_back(std::max(mean - strike_, 0.0)); - } - - return std::vector{payoffs}; + double mean = + std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); + payoffs.push_back(std::max(mean - strike_, 0.0)); } -}; -class AsianPut : public Payoff { - private: - double strike_; + return std::vector{payoffs}; +} - public: - explicit AsianPut(double strike) : strike_(strike) {} +AsianPut::AsianPut(double strike) : strike_(strike) {} - std::expected, std::string> - generate_payoffs(const std::vector& paths) const override { +std::expected, std::string> +AsianPut::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + payoffs.reserve(paths.size()); - Payoffs payoffs; + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in AsianPut"); - payoffs.reserve(paths.size()); - - for (const auto& path : paths) { - if (path.empty()) - return std::unexpected("Empty path encountered in AsianPut"); - - double mean = - std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); - payoffs.push_back(std::max(strike_ - mean, 0.0)); - } - - return std::vector{payoffs}; + double mean = + std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); + payoffs.push_back(std::max(strike_ - mean, 0.0)); } -}; + + return std::vector{payoffs}; +} } // namespace lfmc diff --git a/src/payoff/barrier_payoffs.cpp b/src/payoff/barrier_payoffs.cpp index 25ff2ce..09833e4 100644 --- a/src/payoff/barrier_payoffs.cpp +++ b/src/payoff/barrier_payoffs.cpp @@ -8,68 +8,44 @@ namespace lfmc { -class UpAndOutCall : public Payoff { - private: - double strike_; - double barrier_; +UpAndOutCall::UpAndOutCall(double strike, double barrier) : strike_(strike), barrier_(barrier) {} - public: - UpAndOutCall(double strike, double barrier) : strike_(strike), barrier_(barrier) {} +std::expected, std::string> +UpAndOutCall::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + payoffs.reserve(paths.size()); - std::expected, std::string> - generate_payoffs(const std::vector& paths) const override { + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in UpAndOutCall"); - Payoffs payoffs; - payoffs.reserve(paths.size()); + bool knocked_out = + std::any_of(path.begin(), path.end(), [this](double s) { return s >= barrier_; }); - for (const auto& path : paths) { - if (path.empty()) - return std::unexpected("Empty path encountered in UpAndOutCall"); - - bool knocked_out = - std::any_of(path.begin(), path.end(), [this](double s) { return s >= barrier_; }); - - if (knocked_out) { - payoffs.push_back(0.0); - } else { - payoffs.push_back(std::max(path.back() - strike_, 0.0)); - } - } - - return std::vector{payoffs}; + payoffs.push_back(knocked_out ? 0.0 : std::max(path.back() - strike_, 0.0)); } -}; -class DownAndInPut : public Payoff { - private: - double strike_; - double barrier_; + return std::vector{payoffs}; +} - public: - DownAndInPut(double strike, double barrier) : strike_(strike), barrier_(barrier) {} +DownAndInPut::DownAndInPut(double strike, double barrier) : strike_(strike), barrier_(barrier) {} - std::expected, std::string> - generate_payoffs(const std::vector& paths) const override { - Payoffs payoffs; - payoffs.reserve(paths.size()); +std::expected, std::string> +DownAndInPut::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + payoffs.reserve(paths.size()); - for (const auto& path : paths) { + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in DownAndInPut"); - if (path.empty()) - return std::unexpected("Empty path encountered in DownAndInPut"); + bool knocked_in = + std::any_of(path.begin(), path.end(), [this](double s) { return s <= barrier_; }); - bool knocked_in = - std::any_of(path.begin(), path.end(), [this](double s) { return s <= barrier_; }); - - if (knocked_in) { - payoffs.push_back(std::max(strike_ - path.back(), 0.0)); + payoffs.push_back(knocked_in ? std::max(strike_ - path.back(), 0.0) : 0.0); + } - } else { - payoffs.push_back(0.0); - } - } + return std::vector{payoffs}; +} - return std::vector{payoffs}; - } -}; } // namespace lfmc diff --git a/src/payoff/lookback_payoffs.cpp b/src/payoff/lookback_payoffs.cpp index 40e98aa..25c21d8 100644 --- a/src/payoff/lookback_payoffs.cpp +++ b/src/payoff/lookback_payoffs.cpp @@ -8,44 +8,36 @@ namespace lfmc { -class LookbackCall : public Payoff { - public: - std::expected, std::string> +std::expected, std::string> +LookbackCall::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + payoffs.reserve(paths.size()); - generate_payoffs(const std::vector& paths) const override { - Payoffs payoffs; - payoffs.reserve(paths.size()); + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in LookbackCall"); - for (const auto& path : paths) { - if (path.empty()) - return std::unexpected("Empty path encountered in LookbackCall"); - - double min_price = *std::min_element(path.begin(), path.end()); - payoffs.push_back(path.back() - min_price); - } - - return std::vector{payoffs}; + double min_price = *std::min_element(path.begin(), path.end()); + payoffs.push_back(path.back() - min_price); } -}; -class LookbackPut : public Payoff { - public: - std::expected, std::string> - generate_payoffs(const std::vector& paths) const override { - Payoffs payoffs; - payoffs.reserve(paths.size()); + return std::vector{payoffs}; +} - for (const auto& path : paths) { - if (path.empty()) - return std::unexpected("Empty path encountered in LookbackPut"); +std::expected, std::string> +LookbackPut::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + payoffs.reserve(paths.size()); - double max_price = *std::max_element(path.begin(), path.end()); + for (const auto& path : paths) { + if (path.empty()) + return std::unexpected("Empty path encountered in LookbackPut"); - payoffs.push_back(max_price - path.back()); - } - - return std::vector{payoffs}; + double max_price = *std::max_element(path.begin(), path.end()); + payoffs.push_back(max_price - path.back()); } -}; + + return std::vector{payoffs}; +} } // namespace lfmc diff --git a/src/payoffs/european_payoffs.cpp b/src/payoffs/european_payoffs.cpp new file mode 100644 index 0000000..29b5fef --- /dev/null +++ b/src/payoffs/european_payoffs.cpp @@ -0,0 +1,33 @@ +#include "lfmc/payoff.hpp" + +#include +#include +#include + +namespace lfmc { + +EuropeanCall::EuropeanCall(double strike) : strike_(strike) {} + +std::expected, std::string> +EuropeanCall::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + for (const auto& path : paths) { + double final_price = path.back(); + payoffs.push_back(std::max(final_price - strike_, 0.0)); + } + return std::vector{payoffs}; +} + +EuropeanPut::EuropeanPut(double strike) : strike_(strike) {} + +std::expected, std::string> +EuropeanPut::generate_payoffs(const std::vector& paths) const { + Payoffs payoffs; + for (const auto& path : paths) { + double final_price = path.back(); + payoffs.push_back(std::max(strike_ - final_price, 0.0)); + } + return std::vector{payoffs}; +} + +} // namespace lfmc diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7ae9a63..0bf1a99 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,12 +21,43 @@ set(TEST_SOURCES test_process.cpp test_scheme.cpp test_pipeline.cpp + test_adaptive.cpp + test_strategies.cpp + test_engine.cpp # test_convergency.cpp ) add_executable(tests ${TEST_SOURCES}) target_link_libraries(tests PRIVATE lfmc Catch2::Catch2WithMain) +# Standalone benchmark executable (not a Catch2 test — has its own main) +add_executable(bench bench_asvr_vs_fixed.cpp) +target_link_libraries(bench PRIVATE lfmc) + +# Quick comparison: Engine vs fixed strategies (same budget) +add_executable(quick_compare quick_compare.cpp) +target_link_libraries(quick_compare PRIVATE lfmc) +add_test(NAME quick_compare COMMAND quick_compare) + +# Copy UCRT64 C++ runtime DLLs so the executables find the right runtime on Windows. +# Without this, Git's bundled MinGW64 DLLs (on PATH earlier) take precedence and +# cause exit code 0xc0000139 (STATUS_ENTRYPOINT_NOT_FOUND). +if(MSYS OR MINGW) + set(UCRT64_BIN_DIR "C:/msys64/ucrt64/bin") + foreach(target tests bench) + foreach(dll libgcc_s_seh-1.dll libstdc++-6.dll libwinpthread-1.dll) + if(EXISTS "${UCRT64_BIN_DIR}/${dll}") + add_custom_command(TARGET ${target} POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy + "${UCRT64_BIN_DIR}/${dll}" + "$/${dll}" + VERBATIM + ) + endif() + endforeach() + endforeach() +endif() + catch_discover_tests( tests ADD_TAGS_AS_LABELS diff --git a/tests/bench_asvr_vs_fixed.cpp b/tests/bench_asvr_vs_fixed.cpp new file mode 100644 index 0000000..d061a7a --- /dev/null +++ b/tests/bench_asvr_vs_fixed.cpp @@ -0,0 +1,158 @@ +// bench_asvr_vs_fixed.cpp +// +// Critical publishability benchmark: +// Compare ASVR against EVERY fixed strategy (not just plain MC) +// using the same total sample budget. +// +// If ASVR beats the *best fixed strategy known in advance*, that is a strong result. +// If ASVR only beats plain MC but not antithetic+CV, that is a weak result. +// +// Run standalone: compile and link with lfmc, then ./bench + +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/strategies.hpp" + +#include +#include +#include +#include +#include +#include + +using namespace lfmc; + +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; +static const double CONTROL_MEAN = S0 * std::exp(MU * T); + +static GeometricBrownianMotion GBM{MU, SIGMA, S0}; +static EulerMaruyama EULER; + +// Run RUNS trials of each fixed strategy and ASVR; report average squared error +struct BenchResult { + std::string name; + double mean_sq_error; + double mean_estimate; + double std_error_of_mean; +}; + +BenchResult bench_fixed(const std::string& name, SamplerFn sampler, + double ref_mean, size_t N, size_t RUNS) { + std::vector errors; + errors.reserve(RUNS); + double sum_est = 0.0; + for (size_t r = 0; r < RUNS; ++r) { + auto samples = sampler(N, detail::make_seed(r + 9999, 44)); + double est = std::accumulate(samples.begin(), samples.end(), 0.0) / + static_cast(samples.size()); + double e = est - ref_mean; + errors.push_back(e * e); + sum_est += est; + } + double mse = std::accumulate(errors.begin(), errors.end(), 0.0) / static_cast(RUNS); + // std error of the MSE estimate (rough) + double mean_err2 = mse; + double sq_sum = 0.0; + for (double e2 : errors) { double d = e2 - mean_err2; sq_sum += d*d; } + double sem = std::sqrt(sq_sum / static_cast(RUNS * (RUNS-1))); + return {name, mse, sum_est / RUNS, sem}; +} + +BenchResult bench_asvr(std::shared_ptr payoff, double ref_mean, + size_t N, size_t RUNS) { + std::vector errors; + errors.reserve(RUNS); + double sum_est = 0.0; + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + for (size_t r = 0; r < RUNS; ++r) { + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + auto res = AdaptiveVarianceReduction::run(std::move(strategies), N, cfg); + double e = res.estimate - ref_mean; + errors.push_back(e * e); + sum_est += res.estimate; + } + double mse = std::accumulate(errors.begin(), errors.end(), 0.0) / static_cast(RUNS); + double mean_err2 = mse; + double sq_sum = 0.0; + for (double e2 : errors) { double d = e2 - mean_err2; sq_sum += d*d; } + double sem = std::sqrt(sq_sum / static_cast(RUNS * (RUNS-1))); + return {"ASVR", mse, sum_est / RUNS, sem}; +} + +void run_option_bench(const std::string& opt_name, std::shared_ptr payoff, + size_t N, size_t RUNS) { + printf("\n=== %-35s (N=%zu, runs=%zu) ===\n", opt_name.c_str(), N, RUNS); + + // Reference from 500k plain MC samples + auto plain_ref_fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto ref_samples = plain_ref_fn(500'000, detail::make_seed(0, 99)); + double ref_mean = std::accumulate(ref_samples.begin(), ref_samples.end(), 0.0) + / static_cast(ref_samples.size()); + printf(" Reference mean: %.6f\n", ref_mean); + + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + std::vector results; + + // Benchmark each fixed strategy + for (auto& [name, sampler] : strategies) { + results.push_back(bench_fixed(name, sampler, ref_mean, N, RUNS)); + } + + // Benchmark ASVR + results.push_back(bench_asvr(payoff, ref_mean, N, RUNS)); + + // Sort by MSE ascending + std::sort(results.begin(), results.end(), + [](const BenchResult& a, const BenchResult& b) { + return a.mean_sq_error < b.mean_sq_error; + }); + + double plain_mse = 0.0; + for (auto& r : results) if (r.name == "plain_mc") plain_mse = r.mean_sq_error; + + printf(" %-28s %12s %8s %8s\n", "Strategy", "MSE", "vs PlainMC", "vs ASVR"); + printf(" %-28s %12s %8s %8s\n", + std::string(28,'-').c_str(), std::string(12,'-').c_str(), + std::string(8,'-').c_str(), std::string(8,'-').c_str()); + + double asvr_mse = 0.0; + for (auto& r : results) if (r.name == "ASVR") asvr_mse = r.mean_sq_error; + + for (auto& r : results) { + double vs_plain = plain_mse / r.mean_sq_error; + double vs_asvr = asvr_mse / r.mean_sq_error; + printf(" %-28s %12.6f %8.3fx %8.3fx\n", + r.name.c_str(), r.mean_sq_error, vs_plain, vs_asvr); + } +} + +int main() { + const size_t N = 10'000; // sample budget per trial + const size_t RUNS = 200; // trials for MSE estimate + + printf("ASVR vs Fixed Strategy Benchmark\n"); + printf("=================================\n"); + printf("Budget N=%zu per trial, %zu independent runs\n", N, RUNS); + printf("Strategies compared: all 10 fixed + ASVR (adaptive)\n"); + + run_option_bench("European Call", std::make_shared(K), N, RUNS); + run_option_bench("European Put", std::make_shared(K), N, RUNS); + run_option_bench("Asian Call", std::make_shared(K), N, RUNS); + run_option_bench("Asian Put", std::make_shared(K), N, RUNS); + run_option_bench("Up-And-Out Call B=130", std::make_shared(K, 130.0), N, RUNS); + run_option_bench("Down-And-In Put B=70", std::make_shared(K, 70.0), N, RUNS); + run_option_bench("Lookback Call", std::make_shared(), N, RUNS); + run_option_bench("Lookback Put", std::make_shared(), N, RUNS); + + printf("\nDone.\n"); + return 0; +} diff --git a/tests/quick_compare.cpp b/tests/quick_compare.cpp new file mode 100644 index 0000000..b300912 --- /dev/null +++ b/tests/quick_compare.cpp @@ -0,0 +1,114 @@ +// Quick comparison: Engine (winner-takes-all) vs single fixed strategies +// Uses same total sample budget across all approaches. + +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/engine.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/strategies.hpp" + +#include +#include +#include +#include +#include + +using namespace lfmc; + +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; +static const double CONTROL_MEAN = S0 * std::exp(MU * T); + +static GeometricBrownianMotion GBM{MU, SIGMA, S0}; +static EulerMaruyama EULER; + +// Total budget: Engine uses 4*2000 (explore) + 8*5000 (exploit) = 48000 samples. +// Fixed strategies get the same 48000 samples for a fair comparison. +static constexpr size_t TOTAL = 48000; +static constexpr size_t RUNS = 50; + +static double ref_mean(std::shared_ptr payoff) { + auto fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto s = fn(500'000, detail::make_seed(0, 77)); + double sum = 0.0; + for (double x : s) sum += x; + return sum / static_cast(s.size()); +} + +static double bench_fixed(const std::string& name, SamplerFn fn, double ref, size_t runs) { + double mse = 0.0; + for (size_t r = 0; r < runs; ++r) { + auto s = fn(TOTAL, detail::make_seed(r + 1000, 5)); + double est = std::accumulate(s.begin(), s.end(), 0.0) / static_cast(s.size()); + double e = est - ref; + mse += e * e; + } + return mse / static_cast(runs); +} + +static double bench_engine(std::shared_ptr payoff, double ref, size_t runs) { + Engine> engine{ + GBM, EULER, CONTROL_MEAN, STEPS, T}; + + EngineConfig cfg; + cfg.n_explore_threads = 4; + cfg.explore_samples = 2000; + cfg.n_exploit_threads = 8; + cfg.exploit_samples = 5000; + + double mse = 0.0; + for (size_t r = 0; r < runs; ++r) { + auto res = engine.run(payoff, cfg); + if (!res) continue; + double e = res->estimate - ref; + mse += e * e; + } + return mse / static_cast(runs); +} + +void compare(const char* opt_name, std::shared_ptr payoff) { + printf("\n=== %s ===\n", opt_name); + const double ref = ref_mean(payoff); + printf("Reference: %.4f (budget per run: %zu samples)\n\n", ref, TOTAL); + + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + + printf(" %-28s %10s %8s\n", "Method", "MSE", "vs PlainMC"); + printf(" %-28s %10s %8s\n", "----------------------------", "----------", "--------"); + + double plain_mse = 0.0; + std::vector> rows; + + for (auto& [name, fn] : strategies) { + double mse = bench_fixed(name, fn, ref, RUNS); + rows.push_back({name, mse}); + if (name == "plain_mc") plain_mse = mse; + } + + double engine_mse = bench_engine(payoff, ref, RUNS); + rows.push_back({"[Engine]", engine_mse}); + + std::sort(rows.begin(), rows.end(), [](auto& a, auto& b) { return a.second < b.second; }); + + for (auto& [name, mse] : rows) { + double ratio = plain_mse / mse; + printf(" %-28s %10.6f %8.2fx\n", name.c_str(), mse, ratio); + } +} + +int main() { + printf("Engine vs Fixed Strategy Comparison\n"); + printf("Same total sample budget (%zu), %zu independent runs\n\n", TOTAL, RUNS); + + compare("European Call", std::make_shared(K)); + compare("Asian Call", std::make_shared(K)); + compare("Barrier (Up-Out Call B=130)", std::make_shared(K, 130.0)); + compare("Lookback Call", std::make_shared()); + + printf("\nDone.\n"); +} diff --git a/tests/test_adaptive.cpp b/tests/test_adaptive.cpp new file mode 100644 index 0000000..edd980e --- /dev/null +++ b/tests/test_adaptive.cpp @@ -0,0 +1,231 @@ +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" + +#include +#include +#include +#include + +using namespace lfmc; + +// ─── Black-Scholes closed-form price ───────────────────────────────────────── +// Used as the ground truth to verify unbiasedness. + +static double standard_normal_cdf(double x) { + return 0.5 * std::erfc(-x / std::sqrt(2.0)); +} + +static double bs_call(double S, double K, double r, double sigma, double T) { + const double d1 = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T) / (sigma * std::sqrt(T)); + const double d2 = d1 - sigma * std::sqrt(T); + return S * standard_normal_cdf(d1) - K * std::exp(-r * T) * standard_normal_cdf(d2); +} + +// ─── Shared test parameters ─────────────────────────────────────────────────── +// ATM European call on GBM: S0=100, K=100, mu=0.05, sigma=0.2, T=1 +// +// The library computes the UNDISCOUNTED expected payoff E[max(S_T-K,0)]. +// With mu = risk-free rate r, the physical and risk-neutral measures coincide, +// so the correct reference is: +// +// E[max(S_T-K,0)] = BS_call(S0,K,r,sigma,T) * exp(r*T) +// ≈ 10.45 * exp(0.05) ≈ 10.99 +// +// The Black-Scholes discounted price is NOT the right comparison target here. + +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; // weekly steps +// control: E[S_T] = S0 * exp(mu*T) +static const double CONTROL_MEAN = S0 * std::exp(MU * T); + +// Build the standard four-strategy suite +static std::vector> build_strategies() { + GeometricBrownianMotion gbm{MU, SIGMA, S0}; + EulerMaruyama euler; + + return { + {"plain_mc", + make_plain_mc_sampler(gbm, euler, std::make_shared(K), STEPS, T)}, + + {"antithetic", + make_antithetic_sampler(gbm, euler, std::make_shared(K), STEPS, T)}, + + // Control variate: target = call payoff, control = S_T (EuropeanCall(0) = max(S_T,0) = S_T + // for GBM), E[S_T] = S0 * exp(mu*T) + {"control_variate", + make_control_variate_sampler(gbm, euler, std::make_shared(K), + std::make_shared(0.0), CONTROL_MEAN, STEPS, + T)}, + + {"antithetic_cv", + make_antithetic_cv_sampler(gbm, euler, std::make_shared(K), + std::make_shared(0.0), CONTROL_MEAN, STEPS, T)}, + }; +} + +// ─── Tests ──────────────────────────────────────────────────────────────────── + +TEST_CASE("ASVR: estimate is close to undiscounted expected payoff") { + // ASVR with 20 000 total samples. + // True target: E[max(S_T-K,0)] = BS_call * exp(r*T) ≈ 10.99 + // We use a wide tolerance (±1.0) since the estimator has no discounting + // and discretization error from Euler-Maruyama adds a small systematic bias. + + const double undiscounted_target = bs_call(S0, K, MU, SIGMA, T) * std::exp(MU * T); + + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + + auto result = AdaptiveVarianceReduction::run(build_strategies(), 20'000, cfg); + + INFO("Undiscounted target: " << undiscounted_target); + INFO("ASVR estimate: " << result.estimate << " ± " << result.estimated_stderr); + REQUIRE(result.estimate > 0.0); + // Within ±1.0 of the true undiscounted payoff (very conservative bound) + REQUIRE(std::abs(result.estimate - undiscounted_target) < 1.0); +} + +TEST_CASE("ASVR: variance reduction ratio > 1 (beats plain MC)") { + // The core claim: same total sample budget, lower estimated variance. + // For ATM European calls on GBM, antithetic and CV each achieve >50% variance + // reduction. Even accounting for the 10% exploration overhead, ASVR should + // easily achieve VR ratio > 1. + + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + + auto result = AdaptiveVarianceReduction::run(build_strategies(), 20'000, cfg); + + INFO("VR ratio: " << result.variance_reduction_ratio); + INFO("ASVR variance: " << result.estimated_variance); + INFO("Plain MC variance: " << result.plain_mc_variance_estimate); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR: precision weights sum to 1") { + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + + auto result = AdaptiveVarianceReduction::run(build_strategies(), 10'000, cfg); + + double weight_sum = 0.0; + for (double w : result.precision_weights) + weight_sum += w; + + REQUIRE(weight_sum == Catch::Approx(1.0).epsilon(1e-10)); +} + +TEST_CASE("ASVR: exploration identifies a VR strategy as best") { + // The best strategy (lowest exploration variance) should NOT be plain_mc. + // For a European call, antithetic/CV should always outperform plain MC. + + ASVRConfig cfg; + cfg.exploration_fraction = 0.15; + cfg.n_threads = 4; + + auto result = AdaptiveVarianceReduction::run(build_strategies(), 20'000, cfg); + + // Find which strategy got the highest weight + const auto& weights = result.precision_weights; + const size_t best_k = static_cast( + std::max_element(weights.begin(), weights.end()) - weights.begin()); + + INFO("Best strategy: " << result.exploration_stats[best_k].name); + INFO("Weights: "); + for (size_t k = 0; k < weights.size(); ++k) + INFO(" " << result.exploration_stats[k].name << ": " << weights[k]); + + // plain_mc is index 0; it should not have the highest weight + REQUIRE(best_k != 0); +} + +TEST_CASE("ASVR: sample budget splits correctly") { + const size_t N = 10'000; + const size_t K = 4; + + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + cfg.min_exploration_per_strategy = 0; // allow the fraction to fully control the split + + auto result = AdaptiveVarianceReduction::run(build_strategies(), N, cfg); + + // Exploration: K * floor(0.1 * N / K) = K * floor(250) = 1000 + // Exploitation: 9000 + REQUIRE(result.n_exploration + result.n_exploitation == N); + REQUIRE(result.total_samples == N); + + // Exploitation counts sum to n_exploitation + size_t exploit_sum = 0; + for (size_t c : result.exploitation_counts) + exploit_sum += c; + REQUIRE(exploit_sum == result.n_exploitation); +} + +// ─── Empirical variance comparison ──────────────────────────────────────────── +// +// This test is the paper-quality benchmark: run both ASVR and plain MC many +// times, compute empirical MSE for each, and verify ASVR MSE < plain MC MSE. +// It is slower (~seconds) but produces the core empirical result. + +TEST_CASE("ASVR: empirical MSE lower than plain MC with same sample budget", "[.slow]") { + // Use a fixed reference: the mean of a large plain MC run as the "true" value. + // This avoids the discounting ambiguity and tests variance reduction in isolation. + // We compute a reference from 200k plain MC samples (essentially noiseless). + GeometricBrownianMotion gbm{MU, SIGMA, S0}; + EulerMaruyama euler; + SamplerFn plain = + make_plain_mc_sampler(gbm, euler, std::make_shared(K), STEPS, T); + + // Reference: 200k-sample plain MC mean (variance ≈ 0 relative to RUNS-run experiment) + auto ref_samples = plain(200'000, detail::make_seed(0, 42)); + double ref_mean = 0.0; + for (double x : ref_samples) + ref_mean += x; + ref_mean /= static_cast(ref_samples.size()); + + const size_t N = 10'000; // samples per run + const size_t RUNS = 100; // independent replications + + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + + double mse_asvr = 0.0; + double mse_plain = 0.0; + + for (size_t r = 0; r < RUNS; ++r) { + // ASVR run — N total samples + auto asvr_res = AdaptiveVarianceReduction::run(build_strategies(), N, cfg); + double err_asvr = asvr_res.estimate - ref_mean; + mse_asvr += err_asvr * err_asvr; + + // Plain MC run — same N samples + auto plain_samples = plain(N, detail::make_seed(r + 10000, 99)); + double plain_mean = 0.0; + for (double x : plain_samples) + plain_mean += x; + plain_mean /= static_cast(plain_samples.size()); + double err_plain = plain_mean - ref_mean; + mse_plain += err_plain * err_plain; + } + + mse_asvr /= static_cast(RUNS); + mse_plain /= static_cast(RUNS); + + INFO("Reference value: " << ref_mean); + INFO("Empirical MSE ASVR: " << mse_asvr); + INFO("Empirical MSE plain MC: " << mse_plain); + INFO("Empirical VR ratio: " << mse_plain / mse_asvr); + + REQUIRE(mse_asvr < mse_plain); +} diff --git a/tests/test_engine.cpp b/tests/test_engine.cpp new file mode 100644 index 0000000..39187a8 --- /dev/null +++ b/tests/test_engine.cpp @@ -0,0 +1,231 @@ +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/engine.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/strategies.hpp" + +#include +#include +#include +#include + +using namespace lfmc; + +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; +static const double CONTROL_MEAN = S0 * std::exp(MU * T); + +static Engine> make_engine() { + return {GeometricBrownianMotion{MU, SIGMA, S0}, + EulerMaruyama{}, CONTROL_MEAN, STEPS, T}; +} + +// ─── Smoke tests ────────────────────────────────────────────────────────────── + +TEST_CASE("Engine: returns a valid positive estimate for European call") { + auto engine = make_engine(); + auto result = engine.run(std::make_shared(K)); + + REQUIRE(result.has_value()); + REQUIRE(result->estimate > 0.0); + REQUIRE(result->estimated_stderr > 0.0); +} + +TEST_CASE("Engine: estimate is within reasonable range of Black-Scholes") { + auto engine = make_engine(); + auto result = engine.run(std::make_shared(K)); + + REQUIRE(result.has_value()); + WARN("Estimate: " << result->estimate << " ± " << result->estimated_stderr); + WARN("Leader: " << result->final_leader); + REQUIRE(std::abs(result->estimate - 10.99) < 2.0); +} + +TEST_CASE("Engine: round history has n_rounds entries") { + auto engine = make_engine(); + + IterativeEngineConfig cfg; + cfg.n_rounds = 3; + + auto result = engine.run(std::make_shared(K), cfg); + + REQUIRE(result.has_value()); + REQUIRE(result->round_history.size() == 3); + for (size_t i = 0; i < 3; ++i) + REQUIRE(result->round_history[i].round == i + 1); +} + +TEST_CASE("Engine: precision weights sum to 1 after each round") { + auto engine = make_engine(); + + IterativeEngineConfig cfg; + cfg.n_rounds = 4; + + auto result = engine.run(std::make_shared(K), cfg); + REQUIRE(result.has_value()); + + for (const auto& r : result->round_history) { + double sum = 0.0; + for (double w : r.precision_weights) sum += w; + REQUIRE(sum == Catch::Approx(1.0).epsilon(1e-10)); + } +} + +TEST_CASE("Engine: total samples consistent with config") { + auto engine = make_engine(); + + IterativeEngineConfig cfg; + cfg.n_compete = 4; + cfg.n_exploit = 8; + cfg.n_rounds = 3; + cfg.samples_per_thread = 500; + + auto result = engine.run(std::make_shared(K), cfg); + REQUIRE(result.has_value()); + + // Round 1: (4 compete + 8 exploit spread equally) = 4*(1+2)=12 threads * 500 + // + leftover 8%4=0, so 3 per strategy = 12 threads * 500 = 6000 + // Rounds 2-3: 4 compete (1 each) + 8 exploit on leader = 12 threads * 500 = 6000 each + // Total = 3 * 12 * 500 = 18000 + const size_t expected = cfg.n_rounds * + (cfg.n_compete + cfg.n_exploit) * + cfg.samples_per_thread; + REQUIRE(result->total_samples == expected); +} + +TEST_CASE("Engine: final stats cover all competing strategies") { + auto engine = make_engine(); + auto result = engine.run(std::make_shared(K)); + + REQUIRE(result.has_value()); + REQUIRE(result->final_stats.size() == 4); + for (const auto& s : result->final_stats) { + REQUIRE(s.n_samples > 0); + REQUIRE(s.sample_variance >= 0.0); + } +} + +TEST_CASE("Engine: works with Asian call payoff") { + auto engine = make_engine(); + auto result = engine.run(std::make_shared(K)); + + REQUIRE(result.has_value()); + REQUIRE(result->estimate > 0.0); +} + +TEST_CASE("Engine: works with barrier option payoff") { + auto engine = make_engine(); + auto result = engine.run(std::make_shared(K, 130.0)); + + REQUIRE(result.has_value()); + REQUIRE(result->estimate >= 0.0); +} + +// ─── Leader stability test ──────────────────────────────────────────────────── +// With enough samples per round, the leader should stabilise and not switch in +// late rounds (once we have enough data to distinguish strategies reliably). + +TEST_CASE("Engine: leader stabilises over rounds") { + auto engine = make_engine(); + + IterativeEngineConfig cfg; + cfg.n_compete = 4; + cfg.n_exploit = 8; + cfg.n_rounds = 6; + cfg.samples_per_thread = 1000; + + auto result = engine.run(std::make_shared(K), cfg); + REQUIRE(result.has_value()); + + // Count leader changes in the last half of rounds + size_t late_changes = 0; + const auto& h = result->round_history; + for (size_t i = h.size() / 2; i < h.size(); ++i) + if (h[i].leader_changed) ++late_changes; + + WARN("Leader per round:"); + for (const auto& r : h) + WARN(" Round " << r.round << ": " << r.leader + << (r.leader_changed ? " (SWITCHED)" : "")); + + // Allow at most 1 late switch — the leader might flip once as estimates refine + REQUIRE(late_changes <= 1); +} + +// ─── Empirical MSE comparison ───────────────────────────────────────────────── + +TEST_CASE("Engine: empirical MSE vs fixed strategies", "[bench]") { + static constexpr size_t TOTAL = 48000; + static constexpr size_t RUNS = 15; + + GeometricBrownianMotion gbm{MU, SIGMA, S0}; + EulerMaruyama euler; + + struct Row { std::string name; double mse; }; + + auto run_comparison = [&](const char* opt_name, std::shared_ptr payoff) { + auto plain_ref = make_plain_mc_sampler(gbm, euler, payoff, STEPS, T); + auto ref_s = plain_ref(300'000, detail::make_seed(0, 77)); + double ref = std::accumulate(ref_s.begin(), ref_s.end(), 0.0) / + static_cast(ref_s.size()); + + auto strategies = build_all_strategies(gbm, euler, payoff, CONTROL_MEAN, STEPS, T); + + std::vector rows; + + for (auto& [name, fn] : strategies) { + double mse = 0.0; + for (size_t r = 0; r < RUNS; ++r) { + auto s = fn(TOTAL, detail::make_seed(r + 1000, 5)); + double est = std::accumulate(s.begin(), s.end(), 0.0) / + static_cast(s.size()); + double e = est - ref; + mse += e * e; + } + rows.push_back({name, mse / RUNS}); + } + + // Engine: 4 compete + 8 exploit, 5 rounds, 800 samples/thread + // = 5 * 12 * 800 = 48000 total + Engine> engine{ + gbm, euler, CONTROL_MEAN, STEPS, T}; + IterativeEngineConfig cfg; + cfg.n_compete = 4; + cfg.n_exploit = 8; + cfg.n_rounds = 5; + cfg.samples_per_thread = 800; + { + double mse = 0.0; + for (size_t r = 0; r < RUNS; ++r) { + auto res = engine.run(payoff, cfg); + REQUIRE(res.has_value()); + double e = res->estimate - ref; + mse += e * e; + } + rows.push_back({"[Engine]", mse / RUNS}); + } + + std::sort(rows.begin(), rows.end(), [](auto& a, auto& b) { return a.mse < b.mse; }); + double plain_mse = 0.0; + for (auto& r : rows) if (r.name == "plain_mc") plain_mse = r.mse; + + WARN(opt_name << " (ref=" << ref << ", budget=" << TOTAL << ", runs=" << RUNS << ")"); + for (auto& r : rows) + WARN(" " << r.name << ": MSE=" << r.mse + << " (" << (plain_mse / r.mse) << "x vs plain_mc)"); + + double engine_mse = 0.0; + for (auto& r : rows) if (r.name == "[Engine]") engine_mse = r.mse; + REQUIRE(engine_mse < plain_mse); + }; + + run_comparison("European Call", std::make_shared(K)); + run_comparison("Asian Call", std::make_shared(K)); + run_comparison("Barrier (Up-Out B=130)", std::make_shared(K, 130.0)); + run_comparison("Lookback Call", std::make_shared()); +} diff --git a/tests/test_strategies.cpp b/tests/test_strategies.cpp new file mode 100644 index 0000000..9691602 --- /dev/null +++ b/tests/test_strategies.cpp @@ -0,0 +1,356 @@ +#include "lfmc/adaptive_estimator.hpp" +#include "lfmc/numerical_scheme.hpp" +#include "lfmc/payoff.hpp" +#include "lfmc/stochastic_process.hpp" +#include "lfmc/strategies.hpp" + +#include +#include +#include +#include +#include + +using namespace lfmc; + +// ─── Shared market parameters ───────────────────────────────────────────────── +// ATM European-style: S0=100, K=100, mu=r=0.05, sigma=0.2, T=1y, 52 steps +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double BARRIER = 130.0; // up-and-out / down-and-in barrier +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; + +static const double CONTROL_MEAN = S0 * std::exp(MU * T); // E[S_T] = 105.13 + +static GeometricBrownianMotion GBM{MU, SIGMA, S0}; +static EulerMaruyama EULER; + +// Compute mean and sample variance of a vector +static std::pair stats(const std::vector& v) { + const double n = static_cast(v.size()); + double sum = 0.0; + for (double x : v) sum += x; + const double mean = sum / n; + double sq = 0.0; + for (double x : v) { double d = x - mean; sq += d * d; } + return {mean, sq / (n - 1.0)}; +} + +// ─── Smoke tests: each sampler returns the requested count ──────────────────── + +TEST_CASE("strategies: all 10 samplers return requested sample count") { + auto payoff = std::make_shared(K); + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + + for (const auto& [name, sampler] : strategies) { + INFO("Strategy: " << name); + auto samples = sampler(100, detail::make_seed(0, 0)); + REQUIRE(samples.size() == 100); + } +} + +// ─── Unbiasedness: strategies agree on mean ─────────────────────────────────── + +// Reference mean computed from 200k plain-MC samples (≈ undiscounted E[max(S_T-K,0)]) +static double plain_mc_ref(std::shared_ptr payoff, size_t n = 200'000) { + auto sampler = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto s = sampler(n, detail::make_seed(999, 999)); + double sum = 0.0; + for (double x : s) sum += x; + return sum / static_cast(s.size()); +} + +TEST_CASE("strategies: all 10 strategies are unbiased for European call") { + auto payoff = std::make_shared(K); + const double ref = plain_mc_ref(payoff); + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + + for (const auto& [name, sampler] : strategies) { + auto s = sampler(30'000, detail::make_seed(42, 1)); + auto [mean, var] = stats(s); + const double se = std::sqrt(var / s.size()); + INFO("Strategy: " << name << " mean=" << mean << " ref=" << ref << " 5-sigma=" << 5.0 * se); + // 5-sigma bound: P(fail) < 0.00006% + REQUIRE(std::abs(mean - ref) < 5.0 * se); + } +} + +TEST_CASE("strategies: importance sampling is unbiased across theta values") { + auto payoff = std::make_shared(K); + const double ref = plain_mc_ref(payoff); + + // theta < 0 shifts paths downward for a call, dramatically increasing variance + // (most paths go OTM, rare ITM paths get huge weights → huge std, tiny mean). + // Only test non-negative theta values where the variance remains manageable. + for (double theta : {0.0, 0.5, 1.0}) { + auto sampler = make_importance_sampler(GBM, EULER, payoff, STEPS, T, theta); + auto s = sampler(30'000, detail::make_seed(7, 2)); + auto [mean, var] = stats(s); + const double se = std::sqrt(var / s.size()); + INFO("theta=" << theta << " mean=" << mean << " ref=" << ref); + REQUIRE(std::abs(mean - ref) < 5.0 * se); + } +} + +// ─── Variance reduction: each strategy beats plain MC for European call ──────── + +// Helper: run both sampler and plain MC for n samples, return var_strategy/var_plain. +// n should be large enough (>=50k) for strategies with small VR (~2%) to be reliable. +static double relative_variance(const SamplerFn& strategy, const SamplerFn& plain_mc, size_t n) { + auto [mean_s, var_s] = stats(strategy(n, detail::make_seed(1, 0))); + auto [mean_p, var_p] = stats(plain_mc(n, detail::make_seed(2, 0))); + (void)mean_s; (void)mean_p; + return var_s / var_p; +} + +// For European call with 52 steps: +// - Antithetic, Halton, stratified-antithetic: large VR (>20%), easily detectable +// - Stratified (1st dim), moment-matching, LHS: ~2% VR from 52 dimensions +// Noise at n=80k is ~0.5% relative std — too close to detect reliably with REQUIRE <1.0 +// So we verify these "marginal" strategies don't increase variance (rv < 1.05), +// and rely on the ASVR option-type tests to demonstrate aggregate benefit. + +TEST_CASE("strategies: stratified sampling has lower variance than plain MC") { + auto payoff = std::make_shared(K); + auto stratified = make_stratified_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + double rv = relative_variance(stratified, plain, 80'000); + INFO("Variance ratio (stratified/plain): " << rv); + // ~2% theoretical VR for the first of 52 dimensions; bound to 5% noise tolerance + REQUIRE(rv < 1.05); +} + +TEST_CASE("strategies: Halton QMC has lower variance than plain MC") { + auto payoff = std::make_shared(K); + auto halton = make_halton_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + double rv = relative_variance(halton, plain, 80'000); + INFO("Variance ratio (halton/plain): " << rv); + REQUIRE(rv < 1.0); +} + +TEST_CASE("strategies: moment matching has lower variance than plain MC") { + auto payoff = std::make_shared(K); + auto mm = make_moment_matching_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + double rv = relative_variance(mm, plain, 80'000); + INFO("Variance ratio (moment_matching/plain): " << rv); + // ~2% theoretical VR from matching all 52 step moments simultaneously + REQUIRE(rv < 1.05); +} + +TEST_CASE("strategies: LHS has lower variance than plain MC") { + auto payoff = std::make_shared(K); + auto lhs = make_lhs_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + double rv = relative_variance(lhs, plain, 80'000); + INFO("Variance ratio (lhs/plain): " << rv); + // ~2% theoretical VR from marginal stratification across 52 dimensions + REQUIRE(rv < 1.05); +} + +TEST_CASE("strategies: stratified antithetic has lower variance than plain MC") { + auto payoff = std::make_shared(K); + auto sa = make_stratified_antithetic_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + double rv = relative_variance(sa, plain, 80'000); + INFO("Variance ratio (stratified_antithetic/plain): " << rv); + REQUIRE(rv < 1.0); +} + +// ─── ASVR over all option types ─────────────────────────────────────────────── + +static ASVRResult run_asvr_for(std::shared_ptr payoff, size_t n = 20'000) { + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + return AdaptiveVarianceReduction::run(std::move(strategies), n, cfg); +} + +static void log_result(const ASVRResult& r) { + for (size_t k = 0; k < r.exploration_stats.size(); ++k) { + INFO(" " << r.exploration_stats[k].name + << ": var=" << r.exploration_stats[k].sample_variance + << " weight=" << r.precision_weights[k] + << " n_exploit=" << r.exploitation_counts[k]); + } + INFO("VR ratio: " << r.variance_reduction_ratio); + INFO("Best strategy: " << r.exploration_stats[ + std::max_element(r.precision_weights.begin(), r.precision_weights.end()) + - r.precision_weights.begin()].name); +} + +TEST_CASE("ASVR with 10 strategies: European call VR ratio > 1") { + auto result = run_asvr_for(std::make_shared(K)); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: European put VR ratio > 1") { + auto result = run_asvr_for(std::make_shared(K)); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: Asian call VR ratio > 1") { + auto result = run_asvr_for(std::make_shared(K)); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: Asian put VR ratio > 1") { + auto result = run_asvr_for(std::make_shared(K)); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: barrier UpAndOutCall VR ratio > 1") { + // Barrier at 130 (30% OTM). Paths that breach 130 pay zero. + auto result = run_asvr_for(std::make_shared(K, BARRIER)); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: barrier DownAndInPut VR ratio > 1") { + // Barrier at 70 (30% below spot). Put only activates if path touches 70. + auto result = run_asvr_for(std::make_shared(K, 70.0)); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: LookbackCall VR ratio > 1") { + auto result = run_asvr_for(std::make_shared()); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +TEST_CASE("ASVR with 10 strategies: LookbackPut VR ratio > 1") { + auto result = run_asvr_for(std::make_shared()); + log_result(result); + REQUIRE(result.variance_reduction_ratio > 1.0); +} + +// ─── Best strategy differs by option type ──────────────────────────────────── +// This test captures the paper's central claim: ASVR automatically selects +// different optimal strategies for fundamentally different option types. + +TEST_CASE("ASVR: best strategy for European call vs lookback call", "[.slow]") { + // Run with large budget for stable weight estimates + auto r_eur = run_asvr_for(std::make_shared(K), 100'000); + auto r_lb = run_asvr_for(std::make_shared(), 100'000); + + const size_t best_eur = static_cast( + std::max_element(r_eur.precision_weights.begin(), r_eur.precision_weights.end()) + - r_eur.precision_weights.begin()); + const size_t best_lb = static_cast( + std::max_element(r_lb.precision_weights.begin(), r_lb.precision_weights.end()) + - r_lb.precision_weights.begin()); + + INFO("European call best: " << r_eur.exploration_stats[best_eur].name); + INFO("Lookback call best: " << r_lb.exploration_stats[best_lb].name); + REQUIRE(best_eur != best_lb); +} + +TEST_CASE("ASVR: best strategy for Asian call vs barrier call", "[.slow]") { + auto r_asian = run_asvr_for(std::make_shared(K), 100'000); + auto r_barrier = run_asvr_for(std::make_shared(K, BARRIER), 100'000); + + const size_t best_asian = static_cast( + std::max_element(r_asian.precision_weights.begin(), r_asian.precision_weights.end()) + - r_asian.precision_weights.begin()); + const size_t best_barrier = static_cast( + std::max_element(r_barrier.precision_weights.begin(), r_barrier.precision_weights.end()) + - r_barrier.precision_weights.begin()); + + INFO("Asian call best: " << r_asian.exploration_stats[best_asian].name); + INFO("Barrier call best: " << r_barrier.exploration_stats[best_barrier].name); + REQUIRE(best_asian != best_barrier); +} + +// ─── Empirical MSE comparison across all option types ───────────────────────── +// Paper table: for each option type, compare ASVR MSE vs plain MC MSE with +// the same total sample budget. This is the core empirical contribution. + +TEST_CASE("strategies: empirical MSE comparison across all option types", "[.slow]") { + const size_t N = 10'000; + const size_t RUNS = 100; + + struct OptionConfig { + std::string name; + std::shared_ptr payoff; + }; + + std::vector options = { + {"European Call", std::make_shared(K)}, + {"European Put", std::make_shared(K)}, + {"Asian Call", std::make_shared(K)}, + {"Asian Put", std::make_shared(K)}, + {"Up-And-Out Call (barrier=130)", std::make_shared(K, BARRIER)}, + {"Down-And-In Put (barrier=70)", std::make_shared(K, 70.0)}, + {"Lookback Call", std::make_shared()}, + {"Lookback Put", std::make_shared()}, + }; + + INFO("\n| Option Type | Plain MC MSE | ASVR MSE | VR Ratio | Best Strategy |"); + INFO( "|---------------------------------|-------------|-----------|----------|---------------------|"); + + auto plain_mc = make_plain_mc_sampler(GBM, EULER, nullptr, STEPS, T); // will be rebuilt per option + + for (auto& [opt_name, payoff] : options) { + // Reference value from 300k plain MC samples + auto plain_ref_fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto ref_samples = plain_ref_fn(300'000, detail::make_seed(0, 77)); + double ref_mean = 0.0; + for (double x : ref_samples) ref_mean += x; + ref_mean /= static_cast(ref_samples.size()); + + auto plain_fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + + double mse_asvr = 0.0; + double mse_plain = 0.0; + std::string best_name; + + ASVRConfig cfg; + cfg.exploration_fraction = 0.1; + cfg.n_threads = 4; + + for (size_t r = 0; r < RUNS; ++r) { + // ASVR + auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); + auto asvr_res = AdaptiveVarianceReduction::run(std::move(strategies), N, cfg); + double err_a = asvr_res.estimate - ref_mean; + mse_asvr += err_a * err_a; + + if (r == 0) { + const size_t best_k = static_cast( + std::max_element(asvr_res.precision_weights.begin(), + asvr_res.precision_weights.end()) + - asvr_res.precision_weights.begin()); + best_name = asvr_res.exploration_stats[best_k].name; + } + + // Plain MC + auto plain_samples = plain_fn(N, detail::make_seed(r + 5000, 88)); + double pm = 0.0; + for (double x : plain_samples) pm += x; + pm /= static_cast(plain_samples.size()); + double err_p = pm - ref_mean; + mse_plain += err_p * err_p; + } + + mse_asvr /= static_cast(RUNS); + mse_plain /= static_cast(RUNS); + const double vr = mse_plain / mse_asvr; + + INFO("| " << opt_name + << " | " << mse_plain + << " | " << mse_asvr + << " | " << vr + << " | " << best_name << " |"); + + REQUIRE(mse_asvr < mse_plain); + } +} From 7d7a0d98e58852abcf11d889c42e80487ed41786 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 31 Mar 2026 14:01:48 +0000 Subject: [PATCH 12/12] style: apply clang-format --- include/lfmc/adaptive_estimator.hpp | 19 ++-- include/lfmc/engine.hpp | 110 +++++++++++----------- include/lfmc/strategies.hpp | 86 +++++++++-------- src/adaptive/adaptive_estimator.cpp | 44 ++++----- src/payoff/asian_payoffs.cpp | 6 +- tests/bench_asvr_vs_fixed.cpp | 87 ++++++++++-------- tests/quick_compare.cpp | 29 +++--- tests/test_adaptive.cpp | 17 ++-- tests/test_engine.cpp | 76 ++++++++------- tests/test_strategies.cpp | 137 +++++++++++++++------------- 10 files changed, 313 insertions(+), 298 deletions(-) diff --git a/include/lfmc/adaptive_estimator.hpp b/include/lfmc/adaptive_estimator.hpp index cde392b..9b6679d 100644 --- a/include/lfmc/adaptive_estimator.hpp +++ b/include/lfmc/adaptive_estimator.hpp @@ -101,9 +101,8 @@ namespace detail { // Two calls with the same (k, p) give the same seed; different (k, p) pairs are // statistically independent because they hit different streams of splitmix64. inline uint64_t make_seed(size_t k, size_t phase) noexcept { - uint64_t x = (static_cast(k) * 0x9e3779b97f4a7c15ULL) - ^ (static_cast(phase) * 0x6c62272e07bb0142ULL) - ^ 0xdeadbeefcafe0000ULL; + uint64_t x = (static_cast(k) * 0x9e3779b97f4a7c15ULL) ^ + (static_cast(phase) * 0x6c62272e07bb0142ULL) ^ 0xdeadbeefcafe0000ULL; x ^= x >> 30; x *= 0xbf58476d1ce4e5b9ULL; x ^= x >> 27; @@ -169,8 +168,8 @@ inline std::pair mean_variance(const std::vector& sample // Plain pseudo-random Monte Carlo — baseline strategy template NS> -SamplerFn make_plain_mc_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { +SamplerFn make_plain_mc_sampler(SP process, NS scheme, std::shared_ptr payoff, size_t steps, + double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::vector samples; @@ -191,7 +190,7 @@ SamplerFn make_plain_mc_sampler(SP process, NS scheme, std::shared_ptr p // because the two paths are negatively correlated for monotone payoffs. template NS> SamplerFn make_antithetic_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { + size_t steps, double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::vector samples; @@ -223,8 +222,8 @@ SamplerFn make_antithetic_sampler(SP process, NS scheme, std::shared_ptr // Returns: X_i - beta_hat * (Y_i - E[Y]) for each i template NS> SamplerFn make_control_variate_sampler(SP process, NS scheme, std::shared_ptr target, - std::shared_ptr control, double control_mean, - size_t steps, double T) { + std::shared_ptr control, double control_mean, + size_t steps, double T) { return [process, scheme, target, control, control_mean, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; @@ -272,8 +271,8 @@ SamplerFn make_control_variate_sampler(SP process, NS scheme, std::shared_ptr NS> SamplerFn make_antithetic_cv_sampler(SP process, NS scheme, std::shared_ptr target, - std::shared_ptr control, double control_mean, - size_t steps, double T) { + std::shared_ptr control, double control_mean, + size_t steps, double T) { return [process, scheme, target, control, control_mean, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; diff --git a/include/lfmc/engine.hpp b/include/lfmc/engine.hpp index e79063b..c7b37c3 100644 --- a/include/lfmc/engine.hpp +++ b/include/lfmc/engine.hpp @@ -3,8 +3,8 @@ #include "lfmc/adaptive_estimator.hpp" #include "lfmc/numerical_scheme.hpp" #include "lfmc/payoff.hpp" -#include "lfmc/strategies.hpp" #include "lfmc/stochastic_process.hpp" +#include "lfmc/strategies.hpp" #include "lfmc/types.hpp" #include @@ -52,38 +52,36 @@ struct IterativeEngineConfig { }; struct RoundResult { - size_t round; - std::string leader; // strategy leading after this round - bool leader_changed; // did the leader switch this round? - double estimate; // precision-weighted estimate so far - double estimated_stderr; + size_t round; + std::string leader; // strategy leading after this round + bool leader_changed; // did the leader switch this round? + double estimate; // precision-weighted estimate so far + double estimated_stderr; std::vector precision_weights; // one per strategy, after this round }; struct IterativeEngineResult { - double estimate; - double estimated_stderr; + double estimate; + double estimated_stderr; std::string final_leader; - std::vector round_history; // diagnostics per round - std::vector final_stats; // cumulative stats per strategy - size_t total_samples; + std::vector round_history; // diagnostics per round + std::vector final_stats; // cumulative stats per strategy + size_t total_samples; }; -template NS> -class Engine { - SP process_; - NS scheme_; +template NS> class Engine { + SP process_; + NS scheme_; double control_mean_; size_t steps_; double T_; public: Engine(SP process, NS scheme, double control_mean, size_t steps, double T) - : process_(process), scheme_(scheme), control_mean_(control_mean), - steps_(steps), T_(T) {} + : process_(process), scheme_(scheme), control_mean_(control_mean), steps_(steps), T_(T) {} - std::expected - run(std::shared_ptr payoff, IterativeEngineConfig config = {}) { + std::expected run(std::shared_ptr payoff, + IterativeEngineConfig config = {}) { auto all = build_all_strategies(process_, scheme_, payoff, control_mean_, steps_, T_); const size_t K = std::min(config.n_compete, all.size()); if (K == 0) @@ -93,13 +91,13 @@ class Engine { // We track (sum, sum_sq, n) across all rounds so precision weights can // be recomputed from the full history after each round. struct Accumulator { - double sum = 0.0; + double sum = 0.0; double sum_sq = 0.0; - size_t n = 0; + size_t n = 0; void add(const std::vector& samples) { for (double x : samples) { - sum += x; + sum += x; sum_sq += x * x; ++n; } @@ -111,10 +109,10 @@ class Engine { // Unbiased sample variance (n-1 denominator) double variance() const { - if (n < 2) return std::numeric_limits::infinity(); + if (n < 2) + return std::numeric_limits::infinity(); const double m = mean(); - return (sum_sq - static_cast(n) * m * m) / - static_cast(n - 1); + return (sum_sq - static_cast(n) * m * m) / static_cast(n - 1); } }; @@ -127,7 +125,7 @@ class Engine { for (size_t k = 0; k < K; ++k) { double var = accum[k].variance(); prec[k] = (var > 0.0 && std::isfinite(var)) ? 1.0 / var : 0.0; - total += prec[k]; + total += prec[k]; } std::vector w(K); if (total > 0.0) { @@ -139,8 +137,8 @@ class Engine { return w; }; - size_t leader = 0; // index of current leader - bool first_round = true; + size_t leader = 0; // index of current leader + bool first_round = true; std::vector history; history.reserve(config.n_rounds); @@ -157,7 +155,7 @@ class Engine { std::vector thread_counts(K, 1); // 1 compete thread each if (first_round) { const size_t bonus_each = config.n_exploit / K; - const size_t leftover = config.n_exploit % K; + const size_t leftover = config.n_exploit % K; for (size_t k = 0; k < K; ++k) thread_counts[k] += bonus_each; thread_counts[0] += leftover; @@ -174,7 +172,10 @@ class Engine { std::vector> round_samples(K); { // Flatten to (strategy_idx, thread_idx) pairs for easy launch - struct Job { size_t k; size_t t; }; + struct Job { + size_t k; + size_t t; + }; std::vector jobs; for (size_t k = 0; k < K; ++k) for (size_t t = 0; t < thread_counts[k]; ++t) @@ -186,10 +187,10 @@ class Engine { threads.reserve(jobs.size()); for (size_t j = 0; j < jobs.size(); ++j) { threads.emplace_back([&, j]() { - const size_t k = jobs[j].k; - const size_t t = jobs[j].t; + const size_t k = jobs[j].k; + const size_t t = jobs[j].t; const uint64_t seed = detail::make_seed(k * 10000 + t, round); - thread_results[j] = all[k].second(config.samples_per_thread, seed); + thread_results[j] = all[k].second(config.samples_per_thread, seed); }); } } // jthreads join here @@ -197,7 +198,7 @@ class Engine { // Merge thread results into per-strategy buckets for (size_t j = 0; j < jobs.size(); ++j) { auto& dest = round_samples[jobs[j].k]; - auto& src = thread_results[j]; + auto& src = thread_results[j]; dest.insert(dest.end(), src.begin(), src.end()); } } @@ -209,7 +210,7 @@ class Engine { } // ── Recompute precision weights and current leader ───────────────── - auto weights = precision_weights(); + auto weights = precision_weights(); size_t new_leader = static_cast( std::max_element(weights.begin(), weights.end()) - weights.begin()); @@ -223,26 +224,25 @@ class Engine { for (size_t k = 0; k < K; ++k) { double var_k = accum[k].variance(); if (std::isfinite(var_k) && accum[k].n > 0) - var_est += weights[k] * weights[k] * var_k / - static_cast(accum[k].n); + var_est += weights[k] * weights[k] * var_k / static_cast(accum[k].n); } history.push_back({ - .round = round + 1, - .leader = all[new_leader].first, - .leader_changed = (!first_round && new_leader != leader), - .estimate = est, + .round = round + 1, + .leader = all[new_leader].first, + .leader_changed = (!first_round && new_leader != leader), + .estimate = est, .estimated_stderr = std::sqrt(var_est), .precision_weights = weights, }); - leader = new_leader; + leader = new_leader; first_round = false; } // ── Final estimate ──────────────────────────────────────────────────── auto final_weights = precision_weights(); - double final_est = 0.0; + double final_est = 0.0; for (size_t k = 0; k < K; ++k) final_est += final_weights[k] * accum[k].mean(); @@ -250,29 +250,29 @@ class Engine { for (size_t k = 0; k < K; ++k) { double var_k = accum[k].variance(); if (std::isfinite(var_k) && accum[k].n > 0) - final_var += final_weights[k] * final_weights[k] * var_k / - static_cast(accum[k].n); + final_var += + final_weights[k] * final_weights[k] * var_k / static_cast(accum[k].n); } // ── Pack cumulative StrategyStats for diagnostics ────────────────────── std::vector final_stats(K); for (size_t k = 0; k < K; ++k) { final_stats[k] = { - .name = all[k].first, - .mean = accum[k].mean(), + .name = all[k].first, + .mean = accum[k].mean(), .sample_variance = accum[k].variance(), - .wall_time_ms = 0.0, // not tracked per-strategy across rounds - .n_samples = accum[k].n, + .wall_time_ms = 0.0, // not tracked per-strategy across rounds + .n_samples = accum[k].n, }; } return IterativeEngineResult{ - .estimate = final_est, - .estimated_stderr = std::sqrt(final_var), - .final_leader = all[leader].first, - .round_history = history, - .final_stats = final_stats, - .total_samples = total_samples, + .estimate = final_est, + .estimated_stderr = std::sqrt(final_var), + .final_leader = all[leader].first, + .round_history = history, + .final_stats = final_stats, + .total_samples = total_samples, }; } }; diff --git a/include/lfmc/strategies.hpp b/include/lfmc/strategies.hpp index 0549148..5debb72 100644 --- a/include/lfmc/strategies.hpp +++ b/include/lfmc/strategies.hpp @@ -50,7 +50,7 @@ inline double normal_icdf(double p) noexcept { } const double q = std::sqrt(-2.0 * std::log(1.0 - p)); return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / - ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0); + ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0); } // ─── Halton radical inverse ─────────────────────────────────────────────────── @@ -100,7 +100,7 @@ namespace lfmc { // each stratum and converted to N(0,1) via normal_icdf. template NS> SamplerFn make_stratified_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { + size_t steps, double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::uniform_real_distribution udist{0.0, 1.0}; @@ -139,8 +139,8 @@ SamplerFn make_stratified_sampler(SP process, NS scheme, std::shared_ptr // large dimension counts (steps > ~20). ASVR will automatically assign low // weight to this strategy for path-dependent options where all steps matter. template NS> -SamplerFn make_halton_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { +SamplerFn make_halton_sampler(SP process, NS scheme, std::shared_ptr payoff, size_t steps, + double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::uniform_real_distribution udist{0.0, 1.0}; @@ -195,38 +195,38 @@ SamplerFn make_halton_sampler(SP process, NS scheme, std::shared_ptr pay // exp(-0.5*theta^2) shrinks too fast and variance increases template NS> SamplerFn make_importance_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T, double theta) { - return [process, scheme, payoff, steps, T, theta](size_t n, - uint64_t seed) -> std::vector { - std::mt19937_64 rng{seed}; - // Per-step shift: distribute the total theta equally across sqrt(steps) normal draws - const double delta = theta / std::sqrt(static_cast(steps)); - std::normal_distribution ndist{delta, 1.0}; // tilted draw + size_t steps, double T, double theta) { + return + [process, scheme, payoff, steps, T, theta](size_t n, uint64_t seed) -> std::vector { + std::mt19937_64 rng{seed}; + // Per-step shift: distribute the total theta equally across sqrt(steps) normal draws + const double delta = theta / std::sqrt(static_cast(steps)); + std::normal_distribution ndist{delta, 1.0}; // tilted draw - // Log-weight constant: 0.5 * delta^2 * steps = 0.5 * theta^2 - const double lw_const = 0.5 * theta * theta; + // Log-weight constant: 0.5 * delta^2 * steps = 0.5 * theta^2 + const double lw_const = 0.5 * theta * theta; - std::vector samples; - samples.reserve(n); + std::vector samples; + samples.reserve(n); - for (size_t i = 0; i < n; ++i) { - Normals normals(steps); - double z_sum = 0.0; - for (size_t j = 0; j < steps; ++j) { - normals[j] = ndist(rng); - z_sum += normals[j]; + for (size_t i = 0; i < n; ++i) { + Normals normals(steps); + double z_sum = 0.0; + for (size_t j = 0; j < steps; ++j) { + normals[j] = ndist(rng); + z_sum += normals[j]; + } + // Likelihood ratio dP/dQ = exp(-delta * sum(Z') + 0.5 * delta^2 * steps) + // = exp(-theta/sqrt(steps) * sum(Z') + 0.5 * theta^2) + const double weight = std::exp(-delta * z_sum + lw_const); + + auto path = detail::generate_path(process, scheme, normals, steps, T); + auto result = payoff->generate_payoffs({path}); + if (result && !result->empty() && !(*result)[0].empty()) + samples.push_back((*result)[0][0] * weight); } - // Likelihood ratio dP/dQ = exp(-delta * sum(Z') + 0.5 * delta^2 * steps) - // = exp(-theta/sqrt(steps) * sum(Z') + 0.5 * theta^2) - const double weight = std::exp(-delta * z_sum + lw_const); - - auto path = detail::generate_path(process, scheme, normals, steps, T); - auto result = payoff->generate_payoffs({path}); - if (result && !result->empty() && !(*result)[0].empty()) - samples.push_back((*result)[0][0] * weight); - } - return samples; - }; + return samples; + }; } // ─── Strategy 8: Moment matching ───────────────────────────────────────────── @@ -238,7 +238,7 @@ SamplerFn make_importance_sampler(SP process, NS scheme, std::shared_ptr // Memory: O(n * steps) doubles. For n=10,000 and steps=52 this is ~4 MB. template NS> SamplerFn make_moment_matching_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { + size_t steps, double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::normal_distribution ndist{0.0, 1.0}; @@ -293,8 +293,8 @@ SamplerFn make_moment_matching_sampler(SP process, NS scheme, std::shared_ptr NS> -SamplerFn make_lhs_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { +SamplerFn make_lhs_sampler(SP process, NS scheme, std::shared_ptr payoff, size_t steps, + double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::uniform_real_distribution udist{0.0, 1.0}; @@ -346,7 +346,7 @@ SamplerFn make_lhs_sampler(SP process, NS scheme, std::shared_ptr payoff // correlated higher dimensions within the pair. template NS> SamplerFn make_stratified_antithetic_sampler(SP process, NS scheme, std::shared_ptr payoff, - size_t steps, double T) { + size_t steps, double T) { return [process, scheme, payoff, steps, T](size_t n, uint64_t seed) -> std::vector { std::mt19937_64 rng{seed}; std::uniform_real_distribution udist{0.0, 1.0}; @@ -376,8 +376,8 @@ SamplerFn make_stratified_antithetic_sampler(SP process, NS scheme, std::shared_ auto r_pos = payoff->generate_payoffs({path_pos}); auto r_neg = payoff->generate_payoffs({path_neg}); - if (r_pos && r_neg && !r_pos->empty() && !r_neg->empty() && - !(*r_pos)[0].empty() && !(*r_neg)[0].empty()) + if (r_pos && r_neg && !r_pos->empty() && !r_neg->empty() && !(*r_pos)[0].empty() && + !(*r_neg)[0].empty()) samples.push_back(0.5 * ((*r_pos)[0][0] + (*r_neg)[0][0])); } return samples; @@ -400,12 +400,10 @@ build_all_strategies(SP process, NS scheme, std::shared_ptr payoff, doub return { {"plain_mc", make_plain_mc_sampler(process, scheme, payoff, steps, T)}, {"antithetic", make_antithetic_sampler(process, scheme, payoff, steps, T)}, - {"control_variate", - make_control_variate_sampler(process, scheme, payoff, control_payoff(), control_mean, - steps, T)}, - {"antithetic_cv", - make_antithetic_cv_sampler(process, scheme, payoff, control_payoff(), control_mean, steps, - T)}, + {"control_variate", make_control_variate_sampler(process, scheme, payoff, control_payoff(), + control_mean, steps, T)}, + {"antithetic_cv", make_antithetic_cv_sampler(process, scheme, payoff, control_payoff(), + control_mean, steps, T)}, {"stratified", make_stratified_sampler(process, scheme, payoff, steps, T)}, {"halton_qmc", make_halton_sampler(process, scheme, payoff, steps, T)}, {"importance_sampling", diff --git a/src/adaptive/adaptive_estimator.cpp b/src/adaptive/adaptive_estimator.cpp index 8a1b86f..582f20a 100644 --- a/src/adaptive/adaptive_estimator.cpp +++ b/src/adaptive/adaptive_estimator.cpp @@ -11,8 +11,8 @@ namespace lfmc { // unbiased sample variance. StrategyStats AdaptiveVarianceReduction::run_strategy_batch(const std::string& name, - const SamplerFn& sampler, - size_t n_samples, uint64_t seed) { + const SamplerFn& sampler, + size_t n_samples, uint64_t seed) { auto t0 = std::chrono::high_resolution_clock::now(); auto samples = sampler(n_samples, seed); auto t1 = std::chrono::high_resolution_clock::now(); @@ -57,17 +57,17 @@ AdaptiveVarianceReduction::compute_precision_weights(const std::vector> strategies, - size_t total_samples, ASVRConfig config) { + size_t total_samples, ASVRConfig config) { const size_t K = strategies.size(); assert(K >= 1); const size_t hw = std::max(size_t{1}, config.n_threads); // ── Budget split ───────────────────────────────────────────────────────── - const size_t n_explore_each = std::max( - config.min_exploration_per_strategy, - static_cast(config.exploration_fraction * static_cast(total_samples) / - static_cast(K))); + const size_t n_explore_each = + std::max(config.min_exploration_per_strategy, + static_cast(config.exploration_fraction * + static_cast(total_samples) / static_cast(K))); const size_t n_explore_total = K * n_explore_each; const size_t n_exploit = (n_explore_total < total_samples) ? total_samples - n_explore_total : 0; @@ -77,10 +77,9 @@ ASVRResult AdaptiveVarianceReduction::run(std::vector(config.exploration_fraction * - static_cast(hw)))); + const size_t explore_thread_cap = std::max( + size_t{1}, + std::min(K, static_cast(config.exploration_fraction * static_cast(hw)))); std::vector exploration_stats(K); @@ -91,9 +90,8 @@ ASVRResult AdaptiveVarianceReduction::run(std::vector 0) { size_t assigned = 0; for (size_t k = 0; k < K; ++k) { - exploit_counts[k] = - static_cast(weights[k] * static_cast(n_exploit)); + exploit_counts[k] = static_cast(weights[k] * static_cast(n_exploit)); assigned += exploit_counts[k]; } if (assigned < n_exploit) { - const size_t best_k = - static_cast(std::max_element(weights.begin(), weights.end()) - - weights.begin()); + const size_t best_k = static_cast( + std::max_element(weights.begin(), weights.end()) - weights.begin()); exploit_counts[best_k] += (n_exploit - assigned); } } @@ -139,9 +135,9 @@ ASVRResult AdaptiveVarianceReduction::run(std::vector 0) ? var_exploit - : exploration_stats[0].sample_variance / - static_cast(n_explore_each); + : exploration_stats[0].sample_variance / + static_cast(n_explore_each); // Plain MC baseline for comparison: sigma_0^2 / N using the first strategy // (conventionally plain MC; the caller must ensure strategies[0] is plain MC) diff --git a/src/payoff/asian_payoffs.cpp b/src/payoff/asian_payoffs.cpp index 749a5e9..a441ea6 100644 --- a/src/payoff/asian_payoffs.cpp +++ b/src/payoff/asian_payoffs.cpp @@ -20,8 +20,7 @@ AsianCall::generate_payoffs(const std::vector& paths) const { if (path.empty()) return std::unexpected("Empty path encountered in AsianCall"); - double mean = - std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); + double mean = std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); payoffs.push_back(std::max(mean - strike_, 0.0)); } @@ -39,8 +38,7 @@ AsianPut::generate_payoffs(const std::vector& paths) const { if (path.empty()) return std::unexpected("Empty path encountered in AsianPut"); - double mean = - std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); + double mean = std::reduce(path.begin(), path.end(), 0.0) / static_cast(path.size()); payoffs.push_back(std::max(strike_ - mean, 0.0)); } diff --git a/tests/bench_asvr_vs_fixed.cpp b/tests/bench_asvr_vs_fixed.cpp index d061a7a..c14fc90 100644 --- a/tests/bench_asvr_vs_fixed.cpp +++ b/tests/bench_asvr_vs_fixed.cpp @@ -24,12 +24,12 @@ using namespace lfmc; -static constexpr double S0 = 100.0; -static constexpr double MU = 0.05; -static constexpr double SIGMA = 0.2; -static constexpr double K = 100.0; -static constexpr double T = 1.0; -static constexpr size_t STEPS = 52; +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; static const double CONTROL_MEAN = S0 * std::exp(MU * T); static GeometricBrownianMotion GBM{MU, SIGMA, S0}; @@ -43,8 +43,8 @@ struct BenchResult { double std_error_of_mean; }; -BenchResult bench_fixed(const std::string& name, SamplerFn sampler, - double ref_mean, size_t N, size_t RUNS) { +BenchResult bench_fixed(const std::string& name, SamplerFn sampler, double ref_mean, size_t N, + size_t RUNS) { std::vector errors; errors.reserve(RUNS); double sum_est = 0.0; @@ -60,13 +60,15 @@ BenchResult bench_fixed(const std::string& name, SamplerFn sampler, // std error of the MSE estimate (rough) double mean_err2 = mse; double sq_sum = 0.0; - for (double e2 : errors) { double d = e2 - mean_err2; sq_sum += d*d; } - double sem = std::sqrt(sq_sum / static_cast(RUNS * (RUNS-1))); + for (double e2 : errors) { + double d = e2 - mean_err2; + sq_sum += d * d; + } + double sem = std::sqrt(sq_sum / static_cast(RUNS * (RUNS - 1))); return {name, mse, sum_est / RUNS, sem}; } -BenchResult bench_asvr(std::shared_ptr payoff, double ref_mean, - size_t N, size_t RUNS) { +BenchResult bench_asvr(std::shared_ptr payoff, double ref_mean, size_t N, size_t RUNS) { std::vector errors; errors.reserve(RUNS); double sum_est = 0.0; @@ -83,20 +85,23 @@ BenchResult bench_asvr(std::shared_ptr payoff, double ref_mean, double mse = std::accumulate(errors.begin(), errors.end(), 0.0) / static_cast(RUNS); double mean_err2 = mse; double sq_sum = 0.0; - for (double e2 : errors) { double d = e2 - mean_err2; sq_sum += d*d; } - double sem = std::sqrt(sq_sum / static_cast(RUNS * (RUNS-1))); + for (double e2 : errors) { + double d = e2 - mean_err2; + sq_sum += d * d; + } + double sem = std::sqrt(sq_sum / static_cast(RUNS * (RUNS - 1))); return {"ASVR", mse, sum_est / RUNS, sem}; } -void run_option_bench(const std::string& opt_name, std::shared_ptr payoff, - size_t N, size_t RUNS) { +void run_option_bench(const std::string& opt_name, std::shared_ptr payoff, size_t N, + size_t RUNS) { printf("\n=== %-35s (N=%zu, runs=%zu) ===\n", opt_name.c_str(), N, RUNS); // Reference from 500k plain MC samples auto plain_ref_fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); auto ref_samples = plain_ref_fn(500'000, detail::make_seed(0, 99)); - double ref_mean = std::accumulate(ref_samples.begin(), ref_samples.end(), 0.0) - / static_cast(ref_samples.size()); + double ref_mean = std::accumulate(ref_samples.begin(), ref_samples.end(), 0.0) / + static_cast(ref_samples.size()); printf(" Reference mean: %.6f\n", ref_mean); auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); @@ -111,47 +116,49 @@ void run_option_bench(const std::string& opt_name, std::shared_ptr payof results.push_back(bench_asvr(payoff, ref_mean, N, RUNS)); // Sort by MSE ascending - std::sort(results.begin(), results.end(), - [](const BenchResult& a, const BenchResult& b) { - return a.mean_sq_error < b.mean_sq_error; - }); + std::sort(results.begin(), results.end(), [](const BenchResult& a, const BenchResult& b) { + return a.mean_sq_error < b.mean_sq_error; + }); double plain_mse = 0.0; - for (auto& r : results) if (r.name == "plain_mc") plain_mse = r.mean_sq_error; + for (auto& r : results) + if (r.name == "plain_mc") + plain_mse = r.mean_sq_error; printf(" %-28s %12s %8s %8s\n", "Strategy", "MSE", "vs PlainMC", "vs ASVR"); - printf(" %-28s %12s %8s %8s\n", - std::string(28,'-').c_str(), std::string(12,'-').c_str(), - std::string(8,'-').c_str(), std::string(8,'-').c_str()); + printf(" %-28s %12s %8s %8s\n", std::string(28, '-').c_str(), std::string(12, '-').c_str(), + std::string(8, '-').c_str(), std::string(8, '-').c_str()); double asvr_mse = 0.0; - for (auto& r : results) if (r.name == "ASVR") asvr_mse = r.mean_sq_error; + for (auto& r : results) + if (r.name == "ASVR") + asvr_mse = r.mean_sq_error; for (auto& r : results) { double vs_plain = plain_mse / r.mean_sq_error; - double vs_asvr = asvr_mse / r.mean_sq_error; - printf(" %-28s %12.6f %8.3fx %8.3fx\n", - r.name.c_str(), r.mean_sq_error, vs_plain, vs_asvr); + double vs_asvr = asvr_mse / r.mean_sq_error; + printf(" %-28s %12.6f %8.3fx %8.3fx\n", r.name.c_str(), r.mean_sq_error, vs_plain, + vs_asvr); } } int main() { - const size_t N = 10'000; // sample budget per trial - const size_t RUNS = 200; // trials for MSE estimate + const size_t N = 10'000; // sample budget per trial + const size_t RUNS = 200; // trials for MSE estimate printf("ASVR vs Fixed Strategy Benchmark\n"); printf("=================================\n"); printf("Budget N=%zu per trial, %zu independent runs\n", N, RUNS); printf("Strategies compared: all 10 fixed + ASVR (adaptive)\n"); - run_option_bench("European Call", std::make_shared(K), N, RUNS); - run_option_bench("European Put", std::make_shared(K), N, RUNS); - run_option_bench("Asian Call", std::make_shared(K), N, RUNS); - run_option_bench("Asian Put", std::make_shared(K), N, RUNS); - run_option_bench("Up-And-Out Call B=130", std::make_shared(K, 130.0), N, RUNS); - run_option_bench("Down-And-In Put B=70", std::make_shared(K, 70.0), N, RUNS); - run_option_bench("Lookback Call", std::make_shared(), N, RUNS); - run_option_bench("Lookback Put", std::make_shared(), N, RUNS); + run_option_bench("European Call", std::make_shared(K), N, RUNS); + run_option_bench("European Put", std::make_shared(K), N, RUNS); + run_option_bench("Asian Call", std::make_shared(K), N, RUNS); + run_option_bench("Asian Put", std::make_shared(K), N, RUNS); + run_option_bench("Up-And-Out Call B=130", std::make_shared(K, 130.0), N, RUNS); + run_option_bench("Down-And-In Put B=70", std::make_shared(K, 70.0), N, RUNS); + run_option_bench("Lookback Call", std::make_shared(), N, RUNS); + run_option_bench("Lookback Put", std::make_shared(), N, RUNS); printf("\nDone.\n"); return 0; diff --git a/tests/quick_compare.cpp b/tests/quick_compare.cpp index b300912..7979bdf 100644 --- a/tests/quick_compare.cpp +++ b/tests/quick_compare.cpp @@ -16,11 +16,11 @@ using namespace lfmc; -static constexpr double S0 = 100.0; -static constexpr double MU = 0.05; +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; static constexpr double SIGMA = 0.2; -static constexpr double K = 100.0; -static constexpr double T = 1.0; +static constexpr double K = 100.0; +static constexpr double T = 1.0; static constexpr size_t STEPS = 52; static const double CONTROL_MEAN = S0 * std::exp(MU * T); @@ -30,13 +30,14 @@ static EulerMaruyama EULER; // Total budget: Engine uses 4*2000 (explore) + 8*5000 (exploit) = 48000 samples. // Fixed strategies get the same 48000 samples for a fair comparison. static constexpr size_t TOTAL = 48000; -static constexpr size_t RUNS = 50; +static constexpr size_t RUNS = 50; static double ref_mean(std::shared_ptr payoff) { auto fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); auto s = fn(500'000, detail::make_seed(0, 77)); double sum = 0.0; - for (double x : s) sum += x; + for (double x : s) + sum += x; return sum / static_cast(s.size()); } @@ -57,14 +58,15 @@ static double bench_engine(std::shared_ptr payoff, double ref, size_t ru EngineConfig cfg; cfg.n_explore_threads = 4; - cfg.explore_samples = 2000; + cfg.explore_samples = 2000; cfg.n_exploit_threads = 8; - cfg.exploit_samples = 5000; + cfg.exploit_samples = 5000; double mse = 0.0; for (size_t r = 0; r < runs; ++r) { auto res = engine.run(payoff, cfg); - if (!res) continue; + if (!res) + continue; double e = res->estimate - ref; mse += e * e; } @@ -87,7 +89,8 @@ void compare(const char* opt_name, std::shared_ptr payoff) { for (auto& [name, fn] : strategies) { double mse = bench_fixed(name, fn, ref, RUNS); rows.push_back({name, mse}); - if (name == "plain_mc") plain_mse = mse; + if (name == "plain_mc") + plain_mse = mse; } double engine_mse = bench_engine(payoff, ref, RUNS); @@ -105,10 +108,10 @@ int main() { printf("Engine vs Fixed Strategy Comparison\n"); printf("Same total sample budget (%zu), %zu independent runs\n\n", TOTAL, RUNS); - compare("European Call", std::make_shared(K)); - compare("Asian Call", std::make_shared(K)); + compare("European Call", std::make_shared(K)); + compare("Asian Call", std::make_shared(K)); compare("Barrier (Up-Out Call B=130)", std::make_shared(K, 130.0)); - compare("Lookback Call", std::make_shared()); + compare("Lookback Call", std::make_shared()); printf("\nDone.\n"); } diff --git a/tests/test_adaptive.cpp b/tests/test_adaptive.cpp index edd980e..25e6c5e 100644 --- a/tests/test_adaptive.cpp +++ b/tests/test_adaptive.cpp @@ -35,11 +35,11 @@ static double bs_call(double S, double K, double r, double sigma, double T) { // // The Black-Scholes discounted price is NOT the right comparison target here. -static constexpr double S0 = 100.0; -static constexpr double MU = 0.05; +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; static constexpr double SIGMA = 0.2; -static constexpr double K = 100.0; -static constexpr double T = 1.0; +static constexpr double K = 100.0; +static constexpr double T = 1.0; static constexpr size_t STEPS = 52; // weekly steps // control: E[S_T] = S0 * exp(mu*T) static const double CONTROL_MEAN = S0 * std::exp(MU * T); @@ -60,12 +60,11 @@ static std::vector> build_strategies() { // for GBM), E[S_T] = S0 * exp(mu*T) {"control_variate", make_control_variate_sampler(gbm, euler, std::make_shared(K), - std::make_shared(0.0), CONTROL_MEAN, STEPS, - T)}, + std::make_shared(0.0), CONTROL_MEAN, STEPS, T)}, {"antithetic_cv", make_antithetic_cv_sampler(gbm, euler, std::make_shared(K), - std::make_shared(0.0), CONTROL_MEAN, STEPS, T)}, + std::make_shared(0.0), CONTROL_MEAN, STEPS, T)}, }; } @@ -136,8 +135,8 @@ TEST_CASE("ASVR: exploration identifies a VR strategy as best") { // Find which strategy got the highest weight const auto& weights = result.precision_weights; - const size_t best_k = static_cast( - std::max_element(weights.begin(), weights.end()) - weights.begin()); + const size_t best_k = + static_cast(std::max_element(weights.begin(), weights.end()) - weights.begin()); INFO("Best strategy: " << result.exploration_stats[best_k].name); INFO("Weights: "); diff --git a/tests/test_engine.cpp b/tests/test_engine.cpp index 39187a8..dbc2d02 100644 --- a/tests/test_engine.cpp +++ b/tests/test_engine.cpp @@ -12,17 +12,17 @@ using namespace lfmc; -static constexpr double S0 = 100.0; -static constexpr double MU = 0.05; +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; static constexpr double SIGMA = 0.2; -static constexpr double K = 100.0; -static constexpr double T = 1.0; +static constexpr double K = 100.0; +static constexpr double T = 1.0; static constexpr size_t STEPS = 52; static const double CONTROL_MEAN = S0 * std::exp(MU * T); static Engine> make_engine() { - return {GeometricBrownianMotion{MU, SIGMA, S0}, - EulerMaruyama{}, CONTROL_MEAN, STEPS, T}; + return {GeometricBrownianMotion{MU, SIGMA, S0}, EulerMaruyama{}, + CONTROL_MEAN, STEPS, T}; } // ─── Smoke tests ────────────────────────────────────────────────────────────── @@ -71,7 +71,8 @@ TEST_CASE("Engine: precision weights sum to 1 after each round") { for (const auto& r : result->round_history) { double sum = 0.0; - for (double w : r.precision_weights) sum += w; + for (double w : r.precision_weights) + sum += w; REQUIRE(sum == Catch::Approx(1.0).epsilon(1e-10)); } } @@ -80,9 +81,9 @@ TEST_CASE("Engine: total samples consistent with config") { auto engine = make_engine(); IterativeEngineConfig cfg; - cfg.n_compete = 4; - cfg.n_exploit = 8; - cfg.n_rounds = 3; + cfg.n_compete = 4; + cfg.n_exploit = 8; + cfg.n_rounds = 3; cfg.samples_per_thread = 500; auto result = engine.run(std::make_shared(K), cfg); @@ -92,9 +93,7 @@ TEST_CASE("Engine: total samples consistent with config") { // + leftover 8%4=0, so 3 per strategy = 12 threads * 500 = 6000 // Rounds 2-3: 4 compete (1 each) + 8 exploit on leader = 12 threads * 500 = 6000 each // Total = 3 * 12 * 500 = 18000 - const size_t expected = cfg.n_rounds * - (cfg.n_compete + cfg.n_exploit) * - cfg.samples_per_thread; + const size_t expected = cfg.n_rounds * (cfg.n_compete + cfg.n_exploit) * cfg.samples_per_thread; REQUIRE(result->total_samples == expected); } @@ -134,9 +133,9 @@ TEST_CASE("Engine: leader stabilises over rounds") { auto engine = make_engine(); IterativeEngineConfig cfg; - cfg.n_compete = 4; - cfg.n_exploit = 8; - cfg.n_rounds = 6; + cfg.n_compete = 4; + cfg.n_exploit = 8; + cfg.n_rounds = 6; cfg.samples_per_thread = 1000; auto result = engine.run(std::make_shared(K), cfg); @@ -146,12 +145,12 @@ TEST_CASE("Engine: leader stabilises over rounds") { size_t late_changes = 0; const auto& h = result->round_history; for (size_t i = h.size() / 2; i < h.size(); ++i) - if (h[i].leader_changed) ++late_changes; + if (h[i].leader_changed) + ++late_changes; WARN("Leader per round:"); for (const auto& r : h) - WARN(" Round " << r.round << ": " << r.leader - << (r.leader_changed ? " (SWITCHED)" : "")); + WARN(" Round " << r.round << ": " << r.leader << (r.leader_changed ? " (SWITCHED)" : "")); // Allow at most 1 late switch — the leader might flip once as estimates refine REQUIRE(late_changes <= 1); @@ -161,18 +160,21 @@ TEST_CASE("Engine: leader stabilises over rounds") { TEST_CASE("Engine: empirical MSE vs fixed strategies", "[bench]") { static constexpr size_t TOTAL = 48000; - static constexpr size_t RUNS = 15; + static constexpr size_t RUNS = 15; GeometricBrownianMotion gbm{MU, SIGMA, S0}; EulerMaruyama euler; - struct Row { std::string name; double mse; }; + struct Row { + std::string name; + double mse; + }; auto run_comparison = [&](const char* opt_name, std::shared_ptr payoff) { auto plain_ref = make_plain_mc_sampler(gbm, euler, payoff, STEPS, T); auto ref_s = plain_ref(300'000, detail::make_seed(0, 77)); - double ref = std::accumulate(ref_s.begin(), ref_s.end(), 0.0) / - static_cast(ref_s.size()); + double ref = + std::accumulate(ref_s.begin(), ref_s.end(), 0.0) / static_cast(ref_s.size()); auto strategies = build_all_strategies(gbm, euler, payoff, CONTROL_MEAN, STEPS, T); @@ -182,8 +184,8 @@ TEST_CASE("Engine: empirical MSE vs fixed strategies", "[bench]") { double mse = 0.0; for (size_t r = 0; r < RUNS; ++r) { auto s = fn(TOTAL, detail::make_seed(r + 1000, 5)); - double est = std::accumulate(s.begin(), s.end(), 0.0) / - static_cast(s.size()); + double est = + std::accumulate(s.begin(), s.end(), 0.0) / static_cast(s.size()); double e = est - ref; mse += e * e; } @@ -195,9 +197,9 @@ TEST_CASE("Engine: empirical MSE vs fixed strategies", "[bench]") { Engine> engine{ gbm, euler, CONTROL_MEAN, STEPS, T}; IterativeEngineConfig cfg; - cfg.n_compete = 4; - cfg.n_exploit = 8; - cfg.n_rounds = 5; + cfg.n_compete = 4; + cfg.n_exploit = 8; + cfg.n_rounds = 5; cfg.samples_per_thread = 800; { double mse = 0.0; @@ -212,20 +214,24 @@ TEST_CASE("Engine: empirical MSE vs fixed strategies", "[bench]") { std::sort(rows.begin(), rows.end(), [](auto& a, auto& b) { return a.mse < b.mse; }); double plain_mse = 0.0; - for (auto& r : rows) if (r.name == "plain_mc") plain_mse = r.mse; + for (auto& r : rows) + if (r.name == "plain_mc") + plain_mse = r.mse; WARN(opt_name << " (ref=" << ref << ", budget=" << TOTAL << ", runs=" << RUNS << ")"); for (auto& r : rows) - WARN(" " << r.name << ": MSE=" << r.mse - << " (" << (plain_mse / r.mse) << "x vs plain_mc)"); + WARN(" " << r.name << ": MSE=" << r.mse << " (" << (plain_mse / r.mse) + << "x vs plain_mc)"); double engine_mse = 0.0; - for (auto& r : rows) if (r.name == "[Engine]") engine_mse = r.mse; + for (auto& r : rows) + if (r.name == "[Engine]") + engine_mse = r.mse; REQUIRE(engine_mse < plain_mse); }; - run_comparison("European Call", std::make_shared(K)); - run_comparison("Asian Call", std::make_shared(K)); + run_comparison("European Call", std::make_shared(K)); + run_comparison("Asian Call", std::make_shared(K)); run_comparison("Barrier (Up-Out B=130)", std::make_shared(K, 130.0)); - run_comparison("Lookback Call", std::make_shared()); + run_comparison("Lookback Call", std::make_shared()); } diff --git a/tests/test_strategies.cpp b/tests/test_strategies.cpp index 9691602..abd2c68 100644 --- a/tests/test_strategies.cpp +++ b/tests/test_strategies.cpp @@ -14,13 +14,13 @@ using namespace lfmc; // ─── Shared market parameters ───────────────────────────────────────────────── // ATM European-style: S0=100, K=100, mu=r=0.05, sigma=0.2, T=1y, 52 steps -static constexpr double S0 = 100.0; -static constexpr double MU = 0.05; -static constexpr double SIGMA = 0.2; -static constexpr double K = 100.0; -static constexpr double BARRIER = 130.0; // up-and-out / down-and-in barrier -static constexpr double T = 1.0; -static constexpr size_t STEPS = 52; +static constexpr double S0 = 100.0; +static constexpr double MU = 0.05; +static constexpr double SIGMA = 0.2; +static constexpr double K = 100.0; +static constexpr double BARRIER = 130.0; // up-and-out / down-and-in barrier +static constexpr double T = 1.0; +static constexpr size_t STEPS = 52; static const double CONTROL_MEAN = S0 * std::exp(MU * T); // E[S_T] = 105.13 @@ -31,10 +31,14 @@ static EulerMaruyama EULER; static std::pair stats(const std::vector& v) { const double n = static_cast(v.size()); double sum = 0.0; - for (double x : v) sum += x; + for (double x : v) + sum += x; const double mean = sum / n; double sq = 0.0; - for (double x : v) { double d = x - mean; sq += d * d; } + for (double x : v) { + double d = x - mean; + sq += d * d; + } return {mean, sq / (n - 1.0)}; } @@ -58,7 +62,8 @@ static double plain_mc_ref(std::shared_ptr payoff, size_t n = 200'000) { auto sampler = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); auto s = sampler(n, detail::make_seed(999, 999)); double sum = 0.0; - for (double x : s) sum += x; + for (double x : s) + sum += x; return sum / static_cast(s.size()); } @@ -71,7 +76,8 @@ TEST_CASE("strategies: all 10 strategies are unbiased for European call") { auto s = sampler(30'000, detail::make_seed(42, 1)); auto [mean, var] = stats(s); const double se = std::sqrt(var / s.size()); - INFO("Strategy: " << name << " mean=" << mean << " ref=" << ref << " 5-sigma=" << 5.0 * se); + INFO("Strategy: " << name << " mean=" << mean << " ref=" << ref + << " 5-sigma=" << 5.0 * se); // 5-sigma bound: P(fail) < 0.00006% REQUIRE(std::abs(mean - ref) < 5.0 * se); } @@ -101,7 +107,8 @@ TEST_CASE("strategies: importance sampling is unbiased across theta values") { static double relative_variance(const SamplerFn& strategy, const SamplerFn& plain_mc, size_t n) { auto [mean_s, var_s] = stats(strategy(n, detail::make_seed(1, 0))); auto [mean_p, var_p] = stats(plain_mc(n, detail::make_seed(2, 0))); - (void)mean_s; (void)mean_p; + (void)mean_s; + (void)mean_p; return var_s / var_p; } @@ -115,7 +122,7 @@ static double relative_variance(const SamplerFn& strategy, const SamplerFn& plai TEST_CASE("strategies: stratified sampling has lower variance than plain MC") { auto payoff = std::make_shared(K); auto stratified = make_stratified_sampler(GBM, EULER, payoff, STEPS, T); - auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); double rv = relative_variance(stratified, plain, 80'000); INFO("Variance ratio (stratified/plain): " << rv); // ~2% theoretical VR for the first of 52 dimensions; bound to 5% noise tolerance @@ -125,7 +132,7 @@ TEST_CASE("strategies: stratified sampling has lower variance than plain MC") { TEST_CASE("strategies: Halton QMC has lower variance than plain MC") { auto payoff = std::make_shared(K); auto halton = make_halton_sampler(GBM, EULER, payoff, STEPS, T); - auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); double rv = relative_variance(halton, plain, 80'000); INFO("Variance ratio (halton/plain): " << rv); REQUIRE(rv < 1.0); @@ -133,7 +140,7 @@ TEST_CASE("strategies: Halton QMC has lower variance than plain MC") { TEST_CASE("strategies: moment matching has lower variance than plain MC") { auto payoff = std::make_shared(K); - auto mm = make_moment_matching_sampler(GBM, EULER, payoff, STEPS, T); + auto mm = make_moment_matching_sampler(GBM, EULER, payoff, STEPS, T); auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); double rv = relative_variance(mm, plain, 80'000); INFO("Variance ratio (moment_matching/plain): " << rv); @@ -143,7 +150,7 @@ TEST_CASE("strategies: moment matching has lower variance than plain MC") { TEST_CASE("strategies: LHS has lower variance than plain MC") { auto payoff = std::make_shared(K); - auto lhs = make_lhs_sampler(GBM, EULER, payoff, STEPS, T); + auto lhs = make_lhs_sampler(GBM, EULER, payoff, STEPS, T); auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); double rv = relative_variance(lhs, plain, 80'000); INFO("Variance ratio (lhs/plain): " << rv); @@ -152,9 +159,9 @@ TEST_CASE("strategies: LHS has lower variance than plain MC") { } TEST_CASE("strategies: stratified antithetic has lower variance than plain MC") { - auto payoff = std::make_shared(K); - auto sa = make_stratified_antithetic_sampler(GBM, EULER, payoff, STEPS, T); - auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); + auto payoff = std::make_shared(K); + auto sa = make_stratified_antithetic_sampler(GBM, EULER, payoff, STEPS, T); + auto plain = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); double rv = relative_variance(sa, plain, 80'000); INFO("Variance ratio (stratified_antithetic/plain): " << rv); REQUIRE(rv < 1.0); @@ -172,15 +179,15 @@ static ASVRResult run_asvr_for(std::shared_ptr payoff, size_t n = 20'000 static void log_result(const ASVRResult& r) { for (size_t k = 0; k < r.exploration_stats.size(); ++k) { - INFO(" " << r.exploration_stats[k].name - << ": var=" << r.exploration_stats[k].sample_variance - << " weight=" << r.precision_weights[k] - << " n_exploit=" << r.exploitation_counts[k]); + INFO(" " << r.exploration_stats[k].name << ": var=" + << r.exploration_stats[k].sample_variance << " weight=" << r.precision_weights[k] + << " n_exploit=" << r.exploitation_counts[k]); } INFO("VR ratio: " << r.variance_reduction_ratio); - INFO("Best strategy: " << r.exploration_stats[ - std::max_element(r.precision_weights.begin(), r.precision_weights.end()) - - r.precision_weights.begin()].name); + INFO("Best strategy: " << r.exploration_stats[std::max_element(r.precision_weights.begin(), + r.precision_weights.end()) - + r.precision_weights.begin()] + .name); } TEST_CASE("ASVR with 10 strategies: European call VR ratio > 1") { @@ -240,14 +247,14 @@ TEST_CASE("ASVR with 10 strategies: LookbackPut VR ratio > 1") { TEST_CASE("ASVR: best strategy for European call vs lookback call", "[.slow]") { // Run with large budget for stable weight estimates auto r_eur = run_asvr_for(std::make_shared(K), 100'000); - auto r_lb = run_asvr_for(std::make_shared(), 100'000); + auto r_lb = run_asvr_for(std::make_shared(), 100'000); const size_t best_eur = static_cast( - std::max_element(r_eur.precision_weights.begin(), r_eur.precision_weights.end()) - - r_eur.precision_weights.begin()); + std::max_element(r_eur.precision_weights.begin(), r_eur.precision_weights.end()) - + r_eur.precision_weights.begin()); const size_t best_lb = static_cast( - std::max_element(r_lb.precision_weights.begin(), r_lb.precision_weights.end()) - - r_lb.precision_weights.begin()); + std::max_element(r_lb.precision_weights.begin(), r_lb.precision_weights.end()) - + r_lb.precision_weights.begin()); INFO("European call best: " << r_eur.exploration_stats[best_eur].name); INFO("Lookback call best: " << r_lb.exploration_stats[best_lb].name); @@ -255,15 +262,15 @@ TEST_CASE("ASVR: best strategy for European call vs lookback call", "[.slow]") { } TEST_CASE("ASVR: best strategy for Asian call vs barrier call", "[.slow]") { - auto r_asian = run_asvr_for(std::make_shared(K), 100'000); + auto r_asian = run_asvr_for(std::make_shared(K), 100'000); auto r_barrier = run_asvr_for(std::make_shared(K, BARRIER), 100'000); const size_t best_asian = static_cast( - std::max_element(r_asian.precision_weights.begin(), r_asian.precision_weights.end()) - - r_asian.precision_weights.begin()); + std::max_element(r_asian.precision_weights.begin(), r_asian.precision_weights.end()) - + r_asian.precision_weights.begin()); const size_t best_barrier = static_cast( - std::max_element(r_barrier.precision_weights.begin(), r_barrier.precision_weights.end()) - - r_barrier.precision_weights.begin()); + std::max_element(r_barrier.precision_weights.begin(), r_barrier.precision_weights.end()) - + r_barrier.precision_weights.begin()); INFO("Asian call best: " << r_asian.exploration_stats[best_asian].name); INFO("Barrier call best: " << r_barrier.exploration_stats[best_barrier].name); @@ -275,7 +282,7 @@ TEST_CASE("ASVR: best strategy for Asian call vs barrier call", "[.slow]") { // the same total sample budget. This is the core empirical contribution. TEST_CASE("strategies: empirical MSE comparison across all option types", "[.slow]") { - const size_t N = 10'000; + const size_t N = 10'000; const size_t RUNS = 100; struct OptionConfig { @@ -284,32 +291,36 @@ TEST_CASE("strategies: empirical MSE comparison across all option types", "[.slo }; std::vector options = { - {"European Call", std::make_shared(K)}, - {"European Put", std::make_shared(K)}, - {"Asian Call", std::make_shared(K)}, - {"Asian Put", std::make_shared(K)}, + {"European Call", std::make_shared(K)}, + {"European Put", std::make_shared(K)}, + {"Asian Call", std::make_shared(K)}, + {"Asian Put", std::make_shared(K)}, {"Up-And-Out Call (barrier=130)", std::make_shared(K, BARRIER)}, - {"Down-And-In Put (barrier=70)", std::make_shared(K, 70.0)}, - {"Lookback Call", std::make_shared()}, - {"Lookback Put", std::make_shared()}, + {"Down-And-In Put (barrier=70)", std::make_shared(K, 70.0)}, + {"Lookback Call", std::make_shared()}, + {"Lookback Put", std::make_shared()}, }; - INFO("\n| Option Type | Plain MC MSE | ASVR MSE | VR Ratio | Best Strategy |"); - INFO( "|---------------------------------|-------------|-----------|----------|---------------------|"); + INFO("\n| Option Type | Plain MC MSE | ASVR MSE | VR Ratio | Best " + "Strategy |"); + INFO("|---------------------------------|-------------|-----------|----------|-----------------" + "----|"); - auto plain_mc = make_plain_mc_sampler(GBM, EULER, nullptr, STEPS, T); // will be rebuilt per option + auto plain_mc = + make_plain_mc_sampler(GBM, EULER, nullptr, STEPS, T); // will be rebuilt per option for (auto& [opt_name, payoff] : options) { // Reference value from 300k plain MC samples auto plain_ref_fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); auto ref_samples = plain_ref_fn(300'000, detail::make_seed(0, 77)); double ref_mean = 0.0; - for (double x : ref_samples) ref_mean += x; + for (double x : ref_samples) + ref_mean += x; ref_mean /= static_cast(ref_samples.size()); auto plain_fn = make_plain_mc_sampler(GBM, EULER, payoff, STEPS, T); - double mse_asvr = 0.0; + double mse_asvr = 0.0; double mse_plain = 0.0; std::string best_name; @@ -320,36 +331,34 @@ TEST_CASE("strategies: empirical MSE comparison across all option types", "[.slo for (size_t r = 0; r < RUNS; ++r) { // ASVR auto strategies = build_all_strategies(GBM, EULER, payoff, CONTROL_MEAN, STEPS, T); - auto asvr_res = AdaptiveVarianceReduction::run(std::move(strategies), N, cfg); - double err_a = asvr_res.estimate - ref_mean; - mse_asvr += err_a * err_a; + auto asvr_res = AdaptiveVarianceReduction::run(std::move(strategies), N, cfg); + double err_a = asvr_res.estimate - ref_mean; + mse_asvr += err_a * err_a; if (r == 0) { - const size_t best_k = static_cast( - std::max_element(asvr_res.precision_weights.begin(), - asvr_res.precision_weights.end()) - - asvr_res.precision_weights.begin()); + const size_t best_k = + static_cast(std::max_element(asvr_res.precision_weights.begin(), + asvr_res.precision_weights.end()) - + asvr_res.precision_weights.begin()); best_name = asvr_res.exploration_stats[best_k].name; } // Plain MC auto plain_samples = plain_fn(N, detail::make_seed(r + 5000, 88)); double pm = 0.0; - for (double x : plain_samples) pm += x; + for (double x : plain_samples) + pm += x; pm /= static_cast(plain_samples.size()); double err_p = pm - ref_mean; - mse_plain += err_p * err_p; + mse_plain += err_p * err_p; } - mse_asvr /= static_cast(RUNS); + mse_asvr /= static_cast(RUNS); mse_plain /= static_cast(RUNS); const double vr = mse_plain / mse_asvr; - INFO("| " << opt_name - << " | " << mse_plain - << " | " << mse_asvr - << " | " << vr - << " | " << best_name << " |"); + INFO("| " << opt_name << " | " << mse_plain << " | " << mse_asvr << " | " << vr << " | " + << best_name << " |"); REQUIRE(mse_asvr < mse_plain); }