Cache the Riemann-Liouville quadrature weights in the rough Heston engine (#2758)

Author: lballabio

ql/math/ode/fractionaladams.hpp

@@ -119,6 +119,72 @@ namespace QuantLib {
     };
 
 
+    //! quadrature weights of the product-trapezoidal Riemann-Liouville integral
+    /*!
+        The weights depend only on the integration order \f$ \alpha \f$ and on
+        the number of grid steps. Callers that integrate repeatedly on the same grid,
+        as the rough Heston characteristic function does once per quadrature node, can
+        build once and reuse.
+    */
+    class RiemannLiouvilleWeights {
+      public:
+        RiemannLiouvilleWeights(Real alpha, Size steps);
+
+        Real alpha() const { return alpha_; }
+
+        //! \f$ I^{\alpha} y(t_N) \f$ for grid values y of spacing dt
+        template <class T>
+        T integrate(const std::vector<T>& y, Real dt) const {
+            QL_REQUIRE(y.size() == steps_ + 1, "grid size (" << y.size()
+                                                             << ") does not match the weights ("
+                                                             << steps_ + 1 << " values expected)");
+            QL_REQUIRE(dt > 0.0, "grid spacing dt (" << dt << ") must be positive");
+
+            // y[steps_] carries weight exactly 1 and is not stored, so w_ is
+            // one shorter than y
+            T sum{0};
+
+            for (Size j{0}; j < steps_; ++j) {
+                sum += w_[j] * y[j];
+            }
+            sum += y[steps_];
+
+            return std::pow(dt, alpha_) / gamma_ * sum;
+        }
+
+      private:
+        Real alpha_{0.0};
+        Size steps_{0};
+        Real gamma_{0.0};
+        std::vector<Real> w_;
+    };
+
+    inline RiemannLiouvilleWeights::RiemannLiouvilleWeights(Real alpha, Size steps) {
+        QL_REQUIRE(alpha >= 0.0, "integration order alpha (" << alpha << ") must be non-negative");
+        QL_REQUIRE(steps > 0, "at least one grid step required");
+
+        alpha_ = alpha;
+        steps_ = steps;
+        gamma_ = GammaFunction().value(alpha + 2.0);
+        w_ = std::vector<Real>(steps);
+
+        // powers m ^ (alpha + 1) for m = 0, ..., steps + 1, reused across the
+        // shifted (k - 1, k, k + 1) terms in the weights below
+        std::vector<Real> powAlphaPlus1(steps + 2);
+        for (Size m{0}; m <= steps + 1; ++m) {
+            powAlphaPlus1[m] = std::pow(static_cast<Real>(m), alpha + 1.0);
+        }
+
+        // weight of y_0: (n - 1) ^ (alpha + 1) - (n - 1 - alpha) n ^ alpha
+        w_[0] = powAlphaPlus1[steps - 1] - (steps - 1 - alpha) * std::pow(Real(steps), alpha);
+
+        for (Size j{1}; j < steps; ++j) {
+            const Size k{steps - j};
+            w_[j] = powAlphaPlus1[k + 1] + powAlphaPlus1[k - 1] - 2.0 * powAlphaPlus1[k];
+        }
+    }
+
+
     //! product-trapezoidal Riemann-Liouville fractional integral
     /*! Approximates
         \f[
@@ -130,31 +196,18 @@ namespace QuantLib {
         piecewise linear interpolant of \f$ y \f$ against the kernel.
         For \f$ \alpha = 1 \f$ this is the trapezoidal rule;
         \f$ \alpha \to 0 \f$ recovers the identity.
+
+        \note The weights are rebuilt on every call.  Callers integrating
+              repeatedly on the same \f$ (\alpha, N) \f$ grid should hold a
+              RiemannLiouvilleWeights instead.
     */
     template <class T = Real>
     T riemannLiouvilleIntegral(const std::vector<T>& y, Real alpha, Real dt) {
         QL_REQUIRE(alpha >= 0.0, "integration order alpha (" << alpha << ") must be non-negative");
         QL_REQUIRE(y.size() >= 2, "at least two grid values required");
         QL_REQUIRE(dt > 0.0, "grid spacing dt (" << dt << ") must be positive");
 
-        const Size n{y.size() - 1};
-
-        // powers m ^ (alpha + 1) for m = 0, ..., n + 1, reused across the shifted
-        // (k - 1, k, k + 1) terms in the quadrature weights below
-        std::vector<Real> powAlphaPlus1(n + 2);
-        for (Size m{0}; m <= n + 1; ++m)
-            powAlphaPlus1[m] = std::pow(static_cast<Real>(m), alpha + 1.0);
-
-        // weight of y_0: (n - 1) ^ (alpha + 1) - (n - 1 - alpha) n ^ alpha
-        T sum{(powAlphaPlus1[n - 1] - (n - 1 - alpha) * std::pow(Real(n), alpha)) * y[0]};
-
-        for (Size j{1}; j < n; ++j) {
-            const Size k{n - j};
-            sum += (powAlphaPlus1[k + 1] + powAlphaPlus1[k - 1] - 2.0 * powAlphaPlus1[k]) * y[j];
-        }
-        sum += y[n];
-
-        return std::pow(dt, alpha) / GammaFunction().value(alpha + 2.0) * sum;
+        return RiemannLiouvilleWeights(alpha, y.size() - 1).integrate(y, dt);
     }
 }
 

