30255006850

Committed 27 Jul 2026 09:41AM UTC coverage: 74.898% (-0.05%) from 74.944%

Build # 30255006850

Build Type

push

github

Committed by

lballabio

Commit Message

Avoid deprecated Ubuntu runner

Coverage Stats

60031 of 80150 relevant lines covered (74.9%)

8637062.36 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.61

/ql/termstructures/volatility/sabr.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2006 Ferdinando Ametrano
 Copyright (C) 2006 Mario Pucci
 Copyright (C) 2006 StatPro Italia srl
 Copyright (C) 2015 Peter Caspers
 Copyright (C) 2019 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/termstructures/volatility/sabr.hpp>
#include <ql/utilities/dataformatters.hpp>
#include <ql/math/comparison.hpp>
#include <ql/math/functional.hpp>
#include <ql/errors.hpp>
#include <ql/termstructures/volatility/volatilitytype.hpp>
#if BOOST_VERSION >= 107800
#include <boost/math/special_functions/sign.hpp>
#include <boost/math/tools/cubic_roots.hpp>
#endif

namespace QuantLib {

    Real unsafeSabrLogNormalVolatility(
                              Rate strike,
                              Rate forward,
                              Time expiryTime,
                              Real alpha,
                              Real beta,
                              Real nu,
                              Real rho) {
        const Real oneMinusBeta = 1.0-beta;
        const Real A = std::pow(forward*strike, oneMinusBeta);
        const Real sqrtA= std::sqrt(A);
        Real logM;
        if (!close(forward, strike))
            logM = std::log(forward/strike);
        else {
            const Real epsilon = (forward-strike)/strike;
            logM = epsilon - .5 * epsilon * epsilon ;
        }
        const Real z = (nu/alpha)*sqrtA*logM;
        const Real B = 1.0-2.0*rho*z+z*z;
        const Real C = oneMinusBeta*oneMinusBeta*logM*logM;
        const Real tmp = (std::sqrt(B)+z-rho)/(1.0-rho);
        const Real xx = std::log(tmp);
        const Real D = sqrtA*(1.0+C/24.0+C*C/1920.0);
        const Real d = 1.0 + expiryTime *
            (oneMinusBeta*oneMinusBeta*alpha*alpha/(24.0*A)
                                + 0.25*rho*beta*nu*alpha/sqrtA
                                    +(2.0-3.0*rho*rho)*(nu*nu/24.0));

        Real multiplier;
        // computations become precise enough if the square of z worth
        // slightly more than the precision machine (hence the m)
        static const Real m = 10;
        if (std::fabs(z*z)>QL_EPSILON * m)
            multiplier = z/xx;
        else {
            multiplier = 1.0 - 0.5*rho*z - (3.0*rho*rho-2.0)*z*z/12.0;
        }
        return (alpha/D)*multiplier*d;
    }

    Real unsafeShiftedSabrVolatility(Rate strike,
                              Rate forward,
                              Time expiryTime,
                              Real alpha,
                              Real beta,
                              Real nu,
                              Real rho,
                              Real shift,
                              VolatilityType volatilityType) {
        if (volatilityType == VolatilityType::Normal) {
            return unsafeSabrNormalVolatility(strike + shift, forward + shift, expiryTime, alpha, beta, nu, rho);
        } else {
            return unsafeSabrLogNormalVolatility(strike + shift, forward + shift, expiryTime, alpha, beta, nu, rho);
        }
    }

    Real unsafeSabrNormalVolatility(
        Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho) {
        const Real oneMinusBeta = 1.0 - beta;
        const Real minusBeta = -1.0 * beta;
        const Real A = std::pow(forward * strike, oneMinusBeta);
        const Real sqrtA = std::sqrt(A);
        Real logM;
        if (!close(forward, strike))
            logM = std::log(forward / strike);
        else {
            const Real epsilon = (forward - strike) / strike;
            logM = epsilon - .5 * epsilon * epsilon;
        }
        const Real z = (nu / alpha) * sqrtA * logM;
        const Real B = 1.0 - 2.0 * rho * z + z * z;
        const Real C = oneMinusBeta * oneMinusBeta * logM * logM;
        const Real D = logM * logM;
        const Real tmp = (std::sqrt(B) + z - rho) / (1.0 - rho);
        const Real xx = std::log(tmp);
        const Real E_1 = (1.0 + D / 24.0 + D * D / 1920.0);
        const Real E_2 = (1.0 + C / 24.0 + C * C / 1920.0);
        const Real E = E_1 / E_2;
        const Real d = 1.0 + expiryTime * (minusBeta * (2 - beta) * alpha * alpha / (24.0 * A) +
                                0.25 * rho * beta * nu * alpha / sqrtA +
                                (2.0 - 3.0 * rho * rho) * (nu * nu / 24.0));

        Real multiplier;
        // computations become precise enough if the square of z worth
        // slightly more than the precision machine (hence the m)
        static const Real m = 10;
        if (std::fabs(z * z) > QL_EPSILON * m)
            multiplier = z / xx;
        else {
            multiplier = 1.0 - 0.5 * rho * z - (3.0 * rho * rho - 2.0) * z * z / 12.0;
        }
        const Real F = alpha * std::pow(forward * strike, beta / 2.0);

       return F * E * multiplier * d;
    }

