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 <http://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 30 namespace QuantLib { 31 unsafeSabrVolatility(Rate strike,Rate forward,Time expiryTime,Real alpha,Real beta,Real nu,Real rho)32 Real unsafeSabrVolatility(Rate strike, 33 Rate forward, 34 Time expiryTime, 35 Real alpha, 36 Real beta, 37 Real nu, 38 Real rho) { 39 const Real oneMinusBeta = 1.0-beta; 40 const Real A = std::pow(forward*strike, oneMinusBeta); 41 const Real sqrtA= std::sqrt(A); 42 Real logM; 43 if (!close(forward, strike)) 44 logM = std::log(forward/strike); 45 else { 46 const Real epsilon = (forward-strike)/strike; 47 logM = epsilon - .5 * epsilon * epsilon ; 48 } 49 const Real z = (nu/alpha)*sqrtA*logM; 50 const Real B = 1.0-2.0*rho*z+z*z; 51 const Real C = oneMinusBeta*oneMinusBeta*logM*logM; 52 const Real tmp = (std::sqrt(B)+z-rho)/(1.0-rho); 53 const Real xx = std::log(tmp); 54 const Real D = sqrtA*(1.0+C/24.0+C*C/1920.0); 55 const Real d = 1.0 + expiryTime * 56 (oneMinusBeta*oneMinusBeta*alpha*alpha/(24.0*A) 57 + 0.25*rho*beta*nu*alpha/sqrtA 58 +(2.0-3.0*rho*rho)*(nu*nu/24.0)); 59 60 Real multiplier; 61 // computations become precise enough if the square of z worth 62 // slightly more than the precision machine (hence the m) 63 static const Real m = 10; 64 if (std::fabs(z*z)>QL_EPSILON * m) 65 multiplier = z/xx; 66 else { 67 multiplier = 1.0 - 0.5*rho*z - (3.0*rho*rho-2.0)*z*z/12.0; 68 } 69 return (alpha/D)*multiplier*d; 70 } 71 unsafeShiftedSabrVolatility(Rate strike,Rate forward,Time expiryTime,Real alpha,Real beta,Real nu,Real rho,Real shift)72 Real unsafeShiftedSabrVolatility(Rate strike, 73 Rate forward, 74 Time expiryTime, 75 Real alpha, 76 Real beta, 77 Real nu, 78 Real rho, 79 Real shift) { 80 81 return unsafeSabrVolatility(strike+shift,forward+shift,expiryTime, 82 alpha,beta,nu,rho); 83 84 } 85 validateSabrParameters(Real alpha,Real beta,Real nu,Real rho)86 void validateSabrParameters(Real alpha, 87 Real beta, 88 Real nu, 89 Real rho) { 90 QL_REQUIRE(alpha>0.0, "alpha must be positive: " 91 << alpha << " not allowed"); 92 QL_REQUIRE(beta>=0.0 && beta<=1.0, "beta must be in (0.0, 1.0): " 93 << beta << " not allowed"); 94 QL_REQUIRE(nu>=0.0, "nu must be non negative: " 95 << nu << " not allowed"); 96 QL_REQUIRE(rho*rho<1.0, "rho square must be less than one: " 97 << rho << " not allowed"); 98 } 99 sabrVolatility(Rate strike,Rate forward,Time expiryTime,Real alpha,Real beta,Real nu,Real rho)100 Real sabrVolatility(Rate strike, 101 Rate forward, 102 Time expiryTime, 103 Real alpha, 104 Real beta, 105 Real nu, 106 Real rho) { 107 QL_REQUIRE(strike>0.0, "strike must be positive: " 108 << io::rate(strike) << " not allowed"); 109 QL_REQUIRE(forward>0.0, "at the money forward rate must be " 110 "positive: " << io::rate(forward) << " not allowed"); 111 QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: " 112 << expiryTime << " not allowed"); 113 validateSabrParameters(alpha, beta, nu, rho); 114 return unsafeSabrVolatility(strike, forward, expiryTime, 115 alpha, beta, nu, rho); 116 } 117 shiftedSabrVolatility(Rate strike,Rate forward,Time expiryTime,Real alpha,Real beta,Real nu,Real rho,Real shift)118 Real shiftedSabrVolatility(Rate strike, 119 Rate forward, 120 Time expiryTime, 121 Real alpha, 122 Real beta, 123 Real nu, 124 Real rho, 125 Real shift) { 126 QL_REQUIRE(strike + shift > 0.0, "strike+shift must be positive: " 127 << io::rate(strike) << "+" << io::rate(shift) << " not allowed"); 128 QL_REQUIRE(forward + shift > 0.0, "at the money forward rate + shift must be " 129 "positive: " << io::rate(forward) << " " << io::rate(shift) << " not allowed"); 130 QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: " 131 << expiryTime << " not allowed"); 132 validateSabrParameters(alpha, beta, nu, rho); 133 return unsafeShiftedSabrVolatility(strike, forward, expiryTime, 134 alpha, beta, nu, rho,shift); 135 } 136 137 namespace { 138 struct SabrFlochKennedyVolatility { 139 Real F, alpha, beta, nu, rho, t; 140 yQuantLib::__anon85a9cde90111::SabrFlochKennedyVolatility141 Real y(Real k) const { 142 return -1.0/(1.0-beta)*(std::pow(F,1-beta)-std::pow(k,1-beta)); 143 } 144 DintQuantLib::__anon85a9cde90111::SabrFlochKennedyVolatility145 Real Dint(Real k) const { 146 return 1/nu*std::log( ( std::sqrt(1+2*rho*nu/alpha*y(k) 147 + square<Real>()(nu/alpha*y(k)) ) 148 - rho - nu/alpha*y(k) ) / (1-rho) ); 149 } 150 DQuantLib::__anon85a9cde90111::SabrFlochKennedyVolatility151 Real D(Real k) const { 152 return std::sqrt(alpha*alpha+2*alpha*rho*nu*y(k) 153 + square<Real>()(nu*y(k)))*std::pow(k,beta); 154 } 155 omega0QuantLib::__anon85a9cde90111::SabrFlochKennedyVolatility156 Real omega0(Real k) const { 157 return std::log(F/k)/Dint(k); 158 } 159 operator ()QuantLib::__anon85a9cde90111::SabrFlochKennedyVolatility160 Real operator()(Real k) const { 161 const Real m = F/k; 162 if (m > 1.0025 || m < 0.9975) { 163 return omega0(k)*(1+0.25*rho*nu*alpha* 164 (std::pow(k,beta)-std::pow(F,beta))/(k-F)*t) 165 -omega0(k)/square<Real>()(Dint(k))*(std::log( 166 omega0(k)) + 0.5*std::log((F*k/(D(F)*D(k))) ))*t; 167 } 168 else { 169 return taylorExpansion(k); 170 } 171 } 172 taylorExpansionQuantLib::__anon85a9cde90111::SabrFlochKennedyVolatility173 Real taylorExpansion(Real k) const { 174 const Real F2 = F*F; 175 const Real alpha2 = alpha*alpha; 176 const Real rho2 = rho*rho; 177 return 178 (alpha*std::pow(F,-3 + beta)*(alpha2*square<Real>()(-1 + beta)*std::pow(F,2*beta)*t + 6*alpha*beta*nu*std::pow(F,1 + beta)*rho*t + 179 F2*(24 + nu*nu*(2 - 3*rho2)*t)))/24.0 + 180 (3*alpha2*alpha*std::pow(-1 + beta,3)*std::pow(F,3*beta)*t + 181 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) + 182 alpha*std::pow(F,2 + beta)*(24*(-1 + beta) + nu*nu*(2*(-1 + beta) + 3*(1 + beta)*rho2)*t))/(48.*F2*F2) * (k-F) + 183 (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 - 184 30*alpha*nu*std::pow(F,3 + beta)*rho*(24 + nu*nu*(-4*(1 + beta) + 3*(1 + 2*beta)*rho2)*t) + 185 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) + 186 nu*nu*F2*F2*(480 - 720*rho2 + nu*nu*(-64 + 75*rho2*(4 - 3*rho2))*t)))/(2880*alpha) * (k-F)*(k-F); 187 } 188 }; 189 } 190 sabrFlochKennedyVolatility(Rate strike,Rate forward,Time expiryTime,Real alpha,Real beta,Real nu,Real rho)191 Real sabrFlochKennedyVolatility(Rate strike, 192 Rate forward, 193 Time expiryTime, 194 Real alpha, 195 Real beta, 196 Real nu, 197 Real rho) { 198 const SabrFlochKennedyVolatility v = 199 {forward, alpha, beta, nu, rho, expiryTime}; 200 201 return v(strike); 202 } 203 204 } 205