add exact Bachelier implied vol from Jaeckel paper

lballabio · Jul 19, 2024 · 41f235e · 41f235e
1 parent 8e1f4ac
commit 41f235e
Show file tree

Hide file tree

Showing 3 changed files with 106 additions and 6 deletions.
diff --git a/ql/pricingengines/blackformula.cpp b/ql/pricingengines/blackformula.cpp
@@ -8,7 +8,7 @@
  Copyright (C) 2007 Cristina Duminuco
  Copyright (C) 2007 Chiara Fornarola
  Copyright (C) 2013 Gary Kennedy
- Copyright (C) 2015 Peter Caspers
+ Copyright (C) 2015, 2024 Peter Caspers
  Copyright (C) 2017 Klaus Spanderen
  Copyright (C) 2019 Wojciech Ślusarski
  Copyright (C) 2020 Marcin Rybacki
@@ -31,6 +31,7 @@
 #include <ql/math/functional.hpp>
 #include <ql/math/solvers1d/newtonsafe.hpp>
 #include <ql/math/distributions/normaldistribution.hpp>
+#include <boost/math/distributions/normal.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
 #include <boost/math/special_functions/atanh.hpp>
 #include <boost/math/special_functions/sign.hpp>
@@ -833,8 +834,89 @@ namespace QuantLib {
         return impliedBpvol;
     }
 
+    namespace {
+        static boost::math::normal_distribution<double> normal_dist;
+        Real phi(const Real x) {
+            return boost::math::pdf(normal_dist, x);
+        }
+        Real Phi(const Real x) {
+            return boost::math::cdf(normal_dist, x);
+        }
+
+        Real PhiTilde(const Real x) {
+            return Phi(x) + phi(x) / x;
+        }
+
+        Real inversePhiTilde(const Real PhiTildeStar) {
+            QL_REQUIRE(PhiTildeStar < 0.0,
+                       "inversePhiTilde(" << PhiTildeStar << "): negative argument required");
+            Real xbar;
+            if (PhiTildeStar < -0.001882039271) {
+                Real g = 1.0 / (PhiTildeStar - 0.5);
+                Real xibar =
+                    (0.032114372355 -
+                     g * g *
+                         (0.016969777977 - g * g * (2.6207332461E-3 - 9.6066952861E-5 * g * g))) /
+                    (1.0 -
+                     g * g * (0.6635646938 - g * g * (0.14528712196 - 0.010472855461 * g * g)));
+                xbar = g * (0.3989422804014326 + xibar * g * g);
+            } else {
+                Real h = std::sqrt(-std::log(-PhiTildeStar));
+                xbar =
+                    (9.4883409779 - h * (9.6320903635 - h * (0.58556997323 + 2.1464093351 * h))) /
+                    (1.0 - h * (0.65174820867 + h * (1.5120247828 + 6.6437847132E-5 * h)));
+            }
+            Real q = (PhiTilde(xbar) - PhiTildeStar) / phi(xbar);
+            Real xstar =
+                xbar +
+                3.0 * q * xbar * xbar * (2.0 - q * xbar * (2.0 + xbar * xbar)) /
+                    (6.0 + q * xbar *
+                               (-12.0 +
+                                xbar * (6.0 * q + xbar * (-6.0 + q * xbar * (3.0 + xbar * xbar)))));
+            return xstar;
+        }
+    } // namespace
+
+    Real bachelierBlackFormulaImpliedVolExact(Option::Type optionType,
+                                              Real strike,
+                                              Real forward,
+                                              Real tte,
+                                              Real bachelierPrice,
+                                              Real discount) {
+
+        Real theta = optionType == Option::Call ? 1.0 : -1.0;
+
+        // compound bechelierPrice, so that effectively discount = 1
+
+        bachelierPrice /= discount;
+
+        // handle case strike = forward
+
+        if (std::abs(strike - forward) < 1E-15) {
+            return bachelierPrice / (std::sqrt(tte) * phi(0.0));
+        }
+
+        // handle case strike != forward
+
+        Real timeValue = bachelierPrice - std::max(theta * (forward - strike), 0.0);
+
+        if (std::abs(timeValue) < 1E-15)
+            return 0.0;
+
+        QL_REQUIRE(timeValue > 0.0, "bachelierBlackFormulaImpliedVolExact(theta="
+                                        << theta << ",strike=" << strike << ",forward=" << forward
+                                        << ",tte=" << tte << ",price=" << bachelierPrice
+                                        << "): option price implies negative time value ("
+                                        << timeValue << ")");
+
+        Real PhiTildeStar = -std::abs(timeValue / (strike - forward));
+        Real xstar = inversePhiTilde(PhiTildeStar);
+        Real impliedVol = std::abs((strike - forward) / (xstar * std::sqrt(tte)));
+
+        return impliedVol;
+    }
 