ql/pricingengines/vanilla/analyticroughhestonengine.cpp

@@ -178,13 +178,36 @@ namespace QuantLib {
 
     void AnalyticRoughHestonEngine::update() {
         chFCache_.clear();
-        // The kernel nodes depend on the Hurst exponent
+        // The kernel nodes depend on the Hurst exponent.
         liftedGridCache_.clear();
+        // As does the fractional order of the Riemann-Liouville integral
+        rlWeightsOne_.reset();
+        rlWeightsFractional_.reset();
+
         GenericModelEngine<RoughHestonModel,
                            VanillaOption::arguments,
                            VanillaOption::results>::update();
     }
 
+    const RiemannLiouvilleWeights& AnalyticRoughHestonEngine::unitWeights() const {
+        if (!rlWeightsOne_) {
+            rlWeightsOne_.emplace(1.0, timeSteps_);
+        }
+
+        return *rlWeightsOne_;
+    }
+
+    const RiemannLiouvilleWeights& AnalyticRoughHestonEngine::fractionalWeights(Real alpha) const {
+        // Identity check on a cache key, not a numerical comparison: alpha is
+        // rebuilt from the same expression every call, so a difference means
+        // the Hurst exponent moved.
+        if (!rlWeightsFractional_ || rlWeightsFractional_->alpha() != alpha) {
+            rlWeightsFractional_.emplace(alpha, timeSteps_);
+        }
+
+        return *rlWeightsFractional_;
+    }
+
     std::complex<Real> AnalyticRoughHestonEngine::lnChF(
         const std::complex<Real>& z, Time t) const {
 
@@ -244,8 +267,8 @@ namespace QuantLib {
 
         const Real dt{t / timeSteps_};
 
-        return kappa * theta * riemannLiouvilleIntegral(h, 1.0, dt)
-            + v0 * riemannLiouvilleIntegral(h, 1.0 - a, dt);
+        return kappa * theta * unitWeights().integrate(h, dt) +
+               v0 * fractionalWeights(1.0 - a).integrate(h, dt);
     }
 
     std::complex<Real> AnalyticRoughHestonEngine::lnChFPade(
@@ -266,8 +289,8 @@ namespace QuantLib {
             h[j] = evaluatePade(c, sigma * std::pow(Real(j) * dt, a));
         }
 
-        return kappa * theta * riemannLiouvilleIntegral(h, 1.0, dt)
-            + v0 * riemannLiouvilleIntegral(h, 1.0 - a, dt);
+        return kappa * theta * unitWeights().integrate(h, dt) +
+               v0 * fractionalWeights(1.0 - a).integrate(h, dt);
     }
 
     // (3, 3) global rational approximation of Gatheral-Radoicic

ql/pricingengines/vanilla/analyticroughhestonengine.hpp

@@ -27,10 +27,12 @@
 #include <ql/instruments/vanillaoption.hpp>
 #include <ql/math/array.hpp>
 #include <ql/math/integrals/fourierintegration.hpp>
+#include <ql/math/ode/fractionaladams.hpp>
 #include <ql/models/equity/roughhestonmodel.hpp>
 #include <ql/pricingengines/genericmodelengine.hpp>
 #include <complex>
 #include <map>
+#include <optional>
 #include <tuple>
 #include <vector>
 
@@ -204,6 +206,12 @@ namespace QuantLib {
             const ext::shared_ptr<PlainVanillaPayoff>& payoff,
             Time maturity, Real fwd) const;
 
+        //! Weights of the `alpha = 1` integral, built on first use
+        const RiemannLiouvilleWeights& unitWeights() const;
+
+        //! Weights of the `alpha = 1 - a` integral, rebuilt when `a` moves
+        const RiemannLiouvilleWeights& fractionalWeights(Real alpha) const;
+
         const Size timeSteps_;
         const Integration integration_;
         const Real andersenPiterbargEpsilon_, alpha_;
@@ -214,6 +222,14 @@ namespace QuantLib {
         mutable std::map<std::tuple<Real, Real, Time>, std::complex<Real>>
             chFCache_;
         mutable std::map<Time, LiftedGrid> liftedGridCache_;
+
+        /* 
+            Quadrature weights of the two Riemann-Liouville integrals the chF is
+            built from.  Both depend only on (alpha, timeSteps_); the fractional
+            one moves with the Hurst exponent, so both are dropped in `update()`.
+        */
+        mutable std::optional<RiemannLiouvilleWeights> rlWeightsOne_;
+        mutable std::optional<RiemannLiouvilleWeights> rlWeightsFractional_;
     };
 }
 

test-suite/roughhestonmodel.cpp

@@ -21,6 +21,7 @@
 #include "utilities.hpp"
 #include <ql/exercise.hpp>
 #include <ql/instruments/vanillaoption.hpp>
+#include <ql/math/array.hpp>
 #include <ql/math/distributions/gammadistribution.hpp>
 #include <ql/math/ode/fractionaladams.hpp>
 #include <ql/math/ode/fractionalkernelapproximation.hpp>
@@ -38,6 +39,7 @@
 #include <ql/time/daycounters/actual365fixed.hpp>
 #include <cmath>
 #include <complex>
+#include <iomanip>
 #include <vector>
 
 using namespace QuantLib;
@@ -188,6 +190,144 @@ BOOST_AUTO_TEST_CASE(testRiemannLiouvilleIntegral) {
     }
 }
 
