feat: finished refactor and added control variates

Author: ChronoInferna

include/lfmc/estimator.hpp

@@ -4,12 +4,13 @@
 
 #include <expected>
 #include <string>
+#include <vector>
 
 namespace lfmc {
 
 class Estimator {
   public:
-    virtual std::expected<void, std::string> add_payoffs(const Payoffs& payoffs) = 0;
+    virtual std::expected<void, std::string> add_payoffs(const std::vector<Payoffs>& payoffs) = 0;
     virtual bool converged() const = 0;
     virtual std::expected<double, std::string> 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<void, std::string> add_payoffs(const Payoffs& payoffs) override;
+    std::expected<void, std::string> add_payoffs(const std::vector<Payoffs>& payoffs) override;
+    bool converged() const override;
+    std::expected<double, std::string> 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<void, std::string> add_payoffs(const std::vector<Payoffs>& payoffs) override;
     bool converged() const override;
     std::expected<double, std::string> result() const override;
     // void merge(Estimator const& other) override;

include/lfmc/numerical_scheme.hpp

@@ -1,5 +1,8 @@
 #pragma once
 
+#include "lfmc/stochastic_process.hpp"
+
+#include <cmath>
 #include <concepts>
 
 namespace lfmc {
@@ -10,6 +13,15 @@ concept NumericalScheme =
         { s.step(p, t, x, dt, z) } -> std::convertible_to<double>;
     };
 
+template <StochasticProcess P> 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.
  *

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;

include/lfmc/payoff.hpp

@@ -2,11 +2,16 @@
 
 #include "lfmc/types.hpp"
 
+#include <expected>
+#include <memory>
+#include <vector>
+
 namespace lfmc {
 
 class Payoff {
   public:
-    virtual Payoffs generate_payoffs(const std::vector<Path>& paths) const = 0;
+    virtual std::expected<std::vector<Payoffs>, std::string>
+    generate_payoffs(const std::vector<Path>& 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<Path>& paths) const override;
+    std::expected<std::vector<Payoffs>, std::string>
+    generate_payoffs(const std::vector<Path>& 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<Path>& paths) const override;
+    std::expected<std::vector<Payoffs>, std::string>
+    generate_payoffs(const std::vector<Path>& 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<Payoff> target_payoff_;
+    std::unique_ptr<Payoff> control_payoff_;
+
+  public:
+    explicit ControlVariatePayoff(std::unique_ptr<Payoff> target_payoff,
+                                  std::unique_ptr<Payoff> control_payoff);
+
+    std::expected<std::vector<Payoffs>, std::string>
+    generate_payoffs(const std::vector<Path>& paths) const override;
 };
 
 } // namespace lfmc

include/lfmc/pipeline.hpp

@@ -26,11 +26,25 @@ template <StochasticProcess SP, NumericalScheme<SP> NS> class Pipeline {
 
     std::expected<double, std::string> 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());
             }

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<double> dist_;
+    bool toggle_ = false;
+
+  public:
+    AntitheticRandomSource(unsigned seed = std::random_device{}());
+
+    std::vector<Normals> generate_normals(size_t steps, size_t) override;
+
+    void seed(unsigned seed);
+};
+
 } // namespace lfmc

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

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

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<Normals> AntitheticRandomSource::generate_normals(size_t steps, size_t samples) {
+    std::vector<Normals> 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

src/control_variate_estimator.cpp

@@ -0,0 +1,53 @@
+#include "lfmc/estimator.hpp"
+
+namespace lfmc {
+
+std::expected<void, std::string>
+ControlVariateEstimator::add_payoffs(const std::vector<Payoffs>& 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<double, std::string> 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<double>(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