Merge 006e5c5 into b05e391

lballabio · Jun 9, 2020 · bc4f9c0 · bc4f9c0
2 parents b05e391 + 006e5c5
commit bc4f9c0
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 36 deletions.
diff --git a/ql/methods/finitedifferences/utilities/localvolrndcalculator.cpp b/ql/methods/finitedifferences/utilities/localvolrndcalculator.cpp
@@ -177,23 +177,28 @@ namespace QuantLib {
         calculate();
 
         const Time closeGridTime(timeGrid_->closestTime(t));
+
         if (closeGridTime == 0.0) {
+            const Real stepSize = 0.02*(
+                    xm_[0]->locations().back() - xm_[0]->locations().front());
             return RiskNeutralDensityCalculator::InvCDFHelper(
-                this, std::log(spot_->value()), 0.1*localVolProbEps_, maxIter_)
-                .inverseCDF(p, t);
+                this, std::log(spot_->value()),
+                0.1*localVolProbEps_, maxIter_, stepSize).inverseCDF(p, t);
         }
         else {
             Array xp(xGrid_);
             const Size idx = timeGrid_->index(closeGridTime)-1;
 
             const Array x(xm_[idx]->locations().begin(),
                           xm_[idx]->locations().end());
+            const Real stepSize = 0.005*(x.back() - x.front());
+
             std::transform(x.begin(), x.end(), pm_->row_begin(idx), xp.begin(),
                            std::multiplies<Real>());
 
             const Real xm = DiscreteSimpsonIntegral()(x, xp);
             return RiskNeutralDensityCalculator::InvCDFHelper(
-                this, xm, 0.1*localVolProbEps_, maxIter_).inverseCDF(p, t);
+                this, xm, 0.1*localVolProbEps_, maxIter_, stepSize).inverseCDF(p, t);
         }
     }
 
@@ -231,6 +236,7 @@ namespace QuantLib {
 
         const Volatility stdDevOfFirstStep = vol * std::sqrt(sT);
         const Real normInvEps = InverseCumulativeNormal()(1 - localVolProbEps_);
+
         Real sLowerBound = xm - normInvEps * stdDevOfFirstStep;
         Real sUpperBound = xm + normInvEps * stdDevOfFirstStep;
 
@@ -287,18 +293,12 @@ namespace QuantLib {
                 const Real vm = DiscreteSimpsonIntegral()(x, vols)
                     /(x.back() - x.front());
 
-                const Size scalingSteps
-                    = std::max(3, Integer(0.01*timeGrid_->size()));
-
-                const Real scalingFactor
-                    = std::sqrt(timeGrid_->at(
-                        std::min(timeGrid_->size()-1, i+scalingSteps)))
-                     * vm;
+                const Real scalingFactor = vm*std::sqrt(0.5*timeGrid_->back());
 
                 if (maxLeftValue > localVolProbEps_)
-                    sLowerBound -= scalingFactor*xm;
+                    sLowerBound -= scalingFactor*(oldUpperBound-oldLowerBound);
                 if (maxRightValue > localVolProbEps_)
-                    sUpperBound += scalingFactor*xm;
+                    sUpperBound += scalingFactor*(oldUpperBound-oldLowerBound);
 
                 mesher = ext::shared_ptr<Fdm1dMesher>(
                     new Concentrating1dMesher(sLowerBound, sUpperBound, xGrid_,

diff --git a/ql/methods/finitedifferences/utilities/riskneutraldensitycalculator.cpp b/ql/methods/finitedifferences/utilities/riskneutraldensitycalculator.cpp
@@ -26,39 +26,25 @@
 namespace QuantLib {
     RiskNeutralDensityCalculator::InvCDFHelper::InvCDFHelper(
         const RiskNeutralDensityCalculator* calculator,
-        Real guess, Real accuracy, Size maxEvaluations)
+        Real guess, Real accuracy, Size maxEvaluations,
+        Real stepSize)
     : calculator_(calculator),
       guess_(guess),
       accuracy_(accuracy),
-      maxEvaluations_(maxEvaluations) { }
+      maxEvaluations_(maxEvaluations),
+      stepSize_(stepSize) { }
 
     Real RiskNeutralDensityCalculator::InvCDFHelper::inverseCDF(Real p, Time t)
     const {
         using namespace ext::placeholders;
 
-        const Real guessCDF = calculator_->cdf(guess_, t);
-
-        Size evaluations = maxEvaluations_;
-        Real upper = guess_, lower = guess_;
-
-        if (guessCDF < p)
-            while (calculator_->cdf(upper*=1.5, t) < p && evaluations > 0) {
-                --evaluations;
-            }
-        else
-            while (calculator_->cdf(lower*=0.75, t) > p && evaluations > 0) {
-                --evaluations;
-            }
-
-        QL_REQUIRE(evaluations, "could not calculate interval");
-
         const ext::function<Real(Real)> cdf
             = ext::bind(&RiskNeutralDensityCalculator::cdf,
                           calculator_, _1, t);
 
         Brent solver;
-        solver.setMaxEvaluations(evaluations);
+        solver.setMaxEvaluations(maxEvaluations_);
         return solver.solve(compose(subtract<Real>(p), cdf),
-                            accuracy_, 0.5*(lower + upper), lower, upper);
+                            accuracy_, guess_, stepSize_);
     }
 }
diff --git a/ql/methods/finitedifferences/utilities/riskneutraldensitycalculator.hpp b/ql/methods/finitedifferences/utilities/riskneutraldensitycalculator.hpp
@@ -40,14 +40,16 @@ namespace QuantLib {
         class InvCDFHelper {
           public:
             InvCDFHelper(const RiskNeutralDensityCalculator* calculator,
-                         Real guess, Real accuracy, Size maxEvaluations);
+                         Real guess, Real accuracy, Size maxEvaluations,
+                         Real stepSize=0.01);
 
             Real inverseCDF(Real p, Time t) const;
           private:
             const RiskNeutralDensityCalculator* const calculator_;
             const Real guess_;
             const Real accuracy_;
             const Size maxEvaluations_;
+            const Real stepSize_;
         };
     };
 }

diff --git a/test-suite/hestonslvmodel.cpp b/test-suite/hestonslvmodel.cpp
@@ -1478,7 +1478,7 @@ namespace {
                         Handle<BlackVolTermStructure>(flatVol(lv,
                                                       dc))));
 
-                const Real tol = 0.0005;//testCase.eps;
+                const Real tol = 0.001;//testCase.eps;
                 if (std::fabs((calculated-expected)/vega) > tol) {
                     BOOST_FAIL("failed to reproduce round trip vola "
                               << "\n   strike         " << strike
@@ -1811,7 +1811,7 @@ void HestonSLVModelTest::testBarrierPricingMixedModels() {
 
     const Real localDeltaExpected =  -0.439068;
     const Real hestonDeltaExpected = -0.342059;
-    const Real tol = 0.0005;
+    const Real tol = 0.001;
     if (std::fabs(hestonDeltaExpected - hestonDeltaCalculated) > tol) {
         BOOST_ERROR("Heston Delta does not match"
                 << "\n calculated : " << hestonDeltaCalculated
@@ -2272,7 +2272,7 @@ void HestonSLVModelTest::testForwardSkewSLV() {
                 QuantLib::detail::ImpliedVolatilityHelper::calculate(
                     *fwdOption, *fwdEngine, *vol, npv, 1e-8, 200, 1e-4, 2.0);
 
-            const Real tol = 0.001;
+            const Real tol = 0.002;
             const Volatility volError = std::fabs(implVol - expected[j]);
 
             if (volError > tol) {