QPR-12490 bug fixes and stabilization for sabr model

pcaspers · jenkins · commit 0fa1bc9e3411 · 2024-05-03T06:43:49.000Z
diff --git a/Docker/.env b/Docker/.env
@@ -3,7 +3,7 @@ DEBIAN_TAG=11.7
 # it's recommended to include CMAKE_BUILD_TYPE and BOOST_VARIANT in QL_TAG, ORE_TAG, BOOST_TAG
 # to distinguish a release build from a debug build
 
-QL_TAG=1.31.1_ddda9430
+QL_TAG=1.31.1_f8955cf0
 ORE_TAG=latest
 
 # debug or release
diff --git a/Examples/Example_60/ExpectedOutput/npv.csv b/Examples/Example_60/ExpectedOutput/npv.csv
@@ -1,3 +1,3 @@
 #TradeId,TradeType,Maturity,MaturityTime,NPV,NpvCurrency,NPV(Base),BaseCurrency,Notional,NotionalCurrency,Notional(Base),NettingSet,CounterParty
 Cap_USD_SOFR,Swap,2025-03-21,1.003324,0.000931,USD,0.000931,USD,100000000.00,USD,100000000.00,DUMMY_NS,DUMMY_CP
-Swaption_USD_SOFR,Swaption,2046-01-12,21.817022,250911.447892,USD,250911.447892,USD,100000000.00,USD,100000000.00,DUMMY_NS,DUMMY_CP
+Swaption_USD_SOFR,Swaption,2046-01-12,21.817022,250911.447904,USD,250911.447904,USD,100000000.00,USD,100000000.00,DUMMY_NS,DUMMY_CP
diff --git a/QuantExt/qle/termstructures/parametricvolatility.cpp b/QuantExt/qle/termstructures/parametricvolatility.cpp
@@ -57,20 +57,21 @@ Real ParametricVolatility::convert(const Real inputQuote, const MarketQuoteType
     Real forwardPremium;
     switch (marketModelType_) {
     case MarketModelType::Black76:
-        if (inputMarketQuoteType_ == MarketQuoteType::Price) {
+        if (inputMarketQuoteType == MarketQuoteType::Price) {
             forwardPremium = inputQuote / (discountCurve_.empty() ? 1.0 : discountCurve_->discount(timeToExpiry));
-        } else if (inputMarketQuoteType_ == MarketQuoteType::NormalVolatility) {
+        } else if (inputMarketQuoteType == MarketQuoteType::NormalVolatility) {
             forwardPremium =
                 bachelierBlackFormula(inputOptionType, strike, forward, inputQuote * std::sqrt(timeToExpiry));
-        } else if (inputMarketQuoteType_ == MarketQuoteType::ShiftedLognormalVolatility) {
+        } else if (inputMarketQuoteType == MarketQuoteType::ShiftedLognormalVolatility) {
             if (strike < -inputLognormalShift)
                 forwardPremium = inputOptionType == Option::Call ? forward - strike : 0.0;
-            else
+            else {
                 forwardPremium = blackFormula(inputOptionType, strike, forward, inputQuote * std::sqrt(timeToExpiry),
-                                              inputLognormalShift);
+                                              1.0, inputLognormalShift);
+            }
         } else {
             QL_FAIL("ParametricVolatility::convert(): MarketQuoteType ("
-                    << static_cast<int>(inputMarketQuoteType_) << ") not handled. Internal error, contact dev.");
+                    << static_cast<int>(inputMarketQuoteType) << ") not handled. Internal error, contact dev.");
         }
         break;
     default:
diff --git a/QuantExt/qle/termstructures/sabrparametricvolatility.cpp b/QuantExt/qle/termstructures/sabrparametricvolatility.cpp
@@ -123,27 +123,22 @@ std::vector<std::pair<Real, bool>> SabrParametricVolatility::defaultModelParamet
 std::vector<Real> SabrParametricVolatility::direct(const std::vector<Real>& x, const Real forward,
                                                    const Real lognormalShift) const {
     std::vector<Real> y(4);
-    // beta as in QuantLib::SABRSpecs
-    y[1] = std::fabs(x[1]) < std::sqrt(-std::log(eps1)) ? std::exp(-(x[1] * x[1])) : eps1;
-    // max 0.02 normal vol equivalent
-    Real fbeta = std::pow(forward + lognormalShift, y[1]);
-    y[0] =
-        (std::fabs(x[0]) < std::sqrt(-std::log(fbeta * eps1 / 0.02)) ? 0.02 * std::exp(-(x[0] * x[0])) / fbeta : eps1);
-    // nu max 5.0
-    y[2] = (std::fabs(x[2]) < std::sqrt(-std::log(2.0 * eps1)) ? 2.0 * std::exp(-(x[2] * x[2])) : eps1);
-    // rho as in QuantLib::SABRSpecs
+    y[1] = std::max(eps1, std::min(1.0 - eps1, std::exp(-(x[1] * x[1]))));
+    Real fbeta = std::pow(std::max(forward + lognormalShift, eps1), y[1]);
+    y[0] = std::max(eps1, std::exp(-(x[0] * x[0])) / fbeta * max_nvol_equiv);
+    y[2] = std::max(eps1, std::exp(-(x[2] * x[2])) * max_nu);
     y[3] = std::fabs(x[3]) < 2.5 * M_PI ? eps2 * std::sin(x[3]) : eps2 * (x[3] > 0.0 ? 1.0 : (-1.0));
     return y;
 }
 
 std::vector<Real> SabrParametricVolatility::inverse(const std::vector<Real>& y, const Real forward,
                                                     const Real lognormalShift) const {
     std::vector<Real> x(4);
-    x[1] = std::sqrt(-std::log(y[1]));
-    Real fbeta = std::pow(forward + lognormalShift, y[1]);
-    x[0] = std::sqrt(-std::log(y[0] * fbeta / 0.02));
-    x[2] = std::sqrt(-std::log(y[2] / 2.0));
-    x[3] = std::asin(y[3] / eps2);
+    x[1] = std::sqrt(-std::log(std::min(1.0 - eps1, std::max(eps1, y[1]))));
+    Real fbeta = std::pow(std::max(forward + lognormalShift, eps1), y[1]);
+    x[0] = std::sqrt(-std::log(std::min(1.0 - eps1, std::max(eps1, y[0] * fbeta / max_nvol_equiv))));
+    x[2] = std::sqrt(-std::log(std::min(1.0 - eps1, std::max(eps1, y[2] / max_nu))));
+    x[3] = std::asin(std::max(-eps2, std::min(eps2, y[3])));
     return x;
 }
 
@@ -267,9 +262,9 @@ SabrParametricVolatility::calibrateModelParameters(const MarketSmile& marketSmil
 
     // if there are no free parameters, we just pass back the fixed parameters as the result
 
-    if(noFreeParams == 0) {
+    if (noFreeParams == 0) {
         std::vector<Real> resultParams;
-        for(auto const& p : params)
+        for (auto const& p : params)
             resultParams.push_back(p.first);
         return std::make_tuple(resultParams, 0.0, 0);
     }
@@ -286,6 +281,7 @@ SabrParametricVolatility::calibrateModelParameters(const MarketSmile& marketSmil
         Real lognormalShift_;
         std::vector<Real> strikes_;
         std::vector<Real> marketQuotes_;
+        Real refQuote_;
         std::function<std::vector<Real>(const std::vector<Real>&, const Real, const Real, const Real,
                                         const std::vector<Real>&)>
             evalSabr_;
@@ -310,23 +306,10 @@ SabrParametricVolatility::calibrateModelParameters(const MarketSmile& marketSmil
             Array result(strikes_.size());
             auto sabr = evalSabr(x);
             for (Size i = 0; i < strikes_.size(); ++i) {
-                result[i] = (marketQuotes_[i] - sabr[i]);
+                result[i] = (marketQuotes_[i] - sabr[i]) / refQuote_;
             }
             return result;
         }
-
-        Real normalizedError(const Array& x) const {
-            auto sabr = evalSabr(x);
-            Real maxQuote = *std::max_element(marketQuotes_.begin(), marketQuotes_.end());
-            Real result = 0.0;
-            for (Size i = 0; i < strikes_.size(); ++i) {
-                /* we want a relative error measure, but can't do this point by point, because for far-ootm strikes
-                 * premiums are close to zero */
-                result += std::pow((marketQuotes_[i] - sabr[i]) / maxQuote, 2.0);
-            }
-            result = std::sqrt(result / static_cast<Real>(strikes_.size()));
-            return result;
-        }
     };
 
     /* build the target function and populate the members:
@@ -364,6 +347,8 @@ SabrParametricVolatility::calibrateModelParameters(const MarketSmile& marketSmil
             marketSmile.timeToExpiry, marketSmile.strikes[i], marketSmile.forward, preferredOutputQuoteType(),
             marketSmile.lognormalShift, boost::none));
     }
+    // we use relative errors w.r.t. the max market quote, because far otm quotes are close to zero
+    t.refQuote_ = *std::max_element(t.marketQuotes_.begin(), t.marketQuotes_.end());
 
     // perform the calibration (this step might throw if all minimizations go wrong)
 
@@ -406,7 +391,7 @@ SabrParametricVolatility::calibrateModelParameters(const MarketSmile& marketSmil
             continue;
         }
 
-        Real thisError = t.normalizedError(problem.currentValue());
+        Real thisError = problem.functionValue();
         if (thisError < bestError) {
             bestError = thisError;
             for (Size i = 0, j = 0; i < bestResult.size(); ++i) {
@@ -427,6 +412,17 @@ SabrParametricVolatility::calibrateModelParameters(const MarketSmile& marketSmil
     return std::make_tuple(bestResult, bestError, ++attempt);
 }
 
+namespace {
+void laplaceInterpolationWithErrorHandling(Matrix& m, const std::vector<Real>& x, const std::vector<Real>& y) {
+    try {
+        laplaceInterpolation(m, x, y, 1E-6, 100);
+    } catch (const std::exception& e) {
+        QL_FAIL("Error during laplaceInterpolation() in SabrParametricVolatility: "
+                << e.what() << ", this might be related to the numerical parameters relTol, maxIterMult. Contact dev.");
+    }
+}
+} // namespace
+
 void SabrParametricVolatility::calculate() {
 
     // if no model parameters are given, we provide the default ones
@@ -460,6 +456,7 @@ void SabrParametricVolatility::calculate() {
             calibrationErrors_[key] = error;
             noOfAttempts_[key] = noOfAttempts;
         } catch (const std::exception& e) {
+            std::cerr << "error: " << e.what() << std::endl;
             // all calibration failed -> do not populate params, but interpolate them below
         }
         lognormalShifts_[key] = s.lognormalShift;
@@ -515,10 +512,10 @@ void SabrParametricVolatility::calculate() {
 
     // interpolate the null values
 
-    laplaceInterpolation(alpha_, timeToExpiries_, underlyingLengths_);
-    laplaceInterpolation(beta_, timeToExpiries_, underlyingLengths_);
-    laplaceInterpolation(nu_, timeToExpiries_, underlyingLengths_);
-    laplaceInterpolation(rho_, timeToExpiries_, underlyingLengths_);
+    laplaceInterpolationWithErrorHandling(alpha_, timeToExpiries_, underlyingLengths_);
+    laplaceInterpolationWithErrorHandling(beta_, timeToExpiries_, underlyingLengths_);
+    laplaceInterpolationWithErrorHandling(nu_, timeToExpiries_, underlyingLengths_);
+    laplaceInterpolationWithErrorHandling(rho_, timeToExpiries_, underlyingLengths_);
 
     // sanitize values produced by the the interpolation that are not allowed
 
diff --git a/QuantExt/qle/termstructures/sabrparametricvolatility.hpp b/QuantExt/qle/termstructures/sabrparametricvolatility.hpp
@@ -74,6 +74,8 @@ class SabrParametricVolatility final : public ParametricVolatility {
 private:
     static constexpr double eps1 = .0000001;
     static constexpr double eps2 = .9999;
+    static constexpr double max_nvol_equiv = 0.02;
+    static constexpr double max_nu = 2.0;
 
     void calculate();
 
diff --git a/QuantExt/qle/termstructures/swaptionsabrcube.cpp b/QuantExt/qle/termstructures/swaptionsabrcube.cpp
@@ -46,6 +46,17 @@ SwaptionSabrCube::SwaptionSabrCube(
             registerWith(s);
 }
 
+namespace {
+void laplaceInterpolationWithErrorHandling(Matrix& m, const std::vector<Real>& x, const std::vector<Real>& y) {
+    try {
+        laplaceInterpolation(m, x, y, 1e-6, 100);
+    } catch (const std::exception& e) {
+        QL_FAIL("Error during laplaceInterpolation() in SwaptionSabrCube: "
+                << e.what() << ", this might be related to the numerical parameters relTol, maxIterMult. Contact dev.");
+    }
+}
+} // namespace
+
 void SwaptionSabrCube::performCalculations() const {
 
     SwaptionVolatilityCube::performCalculations();
@@ -81,7 +92,7 @@ void SwaptionSabrCube::performCalculations() const {
     }
 
     for (auto& v : interpolatedVolSpreads) {
-        laplaceInterpolation(v, allOptionTimes, allSwapLengths);
+        laplaceInterpolationWithErrorHandling(v, allOptionTimes, allSwapLengths);
     }
 
     // build market smiles on the grid
diff --git a/QuantLib b/QuantLib
@@ -1 +1 @@
-Subproject commit ddda9430217b4612f55dc1a2fd48ba7c90976fed
+Subproject commit f8955cf02565a45c60978c8ea3558615bb8a6dc7
