Restructure C++ core into cpp module and package bindings.

Move the pricing engine sources out of src/ into cpp/, add the closed-form engine and pybind wiring, and align tests/build targets with the new project layout.

Made-with: Cursor

cpp/BSWrapper.cpp

@@ -0,0 +1,49 @@
+//
+// Created by David Doebel  on 27.03.2026.
+//
+
+#include "BSWrapper.hpp"
+
+#include "BlackScholesClosedFormEngine.hpp"
+#include "BlackScholesProcess.hpp"
+#include "Instrument.hpp"
+#include "Option.hpp"
+#include "FlatVolatilitySurface.hpp"
+#include "FlatYieldCurve.hpp"
+#include <cassert>
+#include <iostream>
+
+class FlatYieldCurve;
+
+double BSWrapper::bs_price_wrapper(double S, double K, double T, double r, double sigma,
+    bool is_call) {
+    std::shared_ptr<FlatYieldCurve> flat_curve = std::make_shared<FlatYieldCurve>(r);
+    auto flat_vol_surface = std::make_shared<FlatVolatilitySurface>(sigma);
+    MarketData data(S,flat_curve, flat_vol_surface);
+    std::unique_ptr<BlackScholesProcess> process = std::make_unique<BlackScholesProcess>(data);
+    std::unique_ptr<BlackScholesClosedFormEngine> pricing_engine =
+        std::make_unique<BlackScholesClosedFormEngine>(std::move(process));
+    std::unique_ptr<Payoff> payoff;
+    if (is_call)
+        payoff = std::make_unique<CallPayoff>(K);
+    else payoff = std::make_unique<PutPayoff>(K);
+    EuropeanExercise exercise(T);
+    VanillaOption option(T,std::make_unique<EuropeanExercise>(exercise),
+        std::move(payoff),std::move(pricing_engine));
+    return option.price();
+}
+
+std::vector<double> BSWrapper::batch_bs_price_wrapper(const std::vector<double> &S, const std::vector<double> &K,
+    const std::vector<double> &T, const std::vector<double> &r, const std::vector<double> &sigma,
+    const std::vector<bool> &is_call) {
+    assert(K.size() == S.size() && K.size() == T.size() && K.size() == r.size() && K.size() ==
+        sigma.size() && K.size() == is_call.size());
+    std::size_t n = K.size();
+    std::vector<double> result(n);
+    for (std::size_t i = 0; i < n; ++i) {
+        result[i] = bs_price_wrapper(S[i], K[i], T[i], r[i], sigma[i], is_call[i]);
+        if (i % 100 == 0)
+            std::cout << "i = " << i << " result = " << result[i] << std::endl; // ( i % 1000 == 0)
+    }
+    return result;
+}

cpp/BSWrapper.hpp

@@ -0,0 +1,24 @@
+/**
+ * @file BSWrapper.hpp
+ * @brief Black–Scholes vanilla price (closed form; used from Python / implied vol).
+ */
+
+#ifndef QUANTENGINE_BSWRAPPER_HPP
+#define QUANTENGINE_BSWRAPPER_HPP
+#include <vector>
+
+/**
+ * @brief Static helpers wrapping scalar and batch pricing.
+ */
+class BSWrapper {
+public:
+    BSWrapper() = delete;
+    static double bs_price_wrapper(double S, double K, double T, double r, double sigma, bool is_call);
+    static std::vector<double> batch_bs_price_wrapper(const std::vector<double>& S, const std::vector<double>& K,
+        const std::vector<double>& T, const std::vector<double>& r, const std::vector<double>& sigma,
+        const std::vector<bool>& is_call);
+
+};
+
+
+#endif //QUANTENGINE_BSWRAPPER_HPP

cpp/BlackScholesClosedFormEngine.cpp

@@ -0,0 +1,69 @@
+/**
+ * @file BlackScholesClosedFormEngine.cpp
+ * @brief Black–Scholes closed-form pricing (calls, puts, cash-or-nothing digital).
+ */
+
+#include "BlackScholesClosedFormEngine.hpp"
+#include "Instrument.hpp"
+#include "Payoff.hpp"
+#include <cmath>
+#include <stdexcept>
+
+namespace {
+
+double norm_cdf(double x) {
+    return 0.5 * (1.0 + std::erf(x / std::sqrt(2.0)));
+}
+
+} // namespace
+
+double BlackScholesClosedFormEngine::calculate(const Instrument &instrument) const {
+    if (instrument.exerciseType() != Exercise::Type::European) {
+        throw std::invalid_argument("BlackScholesClosedFormEngine: European exercise only");
+    }
+
+    const double T = instrument.maturity();
+    const MarketData &md = process_->data();
+    const double S = md.spot();
+    double K = instrument.payoff().strike();
+    const PayoffKind pk = instrument.payoff().kind();
+
+    if (T <= 0.0) {
+        return instrument.payoff()(S);
+    }
+
+    const double r = md.yield_curve().zeroRate(T);
+    const double sigma = md.volatility_surface().sigma(K, T);
+    if (sigma <= 0.0) {
+        throw std::invalid_argument("BlackScholesClosedFormEngine: volatility must be positive");
+    }
+
+    const double disc = md.yield_curve().discount(T);
+    const double sqrtT = std::sqrt(T);
+    const double sig_sqrtT = sigma * sqrtT;
+
+    if (sig_sqrtT < 1e-15) {
+        const double forward = S * std::exp(r * T);
+        switch (pk) {
+        case PayoffKind::Call:
+            return disc * std::max(0.0, forward - K);
+        case PayoffKind::Put:
+            return disc * std::max(0.0, K - forward);
+        case PayoffKind::Digital:
+            return (forward > K) ? disc : 0.0;
+        }
+    }
+
+    const double d1 = (std::log(S / K) + (r + 0.5 * sigma * sigma) * T) / sig_sqrtT;
+    const double d2 = d1 - sig_sqrtT;
+
+    switch (pk) {
+    case PayoffKind::Call:
+        return S * norm_cdf(d1) - K * disc * norm_cdf(d2);
+    case PayoffKind::Put:
+        return K * disc * norm_cdf(-d2) - S * norm_cdf(-d1);
+    case PayoffKind::Digital:
+        return disc * norm_cdf(d2);
+    }
+    throw std::logic_error("BlackScholesClosedFormEngine: unhandled PayoffKind");
+}