-        Real bachelierBlackFormulaStdDevDerivative(Rate strike,
+    Real bachelierBlackFormulaStdDevDerivative(Rate strike,
                                       Rate forward,
                                       Real stdDev,
                                       Real discount)

diff --git a/ql/pricingengines/blackformula.hpp b/ql/pricingengines/blackformula.hpp
@@ -8,7 +8,7 @@
  Copyright (C) 2007 Cristina Duminuco
  Copyright (C) 2007 Chiara Fornarola
  Copyright (C) 2013 Gary Kennedy
- Copyright (C) 2015 Peter Caspers
+ Copyright (C) 2015, 2024 Peter Caspers
  Copyright (C) 2017 Klaus Spanderen
  Copyright (C) 2019 Wojciech Ślusarski
  Copyright (C) 2020 Marcin Rybacki
@@ -372,6 +372,18 @@ namespace QuantLib {
                                          Real bachelierPrice,
                                          Real discount = 1.0);
 
+    /*! Exact Bachelier implied volatility
+
+        It is calculated using the analytic implied volatility formula
+        of Jaeckel (2017), "Implied Normal Volatility"
+    */
+    Real bachelierBlackFormulaImpliedVolExact(Option::Type optionType,
+                                              Real strike,
+                                              Real forward,
+                                              Real tte,
+                                              Real bachelierPrice,
+                                              Real discount = 1.0);
+
     /*! Bachelier formula for standard deviation derivative
         \warning instead of volatility it uses standard deviation, i.e.
                  volatility*sqrt(timeToMaturity), and it returns the

diff --git a/test-suite/blackformula.cpp b/test-suite/blackformula.cpp
@@ -2,7 +2,7 @@
 
 /*
  Copyright (C) 2013 Gary Kennedy
- Copyright (C) 2015 Peter Caspers
+ Copyright (C) 2015, 2024 Peter Caspers
  Copyright (C) 2017 Klaus Spanderen
  Copyright (C) 2020 Marcin Rybacki
 
@@ -35,7 +35,6 @@ BOOST_AUTO_TEST_SUITE(BlackFormulaTests)
 
 BOOST_AUTO_TEST_CASE(testBachelierImpliedVol){
 
-
     BOOST_TEST_MESSAGE("Testing Bachelier implied vol...");
 
     Real forward = 1.0;
@@ -55,9 +54,16 @@ BOOST_AUTO_TEST_CASE(testBachelierImpliedVol){
 
         Real impliedBpVol = bachelierBlackFormulaImpliedVol(optionType,strike, forward, tte, callPrem, discount);
 
-        if (std::fabs(bpvol-impliedBpVol)>1.0e-12){
+        if (std::fabs(bpvol - impliedBpVol) > 1.0e-12) {
             BOOST_ERROR("Failed, expected " << bpvol << " realised " << impliedBpVol );
         }
+
+        Real impliedBpVolExact = bachelierBlackFormulaImpliedVolExact(optionType,strike, forward, tte, callPrem, discount);
+
+        if (std::fabs(bpvol - impliedBpVolExact) > 1.0e-15) {
+            BOOST_ERROR("Failed, expected " << bpvol << " realised " << impliedBpVolExact );
+        }
+
     }
 }