     Real unsafeSabrVolatility(Rate strike,
                              Rate forward,
                              Time expiryTime,
                              Real alpha,
                              Real beta,
                              Real nu,
                              Real rho,
                              VolatilityType volatilityType) {
        if (volatilityType == VolatilityType::Normal) {
            return unsafeSabrNormalVolatility(strike, forward, expiryTime, alpha, beta, nu, rho);
        } else {
            return unsafeSabrLogNormalVolatility(strike, forward, expiryTime, alpha, beta, nu, rho);
        }
     }

    void validateSabrParameters(Real alpha,
                                Real beta,
                                Real nu,
                                Real rho) {
        QL_REQUIRE(alpha>0.0, "alpha must be positive: "
                              << alpha << " not allowed");
        QL_REQUIRE(beta>=0.0 && beta<=1.0, "beta must be in (0.0, 1.0): "
                                         << beta << " not allowed");
        QL_REQUIRE(nu>=0.0, "nu must be non negative: "
                            << nu << " not allowed");
        QL_REQUIRE(rho*rho<1.0, "rho square must be less than one: "
                                << rho << " not allowed");
    }

    Real sabrVolatility(Rate strike,
                        Rate forward,
                        Time expiryTime,
                        Real alpha,
                        Real beta,
                        Real nu,
                        Real rho,
                        VolatilityType volatilityType) {
        QL_REQUIRE(strike>0.0, "strike must be positive: "
                               << io::rate(strike) << " not allowed");
        QL_REQUIRE(forward>0.0, "at the money forward rate must be "
                   "positive: " << io::rate(forward) << " not allowed");
        QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: "
                                   << expiryTime << " not allowed");
        validateSabrParameters(alpha, beta, nu, rho);
        return unsafeSabrVolatility(strike, forward, expiryTime, alpha, beta, nu, rho,
                                             volatilityType);
    }

    Real shiftedSabrVolatility(Rate strike,
                               Rate forward,
                               Time expiryTime,
                               Real alpha,
                               Real beta,
                               Real nu,
                               Real rho,
                               Real shift,
                               VolatilityType volatilityType) {
        QL_REQUIRE(strike + shift > 0.0, "strike+shift must be positive: "
                   << io::rate(strike) << "+" << io::rate(shift) << " not allowed");
        QL_REQUIRE(forward + shift > 0.0, "at the money forward rate + shift must be "
                   "positive: " << io::rate(forward) << " " << io::rate(shift) << " not allowed");
        QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: "
                                   << expiryTime << " not allowed");
        validateSabrParameters(alpha, beta, nu, rho);
        return unsafeShiftedSabrVolatility(strike, forward, expiryTime,
                                             alpha, beta, nu, rho,shift, volatilityType);
    }

    namespace {
        struct SabrFlochKennedyVolatility {
            Real F, alpha, beta, nu, rho, t;

            Real y(Real k) const {
                return -1.0/(1.0-beta)*(std::pow(F,1-beta)-std::pow(k,1-beta));
            }

            Real Dint(Real k) const {
                return 1/nu*std::log( ( std::sqrt(1+2*rho*nu/alpha*y(k)
                    + squared(nu/alpha*y(k)) )
                    - rho - nu/alpha*y(k) ) / (1-rho) );
            }

            Real D(Real k) const {
                return std::sqrt(alpha*alpha+2*alpha*rho*nu*y(k)
                    + squared(nu*y(k)))*std::pow(k,beta);
            }

            Real omega0(Real k) const {
                return std::log(F/k)/Dint(k);
            }

            Real operator()(Real k) const {
                const Real m = F/k;
                if (m > 1.0025 || m < 0.9975) {
                    return omega0(k)*(1+0.25*rho*nu*alpha*
                       (std::pow(k,beta)-std::pow(F,beta))/(k-F)*t)
                       -omega0(k)/squared(Dint(k))*(std::log(
                           omega0(k)) + 0.5*std::log((F*k/(D(F)*D(k))) ))*t;
                }
                else {
                    return taylorExpansion(k);
                }
            }

            Real taylorExpansion(Real k) const {
                const Real F2 = F*F;
                const Real alpha2 = alpha*alpha;
                const Real rho2 = rho*rho;
                return
                    (alpha*std::pow(F,-3 + beta)*(alpha2*squared(-1 + beta)*std::pow(F,2*beta)*t + 6*alpha*beta*nu*std::pow(F,1 + beta)*rho*t +
                        F2*(24 + nu*nu*(2 - 3*rho2)*t)))/24.0 +
                     (3*alpha2*alpha*std::pow(-1 + beta,3)*std::pow(F,3*beta)*t +
                        3*alpha2*(-1 + beta)*(-1 + 5*beta)*nu*std::pow(F,1 + 2*beta)*rho*t + nu*F2*F*rho*(24 + nu*nu*(-4 + 3*rho2)*t) +
                        alpha*std::pow(F,2 + beta)*(24*(-1 + beta) + nu*nu*(2*(-1 + beta) + 3*(1 + beta)*rho2)*t))/(48.*F2*F2) * (k-F) +
                    (std::pow(F,-5 - beta)*(alpha2*alpha2*std::pow(-1 + beta,3)*(-209 + 119*beta)*std::pow(F,4*beta)*t + 30*alpha2*alpha*(-1 + beta)*(9 + beta*(-37 + 18*beta))*nu*std::pow(F,1 + 3*beta)*rho*t -
                        30*alpha*nu*std::pow(F,3 + beta)*rho*(24 + nu*nu*(-4*(1 + beta) + 3*(1 + 2*beta)*rho2)*t) +
                        10*alpha2*std::pow(F,2 + 2*beta)*(24*(-4 + beta)*(-1 + beta) + nu*nu*(2*(-1 + beta)*(-7 + 4*beta) + 3*(-4 + beta*(-7 + 5*beta))*rho2)*t) +
                        nu*nu*F2*F2*(480 - 720*rho2 + nu*nu*(-64 + 75*rho2*(4 - 3*rho2))*t)))/(2880*alpha) * (k-F)*(k-F);
            }
        };
    }

    Real sabrFlochKennedyVolatility(Rate strike,
                                Rate forward,
                                Time expiryTime,
                                Real alpha,
                                Real beta,
                                Real nu,
                                Real rho) {
        const SabrFlochKennedyVolatility v =
            {forward, alpha, beta, nu, rho, expiryTime};

        return v(strike);
    }


    #if BOOST_VERSION >= 107800

    namespace {

        Real smallest_positive_root(Real c1, Real c2, Real c3, Real c4) {
            auto [r1, r2, r3] = boost::math::tools::cubic_roots(c1, c2, c3, c4);
            if (std::isnan(r3)) {
                // single root (or two equal ones), check that it's positive
                QL_REQUIRE(r1 > 0.0, "no positive root");
                return r1;
            } else {
                // three roots in non-decreasing order, return the first positive one
                QL_REQUIRE(r3 > 0.0, "no positive root");
                return r1 > 0.0 ? r1 : (r2 > 0.0 ? r2 : r3);
            }
        }

    }

    std::array<Real, 4> sabrGuess(Real k_m, Volatility vol_m,
                                  Real k_0, Volatility vol_0,
                                  Real k_p, Volatility vol_p,
                                  Rate forward,
                                  Time expiryTime,
                                  Real beta,
                                  Real shift,
                                  VolatilityType volatilityType) {

        // same variable names as in the equations for ease of reference:
        Real f = forward, b = shift, T = expiryTime;

        // change to log-moneyness

        Real z_m = std::log((k_m + b) / (f + b));
        Real z_0 = std::log((k_0 + b) / (f + b));
        Real z_p = std::log((k_p + b) / (f + b));

        // calculate atm, skew, curvature

        Real w_m = 1 / ((z_m - z_0) * (z_m - z_p));  // eq. (42) in the paper
        Real w_0 = 1 / ((z_0 - z_m) * (z_0 - z_p));  // eq. (43)
        Real w_p = 1 / ((z_p - z_m) * (z_p - z_0));  // eq. (44)

        Real sigma_0 = z_0 * z_p * w_m * vol_m + z_m * z_p * w_0 * vol_0 + z_m * z_0 * w_p * vol_p;         // (39)
        Real sigma_1 = - (z_0 + z_p) * w_m * vol_m - (z_m + z_p) * w_0 * vol_0 - (z_m + z_0) * w_p * vol_p; // (40)
        Real sigma_2 = 2 * w_m * vol_m + 2 * w_0 * vol_0 + 2 * w_p * vol_p;                                 // (41)

        switch (volatilityType) {
          case ShiftedLognormal: {

              // equations (32)

              Real alpha = sigma_0 * std::pow(f + b, 1.0-beta); // NOLINT(clang-analyzer-deadcode.DeadStores)
              Real nu2 =
                  3 * sigma_0 * sigma_2
                  - 0.5 * squared(1-beta) * sigma_0 * sigma_0
                  + 1.5 * squared(2*sigma_1 + (1-beta)*sigma_0);
              Real nu, rho;
              if (nu2 > 0.0) {
                  nu = std::sqrt(nu2);
                  rho = (1/nu) * (2*sigma_1 + (1-beta)*sigma_0);
              } else {
                  rho = boost::math::sign(2*sigma_1 + (1-beta)*sigma_0);
                  nu = (1/rho) * (2*sigma_1 + (1-beta)*sigma_0);
              }

              // coefficients of the polynomial in equation (33)

              Real c1 = squared(1 - beta) * T / (24 * std::pow(f + b, 2 - 2 * beta));
              Real c2 = rho * beta * nu * T / (4 * std::pow(f + b, 1 - beta));
              Real c3 = 1 + ((2 - 3 * rho*rho) / 24) * nu*nu * T;
              Real c4 = - sigma_0 * std::pow(f + b, 1-beta);

              try {
                  alpha = smallest_positive_root(c1, c2, c3, c4);
              } catch (Error&) {}

              return { alpha, beta, nu, rho };
          }
          case Normal: {

              // equations (37)

              Real alpha = sigma_0 * std::pow(f + b, -beta); // NOLINT(clang-analyzer-deadcode.DeadStores)
              Real nu2 = squared(1 / (f + b)) * (
                  3 * sigma_0 * sigma_2
                  - 0.5 * (beta*beta + beta) * (sigma_0*sigma_0)
                  - 3 * sigma_0 * (sigma_1 - 0.5 * beta * sigma_0)
                  + 1.5 * squared(2 * sigma_1 - beta * sigma_0)
              );
              Real nu, rho;
              if (nu2 > 0.0) {
                  nu = std::sqrt(nu2);
                  rho = (1 / (nu * (f + b))) * (2 * sigma_1 - beta * sigma_0);
              } else {
                  rho = boost::math::sign((1 / (f + b)) * (2 * sigma_1 - beta * sigma_0));
                  nu = (1 / (rho * (f + b))) * (2 * sigma_1 - beta * sigma_0);
              }

              // coefficients of the polynomial in equation (38)

              Real c1 = (beta * beta - 2 * beta) * T / (24 * std::pow(f + b, 2 - 2 * beta));
              Real c2 = rho * beta * nu * T / (4 * std::pow(f + b, 1 - beta));
              Real c3 = 1 + ((2 - 3 * rho*rho) / 24) * nu*nu * T;
              Real c4 = - sigma_0 * std::pow(f + b, -beta);

              try {
                  alpha = smallest_positive_root(c1, c2, c3, c4);
              } catch (Error&) {}

              return { alpha, beta, nu, rho };
          }
          default:
            QL_FAIL("unknown volatility type: " << Integer(volatilityType));
        }
    }

    #else

    std::array<Real, 4> sabrGuess(Real, Volatility,
                                  Real, Volatility,
                                  Real, Volatility,
                                  Rate,
                                  Time,
                                  Real,
                                  Real,
                                  VolatilityType) {
        QL_FAIL("Boost 1.78 or later is required for the implementation of this functionality");
    }

    #endif

}