cpp/BlackScholesClosedFormEngine.hpp

@@ -0,0 +1,22 @@
+/**
+ * @file BlackScholesClosedFormEngine.hpp
+ * @brief Risk-neutral Black–Scholes formula for European payoffs under GBM (flat or surface inputs via @ref MarketData).
+ */
+
+#ifndef QUANTENGINE_BLACKSCHOLESCLOSEDFORMENGINE_HPP
+#define QUANTENGINE_BLACKSCHOLESCLOSEDFORMENGINE_HPP
+
+#include "PricingEngine.hpp"
+
+/**
+ * @brief Analytic European vanilla / digital prices using @f$r@f$ and @f$\sigma(K,T)@f$ from the embedded process’s @ref MarketData.
+ */
+class BlackScholesClosedFormEngine : public PricingEngine {
+public:
+    explicit BlackScholesClosedFormEngine(std::unique_ptr<StochasticProcess> process)
+        : PricingEngine(std::move(process)) {}
+
+    double calculate(const Instrument &instrument) const override;
+};
+
+#endif // QUANTENGINE_BLACKSCHOLESCLOSEDFORMENGINE_HPP

cpp/BlackScholesProcess.cpp

@@ -0,0 +1,24 @@
+/**
+ * @file BlackScholesProcess.cpp
+ * @brief Black–Scholes GBM drift, diffusion, and step.
+ */
+
+#include "BlackScholesProcess.hpp"
+
+double BlackScholesProcess::drift(double t, double s) {
+    double r = this->data().yield_curve().zeroRate(t);
+    return r * s;
+}
+
+double BlackScholesProcess::diffusion(double t, double s) {
+    double sigma = this->data().volatility_surface().sigma(s,t);
+    return sigma*s;
+}
+
+double BlackScholesProcess::step(double t, double s, double dt, double dW) {
+    double r = this->data().yield_curve().zeroRate(t);
+    double sigma = this->data().volatility_surface().sigma(s,t);
+    return s*exp((r-0.5*sigma*sigma)*dt + sigma*sqrt(dt)*dW);
+}
+
+

cpp/BlackScholesProcess.hpp

@@ -0,0 +1,26 @@
+/**
+ * @file BlackScholesProcess.hpp
+ * @brief Geometric Brownian motion with yield and volatility surfaces.
+ */
+
+#ifndef QUANTENGINE_BLACKSCHOLESPROCESS_HPP
+#define QUANTENGINE_BLACKSCHOLESPROCESS_HPP
+#include "StochasticProcess.hpp"
+
+/**
+ * @brief GBM: drift @f$r_t S@f$, diffusion @f$\sigma(S,t) S@f$, exact log-step.
+ */
+class BlackScholesProcess : public StochasticProcess{
+public:
+    explicit BlackScholesProcess(MarketData data) : StochasticProcess(std::move(data)){}
+
+    double drift(double t, double s) override;
+
+    double diffusion(double t, double s) override;
+
+    double step(double t, double s, double dt, double dW) override;
+
+};
+
+
+#endif //QUANTENGINE_BLACKSCHOLESPROCESS_HPP

cpp/CMakeLists.txt