+BOOST_AUTO_TEST_CASE(testRiemannLiouvilleWeightReuse) {
+    BOOST_TEST_MESSAGE("Testing reusable Riemann-Liouville quadrature weights...");
+
+    // testRiemannLiouvilleIntegral above already validates the arithmetic
+    // against the closed form, and now reaches it through this class, so it is
+    // not repeated here.  What is new is that one object is integrated many
+    // times: it has to stay stateless, so that a cached set of weights gives
+    // exactly what a set built for that one call would have given.
+    //
+    // steps = 1 leaves the interior loop empty, steps = 2 runs it once and
+    // reaches index 0 of the power table, steps = 256 is the engine default.
+    for (const Size steps : {Size(1), Size(2), Size(256)}) {
+        std::vector<Real> real(steps + 1);
+        std::vector<std::complex<Real>> complex(steps + 1);
+
+        for (Size i{0}; i <= steps; ++i) {
+            const Real x{Real(i) / steps};
+            real[i] = std::exp(-0.7 * x) * (1.0 + x * x);
+            complex[i] = std::complex<Real>(real[i], std::sin(37.0 * x));
+        }
+
+        // 0 and 1 are the degenerate orders; 0.25 is the fractional one the
+        // engine asks for at H = 0.25
+        for (const Real alpha : {0.0, 0.25, 1.0}) {
+            const RiemannLiouvilleWeights reused(alpha, steps);
+
+            // the spacing varies with the maturity while the weights do not,
+            // so a reused object has to survive a change of dt
+            for (const Real dt : {1e-3, 0.5}) {
+                const Real calculated{reused.integrate(real, dt)};
+                const Real expected{riemannLiouvilleIntegral(real, alpha, dt)};
+
+                if (value(calculated) != value(expected))
+                    BOOST_ERROR("reused weights do not reproduce freshly built ones"
+                                << "\n    alpha:      " << alpha << "\n    steps:      " << steps
+                                << "\n    dt:         " << dt
+                                << "\n    calculated: " << std::setprecision(17) << calculated
+                                << "\n    expected:   " << std::setprecision(17) << expected);
+
+                // production only ever integrates complex vectors
+                const std::complex<Real> calculatedC{reused.integrate(complex, dt)};
+                const std::complex<Real> expectedC{
+                    riemannLiouvilleIntegral(complex, alpha, dt)};
+
+                if (value(calculatedC.real()) != value(expectedC.real())
+                    || value(calculatedC.imag()) != value(expectedC.imag()))
+                    BOOST_ERROR("reused weights do not reproduce freshly built ones "
+                                "on the complex path"
+                                << "\n    alpha:      " << alpha << "\n    steps:      " << steps
+                                << "\n    dt:         " << dt
+                                << "\n    calculated: " << std::setprecision(17) << calculatedC
+                                << "\n    expected:   " << std::setprecision(17) << expectedC);
+            }
+        }
+    }
+
+    const Size steps{64};
+    const Real dt{0.01};
+
+    std::vector<Real> y(steps + 1);
+    for (Size i{0}; i <= steps; ++i)
+        y[i] = 1.0 + std::sin(Real(i));
+
+    // An oracle for the weights themselves rather than for the integral: at
+    // alpha = 1 they have to be exactly the trapezoidal rule, which pins down
+    // the formula more tightly than agreement with I^a t^p to 5e-5 does.
+    Real trapezoid{0.5 * (y[0] + y[steps])};
+    for (Size i{1}; i < steps; ++i)
+        trapezoid += y[i];
+    trapezoid *= dt;
+
+    QL_CHECK_CLOSE(RiemannLiouvilleWeights(1.0, steps).integrate(y, dt), trapezoid, 1e-10);
+
+    // steps = 0 would underflow (steps - 1) into an out-of-bounds index, and a
+    // grid that does not match the weights is the mistake a caller reusing a
+    // cached object will make
+    BOOST_CHECK_THROW(RiemannLiouvilleWeights(0.6, 0), Error);
+    BOOST_CHECK_THROW(RiemannLiouvilleWeights(-1e-8, steps), Error);
+    BOOST_CHECK_THROW(RiemannLiouvilleWeights(0.6, steps - 1).integrate(y, dt), Error);
+    BOOST_CHECK_THROW(RiemannLiouvilleWeights(0.6, steps).integrate(y, 0.0), Error);
+}
+
+BOOST_AUTO_TEST_CASE(testWeightCacheInvalidation) {
+    BOOST_TEST_MESSAGE("Testing that the rough Heston engine drops cached quadrature weights "
+                       "when the Hurst exponent moves...");
+
+    const Date today(2, July, 2026);
+    Settings::instance().evaluationDate() = today;
+    const DayCounter dc{Actual365Fixed()};
+
+    const Handle<YieldTermStructure> rTS(ext::make_shared<FlatForward>(today, 0.03, dc));
+    const Handle<YieldTermStructure> qTS(ext::make_shared<FlatForward>(today, 0.01, dc));
+    const Handle<Quote> s0(ext::make_shared<SimpleQuote>(100.0));
+
+    const Real v0{0.04}, kappa{1.5}, theta{0.04}, sigma{0.3}, rho{-0.7};
+
+    const auto makeModel{[&](Real hurst) {
+        return ext::make_shared<RoughHestonModel>(
+            ext::make_shared<HestonProcess>(rTS, qTS, s0, v0, kappa, theta, sigma, rho), hurst);
+    }};
+
+    const auto payoff{ext::make_shared<PlainVanillaPayoff>(Option::Call, 100.0)};
+    const Time t{1.0};
+
+    for (const auto approximation :
+         {AnalyticRoughHestonEngine::Approximation::AdamsPredictorCorrector,
+          AnalyticRoughHestonEngine::Approximation::Pade}) {
+
+        const auto model{makeModel(0.1)};
+        const auto engine{
+            ext::make_shared<AnalyticRoughHestonEngine>(model, 128, 64, approximation)};
+
+        const Real firstPrice{engine->priceVanillaPayoff(payoff, t)};
+
+        // The fractional weights are keyed on 1 - (H + 1/2); moving H has to
+        // rebuild them, or the engine keeps pricing the old model.
+        Array params{model->params()};
+        params[5] = 0.3;
+        model->setParams(params);
+
+        const Real movedPrice{engine->priceVanillaPayoff(payoff, t)};
+
+        const auto reference{
+            ext::make_shared<AnalyticRoughHestonEngine>(makeModel(0.3), 128, 64, approximation)};
+        const Real referencePrice{reference->priceVanillaPayoff(payoff, t)};
+
+        if (std::fabs(value(movedPrice) - value(referencePrice)) > 1e-12)
+            BOOST_ERROR("stale quadrature weights after a change of the Hurst exponent"
+                        << "\n    reused engine: " << movedPrice
+                        << "\n    fresh engine:  " << referencePrice);
+
+        if (std::fabs(value(movedPrice) - value(firstPrice)) < 1e-6)
+            BOOST_ERROR("the Hurst exponent did not move the price, so the test "
+                        "cannot detect a stale cache"
+                        << "\n    H = 0.1: " << firstPrice << "\n    H = 0.3: " << movedPrice);
+    }
+}
+
 BOOST_AUTO_TEST_CASE(testEquivalenceWithHestonModel) {
     BOOST_TEST_MESSAGE(
         "Testing rough Heston engine against the classical Heston engine for H = 0.5...");