diff --git a/ql/pricingengines/blackformula.cpp b/ql/pricingengines/blackformula.cpp index 2cf860d693a..6784483b80f 100644 --- 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 #include #include +#include #include #include #include @@ -833,8 +834,89 @@ namespace QuantLib { return impliedBpvol; } + namespace { + static boost::math::normal_distribution 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 index b85f60551c4..21769973402 100644 --- 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 index 99360123205..38ab61bfdae 100644 --- 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 ); + } + } }