1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2006 Ferdinando Ametrano
5	Copyright (C) 2006 Mario Pucci
6	Copyright (C) 2006 StatPro Italia srl
7	Copyright (C) 2015 Peter Caspers
8	Copyright (C) 2019 Klaus Spanderen
9
10	This file is part of QuantLib, a free-software/open-source library
11	for financial quantitative analysts and developers - http://quantlib.org/
12
13	QuantLib is free software: you can redistribute it and/or modify it
14	under the terms of the QuantLib license. You should have received a
15	copy of the license along with this program; if not, please email
16	<quantlib-dev@lists.sf.net>. The license is also available online at
17	<https://www.quantlib.org/license.shtml>.
18
19	This program is distributed in the hope that it will be useful, but WITHOUT
20	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21	FOR A PARTICULAR PURPOSE. See the license for more details.
22	*/
23
24	#include <ql/termstructures/volatility/sabr.hpp>
25	#include <ql/utilities/dataformatters.hpp>
26	#include <ql/math/comparison.hpp>
27	#include <ql/math/functional.hpp>
28	#include <ql/errors.hpp>
29	#include <ql/termstructures/volatility/volatilitytype.hpp>
30	#if BOOST_VERSION >= 107800
31	#include <boost/math/special_functions/sign.hpp>
32	#include <boost/math/tools/cubic_roots.hpp>
33	#endif
34
35	namespace QuantLib {
36
37	Real unsafeSabrLogNormalVolatility(	5,770,040✔
38	Rate strike,
39	Rate forward,
40	Time expiryTime,
41	Real alpha,
42	Real beta,
43	Real nu,
44	Real rho) {
45	const Real oneMinusBeta = 1.0-beta;	5,770,040✔
46	const Real A = std::pow(forward*strike, oneMinusBeta);	5,770,040✔
47	const Real sqrtA= std::sqrt(A);	5,770,040✔
48	Real logM;
49	if (!close(forward, strike))	5,770,040✔
50	logM = std::log(forward/strike);	5,281,632✔
51	else {
52	const Real epsilon = (forward-strike)/strike;	488,408✔
53	logM = epsilon - .5 * epsilon * epsilon ;	488,408✔
54	}
55	const Real z = (nu/alpha)sqrtAlogM;	5,770,040✔
56	const Real B = 1.0-2.0rhoz+z*z;	5,770,040✔
57	const Real C = oneMinusBetaoneMinusBetalogM*logM;	5,770,040✔
58	const Real tmp = (std::sqrt(B)+z-rho)/(1.0-rho);	5,770,040✔
59	const Real xx = std::log(tmp);	5,770,040✔
60	const Real D = sqrtA(1.0+C/24.0+CC/1920.0);	5,770,040✔
61	const Real d = 1.0 + expiryTime *	5,770,040✔
62	(oneMinusBetaoneMinusBetaalphaalpha/(24.0A)	5,770,040✔
63	+ 0.25rhobetanualpha/sqrtA	5,770,040✔
64	+(2.0-3.0rhorho)(nunu/24.0));	5,770,040✔
65
66	Real multiplier;
67	// computations become precise enough if the square of z worth
68	// slightly more than the precision machine (hence the m)
69	static const Real m = 10;
70	if (std::fabs(zz)>QL_EPSILON m)	5,770,040✔
71	multiplier = z/xx;	5,281,421✔
72	else {
73	multiplier = 1.0 - 0.5rhoz - (3.0rhorho-2.0)zz/12.0;	488,619✔
74	}
75	return (alpha/D)multiplierd;	5,770,040✔
76	}
77
78	Real unsafeShiftedSabrVolatility(Rate strike,	23,432,544✔
79	Rate forward,
80	Time expiryTime,
81	Real alpha,
82	Real beta,
83	Real nu,
84	Real rho,
85	Real shift,
86	VolatilityType volatilityType) {
87	if (volatilityType == VolatilityType::Normal) {	23,432,544✔
88	return unsafeSabrNormalVolatility(strike + shift, forward + shift, expiryTime, alpha, beta, nu, rho);	17,692,173✔
89	} else {
90	return unsafeSabrLogNormalVolatility(strike + shift, forward + shift, expiryTime, alpha, beta, nu, rho);	5,740,371✔
91	}
92	}
93
94	Real unsafeSabrNormalVolatility(	17,692,173✔
95	Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho) {
96	const Real oneMinusBeta = 1.0 - beta;	17,692,173✔
97	const Real minusBeta = -1.0 * beta;	17,692,173✔
98	const Real A = std::pow(forward * strike, oneMinusBeta);	17,692,173✔
99	const Real sqrtA = std::sqrt(A);	17,692,173✔
100	Real logM;
101	if (!close(forward, strike))	17,692,173✔
102	logM = std::log(forward / strike);	14,154,354✔
103	else {
104	const Real epsilon = (forward - strike) / strike;	3,537,819✔
105	logM = epsilon - .5 * epsilon * epsilon;	3,537,819✔
106	}
107	const Real z = (nu / alpha) * sqrtA * logM;	17,692,173✔
108	const Real B = 1.0 - 2.0 * rho * z + z * z;	17,692,173✔
109	const Real C = oneMinusBeta * oneMinusBeta * logM * logM;	17,692,173✔
110	const Real D = logM * logM;	17,692,173✔
111	const Real tmp = (std::sqrt(B) + z - rho) / (1.0 - rho);	17,692,173✔
112	const Real xx = std::log(tmp);	17,692,173✔
113	const Real E_1 = (1.0 + D / 24.0 + D * D / 1920.0);	17,692,173✔
114	const Real E_2 = (1.0 + C / 24.0 + C * C / 1920.0);	17,692,173✔
115	const Real E = E_1 / E_2;	17,692,173✔
116	const Real d = 1.0 + expiryTime * (minusBeta * (2 - beta) * alpha * alpha / (24.0 * A) +	17,692,173✔
117	0.25 * rho * beta * nu * alpha / sqrtA +	17,692,173✔
118	(2.0 - 3.0 * rho * rho) * (nu * nu / 24.0));	17,692,173✔
119
120	Real multiplier;
121	// computations become precise enough if the square of z worth
122	// slightly more than the precision machine (hence the m)
123	static const Real m = 10;
124	if (std::fabs(z * z) > QL_EPSILON * m)	17,692,173✔
125	multiplier = z / xx;	14,150,904✔
126	else {
127	multiplier = 1.0 - 0.5 * rho * z - (3.0 * rho * rho - 2.0) * z * z / 12.0;	3,541,269✔
128	}
129	const Real F = alpha * std::pow(forward * strike, beta / 2.0);	17,692,173✔
130
131	return F * E * multiplier * d;	17,692,173✔
132	}
133
134	Real unsafeSabrVolatility(Rate strike,	29,669✔
135	Rate forward,
136	Time expiryTime,
137	Real alpha,
138	Real beta,
139	Real nu,
140	Real rho,
141	VolatilityType volatilityType) {
142	if (volatilityType == VolatilityType::Normal) {	29,669✔
143	return unsafeSabrNormalVolatility(strike, forward, expiryTime, alpha, beta, nu, rho);	×
144	} else {
145	return unsafeSabrLogNormalVolatility(strike, forward, expiryTime, alpha, beta, nu, rho);	29,669✔
146	}
147	}
148
149	void validateSabrParameters(Real alpha,	27,233,752✔
150	Real beta,
151	Real nu,
152	Real rho) {
153	QL_REQUIRE(alpha>0.0, "alpha must be positive: "	27,233,752✔
154	<< alpha << " not allowed");
155	QL_REQUIRE(beta>=0.0 && beta<=1.0, "beta must be in (0.0, 1.0): "	27,233,752✔
156	<< beta << " not allowed");
157	QL_REQUIRE(nu>=0.0, "nu must be non negative: "	27,233,752✔
158	<< nu << " not allowed");
159	QL_REQUIRE(rho*rho<1.0, "rho square must be less than one: "	27,233,752✔
160	<< rho << " not allowed");
161	}	27,233,752✔
162
163	Real sabrVolatility(Rate strike,	29,533✔
164	Rate forward,
165	Time expiryTime,
166	Real alpha,
167	Real beta,
168	Real nu,
169	Real rho,
170	VolatilityType volatilityType) {
171	QL_REQUIRE(strike>0.0, "strike must be positive: "	29,533✔
172	<< io::rate(strike) << " not allowed");
173	QL_REQUIRE(forward>0.0, "at the money forward rate must be "	29,533✔
174	"positive: " << io::rate(forward) << " not allowed");
175	QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: "	29,533✔
176	<< expiryTime << " not allowed");
177	validateSabrParameters(alpha, beta, nu, rho);	29,533✔
178	return unsafeSabrVolatility(strike, forward, expiryTime, alpha, beta, nu, rho,	29,533✔
179	volatilityType);	29,533✔
180	}
181
182	Real shiftedSabrVolatility(Rate strike,	22,102,621✔
183	Rate forward,
184	Time expiryTime,
185	Real alpha,
186	Real beta,
187	Real nu,
188	Real rho,
189	Real shift,
190	VolatilityType volatilityType) {
191	QL_REQUIRE(strike + shift > 0.0, "strike+shift must be positive: "	22,102,621✔
192	<< io::rate(strike) << "+" << io::rate(shift) << " not allowed");
193	QL_REQUIRE(forward + shift > 0.0, "at the money forward rate + shift must be "	22,102,621✔
194	"positive: " << io::rate(forward) << " " << io::rate(shift) << " not allowed");
195	QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: "	22,102,621✔
196	<< expiryTime << " not allowed");
197	validateSabrParameters(alpha, beta, nu, rho);	22,102,621✔
198	return unsafeShiftedSabrVolatility(strike, forward, expiryTime,	22,102,621✔
199	alpha, beta, nu, rho,shift, volatilityType);	22,102,621✔
200	}
201
202	namespace {
203	struct SabrFlochKennedyVolatility {
204	Real F, alpha, beta, nu, rho, t;
205
206	Real y(Real k) const {	528✔
207	return -1.0/(1.0-beta)*(std::pow(F,1-beta)-std::pow(k,1-beta));	528✔
208	}
209
210	Real Dint(Real k) const {	132✔
211	return 1/nustd::log( ( std::sqrt(1+2rhonu/alphay(k)	132✔
212	+ squared(nu/alpha*y(k)) )	132✔
213	- rho - nu/alpha*y(k) ) / (1-rho) );	132✔
214	}
215
216	Real D(Real k) const {	66✔
217	return std::sqrt(alphaalpha+2alpharhonu*y(k)	66✔
218	+ squared(nuy(k)))std::pow(k,beta);	66✔
219	}
220
221	Real omega0(Real k) const {	99✔
222	return std::log(F/k)/Dint(k);	99✔
223	}
224
225	Real operator()(Real k) const {	85✔
226	const Real m = F/k;	85✔
227	if (m > 1.0025 \|\| m < 0.9975) {	85✔
228	return omega0(k)(1+0.25rhonualpha*	33✔
229	(std::pow(k,beta)-std::pow(F,beta))/(k-F)*t)	33✔
230	-omega0(k)/squared(Dint(k))*(std::log(	33✔
231	omega0(k)) + 0.5std::log((Fk/(D(F)D(k))) ))t;	33✔
232	}
233	else {
234	return taylorExpansion(k);	52✔
235	}
236	}
237
238	Real taylorExpansion(Real k) const {	52✔
239	const Real F2 = F*F;	52✔
240	const Real alpha2 = alpha*alpha;	52✔
241	const Real rho2 = rho*rho;	52✔
242	return
243	(alphastd::pow(F,-3 + beta)(alpha2squared(-1 + beta)std::pow(F,2beta)t + 6alphabetanustd::pow(F,1 + beta)rhot +	52✔
244	F2(24 + nunu(2 - 3rho2)*t)))/24.0 +	52✔
245	(3alpha2alphastd::pow(-1 + beta,3)std::pow(F,3beta)t +	52✔
246	3alpha2(-1 + beta)(-1 + 5beta)nustd::pow(F,1 + 2beta)rhot + nuF2Frho(24 + nunu(-4 + 3rho2)*t) +	52✔
247	alphastd::pow(F,2 + beta)(24(-1 + beta) + nunu(2(-1 + beta) + 3(1 + beta)rho2)t))/(48.F2F2) (k-F) +	52✔
248	(std::pow(F,-5 - beta)(alpha2alpha2std::pow(-1 + beta,3)(-209 + 119beta)std::pow(F,4beta)t + 30alpha2alpha(-1 + beta)(9 + beta(-37 + 18beta))nustd::pow(F,1 + 3beta)rho*t -	52✔
249	30alphanustd::pow(F,3 + beta)rho(24 + nunu(-4(1 + beta) + 3(1 + 2beta)rho2)t) +	52✔
250	10alpha2std::pow(F,2 + 2beta)(24(-4 + beta)(-1 + beta) + nunu(2(-1 + beta)(-7 + 4beta) + 3(-4 + beta(-7 + 5beta))rho2)t) +	52✔
251	nunuF2F2(480 - 720rho2 + nunu(-64 + 75rho2(4 - 3rho2))t)))/(2880alpha) * (k-F)*(k-F);	52✔
252	}
253	};
254	}
255
256	Real sabrFlochKennedyVolatility(Rate strike,	85✔
257	Rate forward,
258	Time expiryTime,
259	Real alpha,
260	Real beta,
261	Real nu,
262	Real rho) {
263	const SabrFlochKennedyVolatility v =
264	{forward, alpha, beta, nu, rho, expiryTime};	85✔
265
266	return v(strike);	85✔
267	}
268
269
270	#if BOOST_VERSION >= 107800
271
272	namespace {
273
274	Real smallest_positive_root(Real c1, Real c2, Real c3, Real c4) {	66✔
275	auto [r1, r2, r3] = boost::math::tools::cubic_roots(c1, c2, c3, c4);	66✔
276	if (std::isnan(r3)) {	66✔
277	// single root (or two equal ones), check that it's positive
278	QL_REQUIRE(r1 > 0.0, "no positive root");	33✔
279	return r1;
280	} else {
281	// three roots in non-decreasing order, return the first positive one
282	QL_REQUIRE(r3 > 0.0, "no positive root");	33✔
283	return r1 > 0.0 ? r1 : (r2 > 0.0 ? r2 : r3);	33✔
284	}
285	}
286
287	}
288
289	std::array<Real, 4> sabrGuess(Real k_m, Volatility vol_m,	66✔
290	Real k_0, Volatility vol_0,
291	Real k_p, Volatility vol_p,
292	Rate forward,
293	Time expiryTime,
294	Real beta,
295	Real shift,
296	VolatilityType volatilityType) {
297
298	// same variable names as in the equations for ease of reference:
299	Real f = forward, b = shift, T = expiryTime;
300
301	// change to log-moneyness
302
303	Real z_m = std::log((k_m + b) / (f + b));	66✔
304	Real z_0 = std::log((k_0 + b) / (f + b));	66✔
305	Real z_p = std::log((k_p + b) / (f + b));	66✔
306
307	// calculate atm, skew, curvature
308
309	Real w_m = 1 / ((z_m - z_0) * (z_m - z_p)); // eq. (42) in the paper	66✔
310	Real w_0 = 1 / ((z_0 - z_m) * (z_0 - z_p)); // eq. (43)	66✔
311	Real w_p = 1 / ((z_p - z_m) * (z_p - z_0)); // eq. (44)	66✔
312
313	Real sigma_0 = z_0 * z_p * w_m * vol_m + z_m * z_p * w_0 * vol_0 + z_m * z_0 * w_p * vol_p; // (39)	66✔
314	Real sigma_1 = - (z_0 + z_p) * w_m * vol_m - (z_m + z_p) * w_0 * vol_0 - (z_m + z_0) * w_p * vol_p; // (40)	66✔
315	Real sigma_2 = 2 * w_m * vol_m + 2 * w_0 * vol_0 + 2 * w_p * vol_p; // (41)	66✔
316
317	switch (volatilityType) {	66✔
318	case ShiftedLognormal: {	33✔
319
320	// equations (32)
321
322	Real alpha = sigma_0 * std::pow(f + b, 1.0-beta); // NOLINT(clang-analyzer-deadcode.DeadStores)	33✔
323	Real nu2 =
324	3 * sigma_0 * sigma_2	33✔
325	- 0.5 * squared(1-beta) * sigma_0 * sigma_0	33✔
326	+ 1.5 * squared(2sigma_1 + (1-beta)sigma_0);	33✔
327	Real nu, rho;
328	if (nu2 > 0.0) {	33✔
329	nu = std::sqrt(nu2);	33✔
330	rho = (1/nu) * (2sigma_1 + (1-beta)sigma_0);	33✔
331	} else {
332	rho = boost::math::sign(2sigma_1 + (1-beta)sigma_0);	×
333	nu = (1/rho) * (2sigma_1 + (1-beta)sigma_0);	×
334	}
335
336	// coefficients of the polynomial in equation (33)
337
338	Real c1 = squared(1 - beta) * T / (24 * std::pow(f + b, 2 - 2 * beta));	33✔
339	Real c2 = rho * beta * nu * T / (4 * std::pow(f + b, 1 - beta));	33✔
340	Real c3 = 1 + ((2 - 3 * rhorho) / 24) nunu T;	33✔
341	Real c4 = - sigma_0 * std::pow(f + b, 1-beta);	33✔
342
343	try {
344	alpha = smallest_positive_root(c1, c2, c3, c4);	33✔
345	} catch (Error&) {}	×
346
347	return { alpha, beta, nu, rho };	33✔
348	}
349	case Normal: {	33✔
350
351	// equations (37)
352
353	Real alpha = sigma_0 * std::pow(f + b, -beta); // NOLINT(clang-analyzer-deadcode.DeadStores)	33✔
354	Real nu2 = squared(1 / (f + b)) * (	33✔
355	3 * sigma_0 * sigma_2	33✔
356	- 0.5 * (betabeta + beta) (sigma_0*sigma_0)	33✔
357	- 3 * sigma_0 * (sigma_1 - 0.5 * beta * sigma_0)	33✔
358	+ 1.5 * squared(2 * sigma_1 - beta * sigma_0)	33✔
359	);	33✔
360	Real nu, rho;
361	if (nu2 > 0.0) {	33✔
362	nu = std::sqrt(nu2);	33✔
363	rho = (1 / (nu * (f + b))) * (2 * sigma_1 - beta * sigma_0);	33✔
364	} else {
365	rho = boost::math::sign((1 / (f + b)) * (2 * sigma_1 - beta * sigma_0));	×
366	nu = (1 / (rho * (f + b))) * (2 * sigma_1 - beta * sigma_0);	×
367	}
368
369	// coefficients of the polynomial in equation (38)
370
371	Real c1 = (beta * beta - 2 * beta) * T / (24 * std::pow(f + b, 2 - 2 * beta));	33✔
372	Real c2 = rho * beta * nu * T / (4 * std::pow(f + b, 1 - beta));	33✔
373	Real c3 = 1 + ((2 - 3 * rhorho) / 24) nunu T;	33✔
374	Real c4 = - sigma_0 * std::pow(f + b, -beta);	33✔
375
376	try {
377	alpha = smallest_positive_root(c1, c2, c3, c4);	33✔
378	} catch (Error&) {}	×
379
380	return { alpha, beta, nu, rho };	33✔
381	}
382	default:	×
383	QL_FAIL("unknown volatility type: " << Integer(volatilityType));	×
384	}
385	}
386
387	#else
388
389	std::array<Real, 4> sabrGuess(Real, Volatility,
390	Real, Volatility,
391	Real, Volatility,
392	Rate,
393	Time,
394	Real,
395	Real,
396	VolatilityType) {
397	QL_FAIL("Boost 1.78 or later is required for the implementation of this functionality");
398	}
399
400	#endif
401
402	}

lballabio / QuantLib / 30255006850

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous