Feat: sabr interpolation using Ceres solver

Author: bersen66

cdr/CMakeLists.txt

@@ -8,3 +8,13 @@ add_subdirectory(options)
 add_subdirectory(fx)
 add_subdirectory(market)
 add_subdirectory(model)
+
+
+cdr_cpp_test(
+  NAME ceres_tests
+  SRCS
+    "ceres_tests.cc"
+  DEPS
+    Ceres::ceres
+    GTest::gtest_main
+)

cdr/ceres_tests.cc

@@ -0,0 +1,95 @@
+// Check that ceres works correctly
+#include <ceres/ceres.h>
+#include <gtest/gtest.h>
+
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+template <typename T>
+T CalculateSABRVol(const T& alpha, const T& rho, const T& nu, double F, double K, double T_exp) {
+    using std::abs;
+    using std::log;
+    using std::pow;
+    using std::sqrt;
+
+    const double beta = 1.0;
+    const double eps = 1e-7;
+
+    if (abs(F - K) < eps) {
+        T sub_1 = T((1.0 - beta) * (1.0 - beta) / 24.0) * alpha * alpha / T(pow(F, 2.0 - 2.0 * beta));
+        T sub_2 = T(0.25 * rho * beta * nu) * alpha / T(pow(F, 1.0 - beta));
+        T sub_3 = (T(2.0) - T(3.0) * rho * rho) / T(24.0) * nu * nu;
+        return (alpha / T(pow(F, 1.0 - beta))) * (T(1.0) + (sub_1 + sub_2 + sub_3) * T(T_exp));
+    } else {
+        T logFK = T(log(F / K));  // Явное приведение double -> T
+        T f_k_beta = T(pow(F * K, (1.0 - beta) / 2.0));
+
+        T z = (nu / alpha) * f_k_beta * logFK;
+        T x_z = log((sqrt(T(1.0) - T(2.0) * rho * z + z * z) + z - rho) / (T(1.0) - rho));
+
+        T numer_1 = T((1.0 - beta) * (1.0 - beta) / 24.0) * alpha * alpha / pow(f_k_beta, 2.0);
+        T numer_2 = T(0.25) * rho * T(beta) * nu * alpha / f_k_beta;
+        T numer_3 = (T(2.0) - T(3.0) * rho * rho) / T(24.0) * nu * nu;
+
+        T denum_1 = T((1.0 - beta) * (1.0 - beta) / 24.0) * logFK * logFK;
+        T denum_2 = T(pow((1.0 - beta), 4.0) / 1920.0) * pow(logFK, 4.0);
+
+        return (alpha / (f_k_beta * (T(1.0) + denum_1 + denum_2))) * (z / x_z) *
+               (T(1.0) + (numer_1 + numer_2 + numer_3) * T(T_exp));
+    }
+}
+
+// 2. Функтор для Ceres
+struct SabrCostFunctor {
+    SabrCostFunctor(double market_vol, double F, double K, double T_exp)
+        : m_vol(market_vol), m_F(F), m_K(K), m_T(T_exp) {
+    }
+
+    template <typename T>
+    bool operator()(const T* const alpha, const T* const rho, const T* const nu, T* residual) const {
+        // Вычисляем волатильность по модели
+        T model_vol = CalculateSABRVol(alpha[0], rho[0], nu[0], m_F, m_K, m_T);
+        // Ошибка (residual)
+        residual[0] = model_vol - T(m_vol);
+        return true;
+    }
+
+    const double m_vol, m_F, m_K, m_T;
+};
+
+TEST(SABRProvider, CalibrationWithJet) {
+    const double F = 100.0;
+    const double K = 105.0;  // Немного вне денег
+    const double T = 1.0;
+    const double market_vol = 0.22;
+
+    double alpha = 0.2;
+    double rho = -0.2;
+    double nu = 0.4;
+
+    ceres::Problem problem;
+
+    problem.AddResidualBlock(
+        new ceres::AutoDiffCostFunction<SabrCostFunctor, 1, 1, 1, 1>(new SabrCostFunctor(market_vol, F, K, T)), nullptr,
+        &alpha, &rho, &nu);
+
+    problem.SetParameterLowerBound(&alpha, 0, 0.001);
+    problem.SetParameterLowerBound(&rho, 0, -0.99);
+    problem.SetParameterUpperBound(&rho, 0, 0.99);
+    problem.SetParameterLowerBound(&nu, 0, 0.001);
+
+    ceres::Solver::Options options;
+    options.max_num_iterations = 50;
+    options.linear_solver_type = ceres::DENSE_QR;
+
+    ceres::Solver::Summary summary;
+    ceres::Solve(options, &problem, &summary);
+
+    std::cout << "Solver report:\n" << summary.BriefReport() << "\n";
+    std::cout << "Final params: a=" << alpha << ", r=" << rho << ", n=" << nu << "\n";
+
+    EXPECT_TRUE(summary.IsSolutionUsable());
+    double final_vol = CalculateSABRVol(alpha, rho, nu, F, K, T);
+    EXPECT_NEAR(final_vol, market_vol, 1e-6);
+}

cdr/options/CMakeLists.txt

@@ -10,18 +10,20 @@ cdr_cpp_library(
     "interpolation/quadratic_spline.h"
     "interpolation/lerp.h"
     "interpolation/flat_forward.h"
+    "interpolation/sabr.h"
   SRCS
     "contract.cc"
-    #"volatility.cc"
     "interpolation/cubic_spline.cc"
     "interpolation/quadratic_spline.cc"
     "interpolation/lerp.cc"
     "interpolation/flat_forward.cc"
+    "interpolation/sabr.cc"
   DEPS
     cdr::calendar
     cdr::base
     cdr::types
     cdr::math
+    Ceres::ceres
   PUBLIC
 )
 
@@ -47,6 +49,16 @@ cdr_cpp_test(
     GTest::gtest_main
 )
 
+cdr_cpp_test(
+  NAME sabr_tests
+  SRCS
+    "interpolation/sabr_tests.cc"
+  DEPS
+    cdr::options
+    Ceres::ceres
+    GTest::gtest_main
+)
+
 cdr_cpp_test(
   NAME option_tests
   SRCS

cdr/options/interpolation/sabr.cc

@@ -0,0 +1,52 @@
+#include <cdr/options/interpolation/sabr.h>
+#include <ceres/ceres.h>
+
+namespace cdr {
+
+Expect<void, Error> SABRInterpolator::InitState(void* coefs_ptr, std::span<const f64> xs,
+                                                std::span<const f64> ys) noexcept {
+    auto* p = static_cast<SabrParameters*>(coefs_ptr);
+
+    // Начальное приближение (Initial Guess)
+    // alpha берем как ATM волатильность (примерно посередине вектора ys)
+    double params[3] = {ys[ys.size() / 2], 0.0, 0.5};
+
+    ceres::Problem problem;
+    for (size_t i = 0; i < xs.size(); ++i) {
+        auto* cost_function =
+            new ceres::AutoDiffCostFunction<SabrCostFunctor, 1, 3>(new SabrCostFunctor(xs[i], ys[i], p->F, p->T));
+        problem.AddResidualBlock(cost_function, nullptr, params);
+    }
+
+    // Накладываем ограничения, чтобы параметры не улетели в космос
+    problem.SetParameterLowerBound(params, 0, 1e-6);   // alpha > 0
+    problem.SetParameterLowerBound(params, 1, -0.99);  // rho > -1
+    problem.SetParameterUpperBound(params, 1, 0.99);   // rho < 1
+    problem.SetParameterLowerBound(params, 2, 1e-6);   // nu > 0
+
+    ceres::Solver::Options options;
+    options.linear_solver_type = ceres::DENSE_QR;
+    options.max_num_iterations = 50;
+    // Levenberg-Marquardt по умолчанию в Ceres
+
+    ceres::Solver::Summary summary;
+    ceres::Solve(options, &problem, &summary);
+
+    if (!summary.IsSolutionUsable()) {
+        return ErrorCalibrationFailed();
+    }
+
+    p->alpha = params[0];
+    p->rho = params[1];
+    p->nu = params[2];
+
+    return Ok();
+}
+
+[[nodiscard]] Expect<f64, Error> SABRInterpolator::Evaluate(const void* coefs_ptr,
+                                                            [[maybe_unused]] std::span<const f64> xs, f64 K) noexcept {
+    const auto& p = *static_cast<const SabrParameters*>(coefs_ptr);
+    return Ok(CalculateHagan<f64>(K, p.F, p.T, p.alpha, p.rho, p.nu));
+}
+
+}  // namespace cdr

cdr/options/interpolation/sabr.h

@@ -0,0 +1,112 @@
+#pragma once
+
+#include <cdr/options/internal/export.h>
+#include <cdr/types/errors.h>
+#include <cdr/types/types.h>
+
+#include <span>
+
+namespace cdr {
+
+class CDR_OPTIONS_EXPORT SABRInterpolator {
+public:
+    struct SabrParameters {
+        static constexpr f64 kBeta = 1.0;
+        f64 alpha;
+        f64 rho;
+        f64 nu;
+        f64 F;  // Forward price
+        f64 T;  // Time to maturity
+        f64 rd; // domestic curve
+        f64 rf; // foreign curve
+    };
+
+    using StrikeVolatilityType = SabrParameters;
+    using DeltaVolatilityType = f64;
+
+public:
+    static constexpr size_t StateRequiredMemorySize(size_t /*n*/) noexcept {
+        // SABR needs only one structure of parameters per date
+        return sizeof(SabrParameters);
+    }
+
+    static constexpr size_t StateRequiredMemoryAlignment(size_t /*n*/) noexcept {
+        return alignof(SabrParameters);
+    }
+
+    // Hagan approximation, implemented as template to support double and ceres::Jet
+    // Uses the Hagan approximation for the SABR model with Beta = 1.0
+    //
+    // Formula: sigma(K) = alpha * (z / chi(z)) * [1 + correction_term * T]
+    // where:
+    //   z = (nu / alpha) * log(F / K)
+    //   chi(z) = log((sqrt(1 - 2*rho*z + z^2) + z - rho) / (1 - rho))
+    //
+    // Note: Handle K = F (ATM) case separately to avoid 0/0 singularity,
+    // resulting in: sigma_atm = alpha * [1 + correction_term * T].
+    template <typename T>
+    static T CalculateHagan(T K, T F, T Time, T alpha, T rho, T nu) {
+        using std::log;
+        using std::sqrt;
+        using std::abs;
+        using std::pow;
+
+        if (K <= T(0.0) || F <= T(0.0)) {
+            return T(0.0);
+        }
+
+        constexpr f64 beta = SabrParameters::kBeta;
+        constexpr f64 one_minus_beta = 1.0 - beta;
+        constexpr f64 one_minus_beta2 = one_minus_beta * one_minus_beta;
+
+        const T f_k_beta = pow(F * K, one_minus_beta / 2.0);
+        const T log_f_k = log(F / K);
+
+        T den = f_k_beta * (T(1.0) + (one_minus_beta2 / 24.0) * log_f_k * log_f_k +
+                            (one_minus_beta2 * one_minus_beta2 / 1920.0) * pow(log_f_k, 4.0));
+
+        const T z = (nu / alpha) * f_k_beta * log_f_k;
+
+        const T xz = log((sqrt(T(1.0) - T(2.0) * rho * z + z * z) + z - rho) / (T(1.0) - rho));
+
+        const T term2 = T(1.0) + (
+            (one_minus_beta2 / 24.0) * (alpha * alpha / pow(f_k_beta, 2.0)) +
+            (0.25 * rho * beta * nu * alpha / f_k_beta) +
+            ((T(2.0) - T(3.0) * rho * rho) / 24.0) * nu * nu
+        ) * Time;
+
+        if (abs(log_f_k) < T(1e-8)) {
+            return (alpha / pow(F, one_minus_beta)) * term2;
+        }
+
+        return (alpha / den) * (z / xz) * term2;
+    }
+
+    [[nodiscard]] static Expect<void, Error> InitState(void* coefs_ptr, std::span<const f64> xs,
+                                                       std::span<const f64> ys) noexcept;
+
+    [[nodiscard]] static Expect<f64, Error> Evaluate(const void* coefs_ptr, std::span<const f64> xs, f64 x) noexcept;
+};
+
+struct CDR_OPTIONS_EXPORT SabrCostFunctor {
+    SabrCostFunctor(f64 strike, f64 market_vol, f64 F, f64 T)
+        : strike_(strike), market_vol_(market_vol), F_(F), T_(T)
+    {}
+
+    template <typename T>
+    bool operator()(const T* const params, T* residual) const {
+        T model_vol = SABRInterpolator::CalculateHagan<T>(T(strike_), T(F_), T(T_), params[0], params[1], params[2]);
+
+        residual[0] = model_vol - T(market_vol_);
+
+        return true;
+    }
+
+private:
+    const f64 strike_;
+    const f64 market_vol_;
+    const f64 F_;
+    const f64 T_;
+};
+
+}  // namespace cdr