@@ -0,0 +1,65 @@
+add_library(qengine_core
+        Instrument.cpp
+        Instrument.hpp
+        Payoff.cpp
+        Payoff.hpp
+        Option.cpp
+        Option.hpp
+        PricingEngine.cpp
+        PricingEngine.hpp
+        MonteCarloEngine.cpp
+        MonteCarloEngine.hpp
+        StochasticProcess.cpp
+        StochasticProcess.hpp
+        Exercise.cpp
+        Exercise.hpp
+        MarketData.cpp
+        MarketData.hpp
+        YieldCurve.cpp
+        YieldCurve.hpp
+        VolatilitySurface.cpp
+        VolatilitySurface.hpp
+        RandomGenerator.cpp
+        RandomGenerator.hpp
+        Statistics.cpp
+        Statistics.hpp
+        BlackScholesClosedFormEngine.cpp
+        BlackScholesClosedFormEngine.hpp
+        BlackScholesProcess.cpp
+        BlackScholesProcess.hpp
+        DBIngest.cpp
+        DBIngest.hpp
+        BSWrapper.cpp
+        BSWrapper.hpp
+        NewtonSolver.cpp
+        NewtonSolver.hpp
+)
+
+target_include_directories(qengine_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+target_include_directories(qengine_core PRIVATE
+        /opt/homebrew/include
+)
+
+find_library(PQXX_LIB pqxx PATHS /opt/homebrew/lib /usr/local/lib /usr/lib)
+find_library(PQ_LIB pq PATHS /opt/homebrew/opt/libpq/lib /opt/homebrew/lib /usr/local/lib /usr/lib)
+if(NOT PQXX_LIB OR NOT PQ_LIB)
+    message(FATAL_ERROR "Could not find libpqxx and/or libpq (install via Homebrew: brew install libpqxx libpq)")
+endif()
+
+target_link_libraries(qengine_core Eigen3::Eigen)
+target_link_libraries(qengine_core ${PQXX_LIB} ${PQ_LIB})
+
+# Python import path: package qengine, extension submodule qengine (file qengine/qengine*.so)
+pybind11_add_module(qengine_cpp MODULE ImpliedVolatility/Pybind.cpp)
+set_target_properties(qengine_cpp PROPERTIES OUTPUT_NAME qengine)
+target_link_libraries(qengine_cpp PRIVATE qengine_core)
+
+# Place the module next to qengine/__init__.py so `import qengine` works from the repo root
+set(_qengine_py_pkg "${CMAKE_SOURCE_DIR}/qengine")
+set_target_properties(qengine_cpp PROPERTIES
+        LIBRARY_OUTPUT_DIRECTORY "${_qengine_py_pkg}"
+        LIBRARY_OUTPUT_DIRECTORY_RELEASE "${_qengine_py_pkg}"
+        LIBRARY_OUTPUT_DIRECTORY_DEBUG "${_qengine_py_pkg}"
+        RUNTIME_OUTPUT_DIRECTORY "${_qengine_py_pkg}"
+        RUNTIME_OUTPUT_DIRECTORY_RELEASE "${_qengine_py_pkg}"
+        RUNTIME_OUTPUT_DIRECTORY_DEBUG "${_qengine_py_pkg}")

cpp/DBIngest.cpp

@@ -0,0 +1,64 @@
+/**
+ * @file DBIngest.cpp
+ * @brief Database connection and placeholder update routines.
+ */
+
+#include "DBIngest.hpp"
+
+#include <cstdlib>
+#include <iostream>
+#include <sstream>
+
+bool DBIngest::connect() {
+    const char* db_name = std::getenv("DB_NAME");
+    const char* db_user = std::getenv("DB_USER");
+    const char* db_password = std::getenv("DB_PASSWORD");
+    const char* db_host = std::getenv("DB_HOST");
+    const char* db_port = std::getenv("DB_PORT");
+
+    std::ostringstream conn_str;
+    conn_str
+        << "dbname=" << (db_name ? db_name : "options_db")
+        << " user=" << (db_user ? db_user : "quant_user")
+        << " host=" << (db_host ? db_host : "localhost")
+        << " port=" << (db_port ? db_port : "5432")
+        << " password=" << (db_password ? db_password : "");
+
+    connection_ = pqxx::connection(conn_str.str());
+
+    if(connection_.is_open()) {
+        std::cout << "Connected\n";
+        return true;
+    }
+    std::cout << "Not connected\n";
+    return false;
+}
+
+bool DBIngest::disconnect() {
+    connection_.close();
+    return true;
+}
+
+bool DBIngest::update(VolatilitySurface &surface) {
+    std::string vol_surface_query = "SELECT c.strike, c.expiration_date, q.mid, u.price "
+                                    "FROM option_quotes q"
+                                    "JOIN option_contracts c "
+                                    "ON q.contract_id = c.id "
+                                    "JOIN underlying_prices u"
+                                    "ON u.underlying_id = c.underlying_id"
+                                    "WHERE q.timestamp = ("
+                                    "SELECT MAX(timestamp) FROM option_quotes"
+                                    ")";
+    pqxx::work work(connection_);
+    pqxx::result result = work.exec(vol_surface_query);
+    for (auto row : result) {
+        std::cout << row[0] << " " << row[1] << " " << row[2] << " " << row[3] << std::endl;
+    }
+    (void)surface;
+    return false;
+}
+
+bool DBIngest::update(YieldCurve &yield_curve) {
+    (void)yield_curve;
+    return false;
+}

cpp/DBIngest.hpp

@@ -0,0 +1,28 @@
+/**
+ * @file DBIngest.hpp
+ * @brief PostgreSQL helpers to load market objects (work in progress).
+ */
+
+#ifndef QUANTENGINE_DBINGEST_HPP
+#define QUANTENGINE_DBINGEST_HPP
+
+#include <pqxx/pqxx>
+
+#include "VolatilitySurface.hpp"
+#include "YieldCurve.hpp"
+
+/**
+ * @brief Connects to Postgres via libpqxx and queries quotes for surface building.
+ */
+class DBIngest {
+
+    bool connect();
+    bool disconnect();
+    bool update(VolatilitySurface& surface);
+    bool update(YieldCurve& yield_curve);
+private:
+    pqxx::connection connection_;
+};
+
+
+#endif //QUANTENGINE_DBINGEST_HPP

cpp/Exercise.cpp

@@ -0,0 +1,6 @@
+/**
+ * @file Exercise.cpp
+ * @brief @ref Exercise translation unit (interface only).
+ */
+
+#include "Exercise.hpp"

cpp/Exercise.hpp

@@ -0,0 +1,61 @@
+/**
+ * @file Exercise.hpp
+ * @brief Exercise style (European, American, Bermudan) and exercise times.
+ */
+
+#ifndef QUANTENGINE_EXERCISE_HPP
+#define QUANTENGINE_EXERCISE_HPP
+#include <vector>
+
+/**
+ * @brief Describes when the holder may exercise (metadata for pricing engines).
+ */
+class Exercise {
+public:
+    Exercise() = default;
+    virtual ~Exercise() = default;
+    enum class Type {
+        European,
+        American,
+        Bermudan
+    };
+
+    virtual Type type() const = 0;
+protected:
+    std::vector<double> exercise_times_;
+
+};
+
+/** @brief Single exercise at maturity. */
+class EuropeanExercise : public Exercise {
+public:
+    EuropeanExercise() : type_(Type::European) {};
+    EuropeanExercise(double maturity) : type_(Type::European){
+        exercise_times_.push_back(maturity);
+    }
+    ~EuropeanExercise() override = default;
+    [[nodiscard]] Type type() const override {
+        return type_;
+    }
+private:
+    Type type_;
+};
+
+/** @brief Continuous American exercise from @f$t=0@f$ to maturity (placeholder grid). */
+class AmericanExercise : public Exercise{
+public:
+    AmericanExercise() : type_(Type::American) {};
+    AmericanExercise(double maturity) : type_(Type::American) {
+        exercise_times_.push_back(0);
+        exercise_times_.push_back(maturity);
+    }
+    [[nodiscard]] Type type() const override {
+        return type_;
+    }
+
+private:
+    Type type_;
+};
+
+
+#endif //QUANTENGINE_EXERCISE_HPP

cpp/FlatVolatilitySurface.cpp

@@ -0,0 +1,5 @@
+/**
+ * @file FlatVolatilitySurface.cpp
+ * @brief Ensures link visibility for @ref FlatVolatilitySurface.
+ */
+#include "FlatVolatilitySurface.hpp"

cpp/FlatVolatilitySurface.hpp

@@ -0,0 +1,21 @@
+/**
+ * @file FlatVolatilitySurface.hpp
+ * @brief Constant implied volatility surface.
+ */
+#ifndef QUANTENGINE_FLATVOLATILITYSURFACE_HPP
+#define QUANTENGINE_FLATVOLATILITYSURFACE_HPP
+#include "VolatilitySurface.hpp"
+
+/**
+ * @brief @f$\sigma(K,T)\equiv\sigma_0@f$.
+ */
+class FlatVolatilitySurface : public VolatilitySurface {
+public:
+    explicit FlatVolatilitySurface(double sigma = 0.2) : sigma_(sigma) {}
+
+    double sigma(double K, double T) const override {return sigma_;}
+
+private:
+    double sigma_;
+};
+#endif

cpp/FlatYieldCurve.cpp

@@ -0,0 +1,5 @@
+/**
+ * @file FlatYieldCurve.cpp
+ * @brief Ensures link visibility for @ref FlatYieldCurve (inline methods in header).
+ */
+#include "FlatYieldCurve.hpp"

cpp/FlatYieldCurve.hpp

@@ -0,0 +1,22 @@
+/**
+ * @file FlatYieldCurve.hpp
+ * @brief Constant zero rate yield curve.
+ */
+#ifndef QUANTENGINE_FLATYIELDCURVE_HPP
+#define QUANTENGINE_FLATYIELDCURVE_HPP
+#include "YieldCurve.hpp"
+#include <cmath>
+
+/**
+ * @brief @f$P(t)=e^{-r t}@f$, @f$f(t)\equiv r@f$.
+ */
+class FlatYieldCurve : public YieldCurve{
+public:
+    explicit FlatYieldCurve(double rate = 0.01) : rate_(rate) {}
+
+    double discount(double t) const override {return std::exp(-rate_ * t); };
+    double zeroRate(double t) const override {return rate_; }
+private:
+    double rate_ = 0.01;
+};
+#endif

cpp/ImpliedVolatility/Pybind.cpp

@@ -0,0 +1,93 @@
+/**
+ * @file Pybind.cpp
+ * @brief pybind11 module @c qengine exposing @ref BSWrapper::bs_price_wrapper overloads.
+ */
+
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <cstdint>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include "BSWrapper.hpp"
+
+namespace py = pybind11;
+
+namespace {
+
+std::vector<double> to_vector_double(const py::array_t<double> &a) {
+    py::buffer_info info = a.request();
+    if (info.ndim != 1) {
+        throw std::runtime_error("expected 1-D ndarray for S, K, T, r, sigma");
+    }
+    const auto *p = static_cast<const double *>(info.ptr);
+    const ssize_t n = info.shape[0];
+    return std::vector<double>(p, p + n);
+}
+
+std::vector<bool> to_vector_bool_1d(const py::array_t<bool> &a) {
+    py::buffer_info info = a.request();
+    if (info.ndim != 1) {
+        throw std::runtime_error("expected 1-D ndarray for is_call");
+    }
+    if (info.itemsize != 1) {
+        throw std::runtime_error("is_call: expected a boolean ndarray (dtype=bool)");
+    }
+    const ssize_t n = info.shape[0];
+    const auto *p = static_cast<const std::uint8_t *>(info.ptr);
+    std::vector<bool> out(static_cast<size_t>(n));
+    for (ssize_t i = 0; i < n; ++i) {
+        out[static_cast<size_t>(i)] = (p[i] != 0);
+    }
+    return out;
+}
+
+void check_same_length(size_t n, size_t k, const char *name) {
+    if (n != k) {
+        throw std::runtime_error(std::string("length mismatch for ") + name);
+    }
+}
+
+} // namespace
+
+PYBIND11_MODULE(qengine, m) {
+    m.doc() = "Binding for the Black Scholes model";
+
+    m.def(
+        "bs_price",
+        [](double S, double K, double T, double r, double sigma, bool is_call) {
+            return BSWrapper::bs_price_wrapper(S, K, T, r, sigma, is_call);
+        },
+        py::arg("S"), py::arg("K"), py::arg("T"), py::arg("r"), py::arg("sigma"), py::arg("is_call"));
+
+    m.def(
+        "bs_price",
+        [](py::array_t<double> S, py::array_t<double> K, py::array_t<double> T, py::array_t<double> r,
+            py::array_t<double> sigma, py::array_t<bool> is_call) {
+            std::vector<double> vS = to_vector_double(S);
+            std::vector<double> vK = to_vector_double(K);
+            std::vector<double> vT = to_vector_double(T);
+            std::vector<double> vr = to_vector_double(r);
+            std::vector<double> vsig = to_vector_double(sigma);
+            std::vector<bool> vC = to_vector_bool_1d(is_call);
+            const size_t n = vS.size();
+            check_same_length(n, vK.size(), "K");
+            check_same_length(n, vT.size(), "T");
+            check_same_length(n, vr.size(), "r");
+            check_same_length(n, vsig.size(), "sigma");
+            check_same_length(n, vC.size(), "is_call");
+            return BSWrapper::batch_bs_price_wrapper(vS, vK, vT, vr, vsig, vC);
+        },
+        py::arg("S"), py::arg("K"), py::arg("T"), py::arg("r"), py::arg("sigma"), py::arg("is_call"));
+
+    m.def(
+        "bs_price",
+        [](const std::vector<double> &S, const std::vector<double> &K, const std::vector<double> &T,
+            const std::vector<double> &r, const std::vector<double> &sigma, const std::vector<bool> &is_call) {
+            return BSWrapper::batch_bs_price_wrapper(S, K, T, r, sigma, is_call);
+        },
+        py::arg("S"), py::arg("K"), py::arg("T"), py::arg("r"), py::arg("sigma"), py::arg("is_call"));
+}

cpp/Instrument.cpp

@@ -0,0 +1,18 @@
+/**
+ * @file Instrument.cpp
+ * @brief @ref Instrument implementation.
+ */
+
+#include "Instrument.hpp"
+
+Instrument::Instrument(double maturity, std::unique_ptr<Payoff> payoff,
+    std::unique_ptr<PricingEngine> engine) : maturity_(maturity), payoff_(std::move(payoff)), engine_
+(std::move(engine)){
+}
+
+
+
+double Instrument::price() const {
+    return engine_->calculate(*this);
+}
+

cpp/Instrument.hpp

@@ -0,0 +1,42 @@
+/**
+ * @file Instrument.hpp
+ * @brief Generic derivative instrument: payoff plus pricing engine.
+ */
+
+#ifndef QUANTENGINE_INSTRUMENT_HPP
+#define QUANTENGINE_INSTRUMENT_HPP
+#include "Exercise.hpp"
+#include "Payoff.hpp"
+#include "PricingEngine.hpp"
+#include <memory>
+
+class PricingEngine;
+
+/**
+ * @brief Represents a tradeable claim priced via a @ref PricingEngine.
+ */
+class Instrument {
+public:
+    Instrument() = default;
+    Instrument(double maturity, std::unique_ptr<Payoff> payoff, std::unique_ptr<PricingEngine> engine);
+    double price() const;
+
+    [[nodiscard]] double maturity() const {
+        return maturity_;
+    }
+
+    [[nodiscard]] Payoff& payoff() const {
+        return *payoff_;
+    }
+
+    /** @brief Base @ref Instrument is treated as European unless overridden by @ref Option. */
+    [[nodiscard]] virtual Exercise::Type exerciseType() const { return Exercise::Type::European; }
+
+protected:
+    double maturity_;
+    std::unique_ptr<Payoff> payoff_;
+    std::unique_ptr<PricingEngine> engine_;
+};
+
+
+#endif //QUANTENGINE_INSTRUMENT_HPP

cpp/MarketData.cpp

@@ -0,0 +1,10 @@
+/**
+ * @file MarketData.cpp
+ * @brief @ref MarketData accessors.
+ */
+
+#include "MarketData.hpp"
+
+double MarketData::spot() const { return spot_; }
+const YieldCurve& MarketData::yield_curve() const { return *yield_curve_; }
+const VolatilitySurface& MarketData::volatility_surface() const { return *volatility_surface_; }

cpp/MarketData.hpp

@@ -0,0 +1,37 @@
+/**
+ * @file MarketData.hpp
+ * @brief Spot, discount curve, and volatility surface bundle.
+ */
+
+#ifndef QUANTENGINE_MARKETDATA_HPP
+#define QUANTENGINE_MARKETDATA_HPP
+#include "YieldCurve.hpp"
+#include "VolatilitySurface.hpp"
+#include <memory>
+
+/**
+ * @brief Immutable snapshot of inputs needed to simulate or price.
+ */
+class MarketData {
+public:
+    MarketData() = delete;
+
+    MarketData(double spot, std::shared_ptr<const YieldCurve> yield_curve,
+         std::shared_ptr<const VolatilitySurface> volatility_surface)
+        : spot_(spot),
+          yield_curve_(std::move(yield_curve)),
+          volatility_surface_(std::move(volatility_surface)) {
+    }
+
+    double spot() const;
+    const YieldCurve& yield_curve() const;
+    const VolatilitySurface& volatility_surface() const;
+
+private:
+    double spot_;
+    std::shared_ptr<const YieldCurve> yield_curve_;
+    std::shared_ptr<const VolatilitySurface> volatility_surface_;
+};
+
+
+#endif //QUANTENGINE_MARKETDATA_HPP

cpp/MonteCarloEngine.cpp

@@ -0,0 +1,25 @@
+/**
+ * @file MonteCarloEngine.cpp
+ * @brief Monte Carlo mean estimator with discounting.
+ */
+
+#include "MonteCarloEngine.hpp"
+#include <iostream>
+#include "Instrument.hpp"
+#include "Statistics.hpp"
+
+double MonteCarloEngine::calculate(const Instrument &instrument) const {
+    // parameters
+    double T = instrument.maturity();
+    double spot = process_->data().spot();
+    Statistics stats;
+
+    auto rNumbers = rng_->nextGaussianVector(numPaths_);
+    std::vector<double> payoffs(numPaths_);
+    for (std::size_t i = 0; i < numPaths_; ++i) {
+        double terminalPrice = process_->step(0.0,spot,T,rNumbers[i]);
+        double payoff = instrument.payoff()(terminalPrice);
+        stats.dump(payoff);
+    }
+    return stats.mean() * process_->data().yield_curve().discount(T);
+}

cpp/MonteCarloEngine.hpp

@@ -0,0 +1,26 @@
+/**
+ * @file MonteCarloEngine.hpp
+ * @brief Monte Carlo pricing using a @ref StochasticProcess and @ref RandomGenerator.
+ */
+
+#ifndef QUANTENGINE_MONTECARLOENGINE_HPP
+#define QUANTENGINE_MONTECARLOENGINE_HPP
+#include "PricingEngine.hpp"
+#include "RandomGenerator.hpp"
+
+/**
+ * @brief Simple path simulation: one Euler/exact step to horizon, average discounted payoff.
+ */
+class MonteCarloEngine : public PricingEngine{
+public:
+    MonteCarloEngine() = default;
+    MonteCarloEngine(int numPaths, std::unique_ptr<StochasticProcess> process, std::shared_ptr<RandomGenerator> rng):
+    numPaths_(numPaths), PricingEngine(std::move(process)), rng_(std::move(rng)) {}
+    double calculate(const Instrument& instrument) const override;
+private:
+    int numPaths_;
+    std::shared_ptr<RandomGenerator> rng_;
+};
+
+
+#endif //QUANTENGINE_MONTECARLOENGINE_HPP

cpp/NewtonSolver.cpp

@@ -0,0 +1,8 @@
+/**
+ * @file NewtonSolver.cpp
+ * @brief Placeholder translation unit for @ref NewtonSolver.
+ */
+
+#include "NewtonSolver.hpp"
+
+

cpp/NewtonSolver.hpp

@@ -0,0 +1,30 @@
+/**
+ * @file NewtonSolver.hpp
+ * @brief Generic Newton iteration helper (incomplete / reserved for solvers).
+ */
+
+#ifndef QUANTENGINE_GAUSSSOLVER_HPP
+#define QUANTENGINE_GAUSSSOLVER_HPP
+
+#include <functional>
+
+/**
+ * @brief Template Newton step loop with relative/absolute tolerances.
+ */
+class NewtonSolver {
+    template<typename F, typename DFinv, typename T>
+    bool solve(F&& func, DFinv&& dfinv,T x0 , double rtol, double atol) {
+        T x = x0;
+        int i = 0;
+        T increment;
+        do {
+            increment = dfinv(x) * func(x);
+            x -= increment;
+            ++i;
+        } while (i < 1000 && std::abs(increment)/ std::abs(x) > rtol && std::abs(increment) > atol);
+
+    }
+};
+
+
+#endif //QUANTENGINE_GAUSSSOLVER_HPP

cpp/Option.cpp

@@ -0,0 +1,11 @@
+/**
+ * @file Option.cpp
+ * @brief @ref Option implementation.
+ */
+
+#include "Option.hpp"
+
+Option::Option(double maturity, std::unique_ptr<Exercise> exercise, std::unique_ptr<Payoff> payoff,
+    std::unique_ptr<PricingEngine> engine) : Instrument(maturity, std::move(payoff),
+        std::move(engine)), exercise_(std::move(exercise)){
+}

cpp/Option.hpp

@@ -0,0 +1,40 @@
+/**
+ * @file Option.hpp
+ * @brief Option instrument with exercise style (@ref Exercise).
+ */
+
+#ifndef QUANTENGINE_OPTION_HPP
+#define QUANTENGINE_OPTION_HPP
+#include "Instrument.hpp"
+#include "Exercise.hpp"
+
+/**
+ * @brief Extends @ref Instrument with exercise schedule / style metadata.
+ */
+class Option : public Instrument{
+public:
+    Option() = default;
+    virtual ~Option() = default;
+    Option(double maturity, std::unique_ptr<Exercise> exercise,
+        std::unique_ptr<Payoff> payoff, std::unique_ptr<PricingEngine> engine);
+    [[nodiscard]] Exercise& exercise() const {
+        return *exercise_;
+    }
+
+    [[nodiscard]] Exercise::Type exerciseType() const override { return exercise_->type(); }
+
+protected:
+    std::unique_ptr<Exercise> exercise_;
+};
+
+/** @brief Plain-vanilla option using the base @ref Option constructor. */
+class VanillaOption : public Option {
+public:
+    using Option::Option;
+};
+
+
+
+
+
+#endif //QUANTENGINE_OPTION_HPP

cpp/Payoff.cpp

@@ -0,0 +1,19 @@
+/**
+ * @file Payoff.cpp
+ * @brief Payoff function implementations.
+ */
+
+#include "Payoff.hpp"
+#include <algorithm>
+
+double CallPayoff::operator()(double S) {
+    return std::max(0., S - strike_);
+}
+
+double PutPayoff::operator()(double S) {
+    return std::max(0., strike_ - S);
+}
+
+double DigitalPayoff::operator()(double S) {
+    return S > strike_ ? 1. : 0.;
+}

cpp/Payoff.hpp

@@ -0,0 +1,66 @@
+/**
+ * @file Payoff.hpp
+ * @brief Payoff interface and standard European payoffs (call, put, digital).
+ */
+
+#ifndef QUANTENGINE_PAYOFF_HPP
+#define QUANTENGINE_PAYOFF_HPP
+
+/**
+ * @brief Standard payoff shapes for routing (e.g. analytic vs Monte Carlo).
+ */
+enum class PayoffKind { Call, Put, Digital };
+
+/**
+ * @brief Terminal payoff as a function of spot @f$S_T@f$.
+ */
+class Payoff {
+public:
+
+
+    Payoff() = default;
+    virtual ~Payoff() = default;
+    virtual double operator()(double S) = 0;
+    virtual double strike() = 0;
+    [[nodiscard]] virtual PayoffKind kind() const = 0;
+};
+
+/** @brief Standard European call @f$\max(S-K,0)@f$. */
+class CallPayoff : public Payoff {
+public:
+    CallPayoff() = default;
+    CallPayoff(double strike) : strike_(strike) {}
+    double operator()(double S) override;
+    double strike() override {return strike_;}
+    [[nodiscard]] PayoffKind kind() const override { return PayoffKind::Call; }
+
+private:
+    double strike_;
+};
+
+/** @brief Standard European put @f$\max(K-S,0)@f$. */
+class PutPayoff : public Payoff {
+public:
+    PutPayoff() = default;
+    PutPayoff(double strike) : strike_(strike) {}
+    double operator()(double S) override;
+    double strike() override {return strike_;}
+    [[nodiscard]] PayoffKind kind() const override { return PayoffKind::Put; }
+private:
+    double strike_;
+};
+
+/** @brief Digital (cash-or-nothing style) payoff @f$1_{S>K}@f$. */
+class DigitalPayoff : public Payoff {
+    public:
+    DigitalPayoff() = default;
+    DigitalPayoff(double strike) : strike_(strike) {}
+    double operator()(double S) override;
+    double strike() override {return strike_;}
+    [[nodiscard]] PayoffKind kind() const override { return PayoffKind::Digital; }
+private:
+    double strike_;
+};
+
+
+#endif //QUANTENGINE_PAYOFF_HPP

cpp/PricingEngine.cpp

@@ -0,0 +1,6 @@
+/**
+ * @file PricingEngine.cpp
+ * @brief @ref PricingEngine translation unit (interface only).
+ */
+
+#include "PricingEngine.hpp"

cpp/PricingEngine.hpp

@@ -0,0 +1,30 @@
+/**
+ * @file PricingEngine.hpp
+ * @brief Abstract pricer for @ref Instrument given a stochastic model.
+ */
+
+#ifndef QUANTENGINE_PRICINGENGINE_HPP
+#define QUANTENGINE_PRICINGENGINE_HPP
+#include <memory>
+
+#include "StochasticProcess.hpp"
+
+class Instrument;
+
+/**
+ * @brief Computes model price of an instrument (e.g. Monte Carlo, PDE, closed form).
+ */
+class PricingEngine {
+public:
+    PricingEngine() = default;
+    PricingEngine(std::unique_ptr<StochasticProcess> process) : process_(std::move(process)){}
+
+    virtual ~PricingEngine() = default;
+    virtual double calculate(const Instrument& instrument) const = 0;
+protected:
+    std::unique_ptr<StochasticProcess> process_;
+
+};
+
+
+#endif //QUANTENGINE_PRICINGENGINE_HPP

cpp/RandomGenerator.cpp

@@ -0,0 +1,19 @@
+/**
+ * @file RandomGenerator.cpp
+ * @brief @ref MersenneTwister implementation.
+ */
+
+#include "RandomGenerator.hpp"
+
+
+double MersenneTwister::nextGaussian() {
+    return distr_(generator_);
+}
+
+std::vector<double> MersenneTwister::nextGaussianVector(std::size_t n) {
+    std::vector<double> v(n);
+    for (auto& e : v) {
+        e = nextGaussian();
+    }
+    return v;
+}

cpp/RandomGenerator.hpp

@@ -0,0 +1,31 @@
+/**
+ * @file RandomGenerator.hpp
+ * @brief Random numbers for Monte Carlo (Gaussian draws).
+ */
+
+#ifndef QUANTENGINE_RANDOMGENERATOR_HPP
+#define QUANTENGINE_RANDOMGENERATOR_HPP
+#include <random>
+
+/** @brief Interface for standard normal variates. */
+class RandomGenerator {
+public:
+    RandomGenerator() = default;
+    virtual ~RandomGenerator() = default;
+    virtual double nextGaussian() = 0;
+    virtual std::vector<double> nextGaussianVector(std::size_t n) = 0;
+};
+
+/** @brief @c std::mt19937 with normal distribution. */
+class MersenneTwister : public RandomGenerator {
+public:
+    MersenneTwister() = default;
+    double nextGaussian() override;
+    std::vector<double> nextGaussianVector(std::size_t n) override;
+private:
+    std::mt19937 generator_;
+    std::normal_distribution<> distr_ {0.0, 1.0};
+};
+
+
+#endif //QUANTENGINE_RANDOMGENERATOR_HPP

cpp/Statistics.cpp

@@ -0,0 +1,53 @@
+/**
+ * @file Statistics.cpp
+ * @brief Streaming moment and extrema updates.
+ */
+
+#include "Statistics.hpp"
+
+void Statistics::dump(double value) {
+    for (std::size_t i = 0; i < 3; ++i) {
+        moments_[i] += std::pow(value, i+1);
+    }
+    ++n;
+    max_ = std::max(max_, value);
+    min_ = std::min(min_, value);
+}
+
+void Statistics::clear() {
+    n = 0;
+    moments_ = {0.,0.,0.};
+}
+
+double Statistics::mean() {
+    return moments_[0]/n;
+}
+
+double Statistics::variance() {
+    return moments_[1]/n - std::pow(mean(), 2);
+}
+
+double Statistics::standardDeviation() {
+    return std::sqrt(variance());
+}
+
+double Statistics::skewness() {
+    return moments_[2]/std::pow(n, 3);
+}
+
+
+double Statistics::max() {
+    return max_;
+}
+
+double Statistics::min() {
+    return min_;
+}
+
+double Statistics::sum() {
+    return moments_[0];
+}
+
+double Statistics::count() {
+    return n;
+}

cpp/Statistics.hpp

@@ -0,0 +1,33 @@
+/**
+ * @file Statistics.hpp
+ * @brief Online sample moments for Monte Carlo diagnostics.
+ */
+
+#ifndef QUANTENGINE_STATISTICS_HPP
+#define QUANTENGINE_STATISTICS_HPP
+#include <vector>
+
+/**
+ * @brief Accumulates count, mean/variance-related sums, and running min/max.
+ */
+class Statistics {
+public:
+    Statistics() : moments_({0., 0., 0.}), n(0), max_(0.), min_(0.) {}
+    void dump(double value);
+    void clear();
+    double mean();
+    double variance();
+    double standardDeviation();
+    double skewness();
+    double max();
+    double min();
+    double sum();
+    double count();
+private:
+    std::vector<double> moments_;
+    std::size_t n;
+    double max_, min_;
+};
+
+
+#endif //QUANTENGINE_STATISTICS_HPP

cpp/StochasticProcess.cpp

@@ -0,0 +1,6 @@
+/**
+ * @file StochasticProcess.cpp
+ * @brief @ref StochasticProcess translation unit (interface only).
+ */
+
+#include "StochasticProcess.hpp"

cpp/StochasticProcess.hpp

@@ -0,0 +1,32 @@
+/**
+ * @file StochasticProcess.hpp
+ * @brief Interface for SDE drift, diffusion, and time stepping.
+ */
+
+#ifndef QUANTENGINE_STOCHASTICPROCESS_HPP
+#define QUANTENGINE_STOCHASTICPROCESS_HPP
+#include "MarketData.hpp"
+#include <memory>
+
+/**
+ * @brief Stochastic model for the underlying, driven by @ref MarketData.
+ */
+class StochasticProcess {
+public:
+    StochasticProcess() = delete;
+    explicit StochasticProcess(MarketData data) : data_(std::move(data)){}
+
+    virtual ~StochasticProcess() = default;
+    virtual double drift(double t, double s) = 0;
+    virtual double diffusion(double t, double s) = 0;
+    virtual double step(double t, double s, double dt, double dW) = 0;
+    const MarketData& data() const {return data_;}
+
+
+private:
+    MarketData data_;
+
+};
+
+
+#endif //QUANTENGINE_STOCHASTICPROCESS_HPP

cpp/VolatilitySurface.cpp

@@ -0,0 +1,6 @@
+/**
+ * @file VolatilitySurface.cpp
+ * @brief @ref VolatilitySurface translation unit (interface only).
+ */
+
+#include "VolatilitySurface.hpp"

cpp/VolatilitySurface.hpp

@@ -0,0 +1,28 @@
+/**
+ * @file VolatilitySurface.hpp
+ * @brief Implied volatility as a function of strike and expiry.
+ */
+
+#ifndef QUANTENGINE_VOLATILITYSURFACE_HPP
+#define QUANTENGINE_VOLATILITYSURFACE_HPP
+#include <vector>
+
+/**
+ * @brief Local/vol surface @f$\sigma(K,T)@f$ used by simulation.
+ */
+class VolatilitySurface {
+public:
+    virtual ~VolatilitySurface() = default;
+    virtual double sigma(double K, double T) const = 0;
+private:
+
+};
+
+class SVI : public VolatilitySurface {
+public:
+    SVI() = default;
+    SVI(std::vector<double> K, std::vector<double> rho, std::vector<double> S, std::vector<double> T);
+};
+
+
+#endif //QUANTENGINE_VOLATILITYSURFACE_HPP

cpp/YieldCurve.cpp

@@ -0,0 +1,6 @@
+/**
+ * @file YieldCurve.cpp
+ * @brief @ref YieldCurve translation unit (interface only).
+ */
+
+#include "YieldCurve.hpp"

cpp/YieldCurve.hpp

@@ -0,0 +1,40 @@
+/**
+ * @file YieldCurve.hpp
+ * @brief Abstract yield curve: discount factors and zero rates.
+ */
+
+#ifndef QUANTENGINE_YIELDCURVE_HPP
+#define QUANTENGINE_YIELDCURVE_HPP
+
+/**
+ * @brief Risk-free rate term structure for discounting and risk-neutral drift.
+ */
+class YieldCurve {
+public:
+    YieldCurve() = default;
+
+    YieldCurve(const YieldCurve &other) {
+    }
+
+    YieldCurve(YieldCurve &&other) noexcept {
+    }
+
+    YieldCurve & operator=(const YieldCurve &other) {
+        if (this == &other)
+            return *this;
+        return *this;
+    }
+
+    YieldCurve & operator=(YieldCurve &&other) noexcept {
+        if (this == &other)
+            return *this;
+        return *this;
+    }
+    virtual ~YieldCurve() = default;
+    virtual double discount(double t) const = 0;
+    virtual double zeroRate(double t) const = 0;
+
+};
+
+
+#endif //QUANTENGINE_YIELDCURVE_HPP