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#include <ql/termstructures/volatility/volatilitytype.hpp>
30
31namespace QuantLib {
32
33 Real unsafeSabrLogNormalVolatility(
34 Rate strike,
35 Rate forward,
36 Time expiryTime,
37 Real alpha,
38 Real beta,
39 Real nu,
40 Real rho) {
41 const Real oneMinusBeta = 1.0-beta;
42 const Real A = std::pow(x: forward*strike, y: oneMinusBeta);
43 const Real sqrtA= std::sqrt(x: A);
44 Real logM;
45 if (!close(x: forward, y: strike))
46 logM = std::log(x: forward/strike);
47 else {
48 const Real epsilon = (forward-strike)/strike;
49 logM = epsilon - .5 * epsilon * epsilon ;
50 }
51 const Real z = (nu/alpha)*sqrtA*logM;
52 const Real B = 1.0-2.0*rho*z+z*z;
53 const Real C = oneMinusBeta*oneMinusBeta*logM*logM;
54 const Real tmp = (std::sqrt(x: B)+z-rho)/(1.0-rho);
55 const Real xx = std::log(x: tmp);
56 const Real D = sqrtA*(1.0+C/24.0+C*C/1920.0);
57 const Real d = 1.0 + expiryTime *
58 (oneMinusBeta*oneMinusBeta*alpha*alpha/(24.0*A)
59 + 0.25*rho*beta*nu*alpha/sqrtA
60 +(2.0-3.0*rho*rho)*(nu*nu/24.0));
61
62 Real multiplier;
63 // computations become precise enough if the square of z worth
64 // slightly more than the precision machine (hence the m)
65 static const Real m = 10;
66 if (std::fabs(x: z*z)>QL_EPSILON * m)
67 multiplier = z/xx;
68 else {
69 multiplier = 1.0 - 0.5*rho*z - (3.0*rho*rho-2.0)*z*z/12.0;
70 }
71 return (alpha/D)*multiplier*d;
72 }
73
74 Real unsafeShiftedSabrVolatility(Rate strike,
75 Rate forward,
76 Time expiryTime,
77 Real alpha,
78 Real beta,
79 Real nu,
80 Real rho,
81 Real shift,
82 VolatilityType volatilityType) {
83 if (volatilityType == VolatilityType::Normal) {
84 return unsafeSabrNormalVolatility(strike: strike + shift, forward: forward + shift, expiryTime, alpha, beta, nu, rho);
85 } else {
86 return unsafeSabrLogNormalVolatility(strike: strike + shift, forward: forward + shift, expiryTime, alpha, beta, nu, rho);
87 }
88 }
89
90 Real unsafeSabrNormalVolatility(
91 Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho) {
92 const Real oneMinusBeta = 1.0 - beta;
93 const Real minusBeta = -1.0 * beta;
94 const Real A = std::pow(x: forward * strike, y: oneMinusBeta);
95 const Real sqrtA = std::sqrt(x: A);
96 Real logM;
97 if (!close(x: forward, y: strike))
98 logM = std::log(x: forward / strike);
99 else {
100 const Real epsilon = (forward - strike) / strike;
101 logM = epsilon - .5 * epsilon * epsilon;
102 }
103 const Real z = (nu / alpha) * sqrtA * logM;
104 const Real B = 1.0 - 2.0 * rho * z + z * z;
105 const Real C = oneMinusBeta * oneMinusBeta * logM * logM;
106 const Real D = logM * logM;
107 const Real tmp = (std::sqrt(x: B) + z - rho) / (1.0 - rho);
108 const Real xx = std::log(x: tmp);
109 const Real E_1 = (1.0 + D / 24.0 + D * D / 1920.0);
110 const Real E_2 = (1.0 + C / 24.0 + C * C / 1920.0);
111 const Real E = E_1 / E_2;
112 const Real d = 1.0 + expiryTime * (minusBeta * (2 - beta) * alpha * alpha / (24.0 * A) +
113 0.25 * rho * beta * nu * alpha / sqrtA +
114 (2.0 - 3.0 * rho * rho) * (nu * nu / 24.0));
115
116 Real multiplier;
117 // computations become precise enough if the square of z worth
118 // slightly more than the precision machine (hence the m)
119 static const Real m = 10;
120 if (std::fabs(x: z * z) > QL_EPSILON * m)
121 multiplier = z / xx;
122 else {
123 multiplier = 1.0 - 0.5 * rho * z - (3.0 * rho * rho - 2.0) * z * z / 12.0;
124 }
125 const Real F = alpha * std::pow(x: forward * strike, y: beta / 2.0);
126
127 return F * E * multiplier * d;
128 }
129
130 Real unsafeSabrVolatility(Rate strike,
131 Rate forward,
132 Time expiryTime,
133 Real alpha,
134 Real beta,
135 Real nu,
136 Real rho,
137 VolatilityType volatilityType) {
138 if (volatilityType == VolatilityType::Normal) {
139 return unsafeSabrNormalVolatility(strike, forward, expiryTime, alpha, beta, nu, rho);
140 } else {
141 return unsafeSabrLogNormalVolatility(strike, forward, expiryTime, alpha, beta, nu, rho);
142 }
143 }
144
145 void validateSabrParameters(Real alpha,
146 Real beta,
147 Real nu,
148 Real rho) {
149 QL_REQUIRE(alpha>0.0, "alpha must be positive: "
150 << alpha << " not allowed");
151 QL_REQUIRE(beta>=0.0 && beta<=1.0, "beta must be in (0.0, 1.0): "
152 << beta << " not allowed");
153 QL_REQUIRE(nu>=0.0, "nu must be non negative: "
154 << nu << " not allowed");
155 QL_REQUIRE(rho*rho<1.0, "rho square must be less than one: "
156 << rho << " not allowed");
157 }
158
159 Real sabrVolatility(Rate strike,
160 Rate forward,
161 Time expiryTime,
162 Real alpha,
163 Real beta,
164 Real nu,
165 Real rho,
166 VolatilityType volatilityType) {
167 QL_REQUIRE(strike>0.0, "strike must be positive: "
168 << io::rate(strike) << " not allowed");
169 QL_REQUIRE(forward>0.0, "at the money forward rate must be "
170 "positive: " << io::rate(forward) << " not allowed");
171 QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: "
172 << expiryTime << " not allowed");
173 validateSabrParameters(alpha, beta, nu, rho);
174 return unsafeSabrVolatility(strike, forward, expiryTime, alpha, beta, nu, rho,
175 volatilityType);
176 }
177
178 Real shiftedSabrVolatility(Rate strike,
179 Rate forward,
180 Time expiryTime,
181 Real alpha,
182 Real beta,
183 Real nu,
184 Real rho,
185 Real shift,
186 VolatilityType volatilityType) {
187 QL_REQUIRE(strike + shift > 0.0, "strike+shift must be positive: "
188 << io::rate(strike) << "+" << io::rate(shift) << " not allowed");
189 QL_REQUIRE(forward + shift > 0.0, "at the money forward rate + shift must be "
190 "positive: " << io::rate(forward) << " " << io::rate(shift) << " not allowed");
191 QL_REQUIRE(expiryTime>=0.0, "expiry time must be non-negative: "
192 << expiryTime << " not allowed");
193 validateSabrParameters(alpha, beta, nu, rho);
194 return unsafeShiftedSabrVolatility(strike, forward, expiryTime,
195 alpha, beta, nu, rho,shift, volatilityType);
196 }
197
198 namespace {
199 struct SabrFlochKennedyVolatility {
200 Real F, alpha, beta, nu, rho, t;
201
202 Real y(Real k) const {
203 return -1.0/(1.0-beta)*(std::pow(x: F,y: 1-beta)-std::pow(x: k,y: 1-beta));
204 }
205
206 Real Dint(Real k) const {
207 return 1/nu*std::log( x: ( std::sqrt(x: 1+2*rho*nu/alpha*y(k)
208 + squared(x: nu/alpha*y(k)) )
209 - rho - nu/alpha*y(k) ) / (1-rho) );
210 }
211
212 Real D(Real k) const {
213 return std::sqrt(x: alpha*alpha+2*alpha*rho*nu*y(k)
214 + squared(x: nu*y(k)))*std::pow(x: k,y: beta);
215 }
216
217 Real omega0(Real k) const {
218 return std::log(x: F/k)/Dint(k);
219 }
220
221 Real operator()(Real k) const {
222 const Real m = F/k;
223 if (m > 1.0025 || m < 0.9975) {
224 return omega0(k)*(1+0.25*rho*nu*alpha*
225 (std::pow(x: k,y: beta)-std::pow(x: F,y: beta))/(k-F)*t)
226 -omega0(k)/squared(x: Dint(k))*(std::log(
227 x: omega0(k)) + 0.5*std::log(x: (F*k/(D(k: F)*D(k))) ))*t;
228 }
229 else {
230 return taylorExpansion(k);
231 }
232 }
233
234 Real taylorExpansion(Real k) const {
235 const Real F2 = F*F;
236 const Real alpha2 = alpha*alpha;
237 const Real rho2 = rho*rho;
238 return
239 (alpha*std::pow(x: F,y: -3 + beta)*(alpha2*squared(x: -1 + beta)*std::pow(x: F,y: 2*beta)*t + 6*alpha*beta*nu*std::pow(x: F,y: 1 + beta)*rho*t +
240 F2*(24 + nu*nu*(2 - 3*rho2)*t)))/24.0 +
241 (3*alpha2*alpha*std::pow(x: -1 + beta,y: 3)*std::pow(x: F,y: 3*beta)*t +
242 3*alpha2*(-1 + beta)*(-1 + 5*beta)*nu*std::pow(x: F,y: 1 + 2*beta)*rho*t + nu*F2*F*rho*(24 + nu*nu*(-4 + 3*rho2)*t) +
243 alpha*std::pow(x: F,y: 2 + beta)*(24*(-1 + beta) + nu*nu*(2*(-1 + beta) + 3*(1 + beta)*rho2)*t))/(48.*F2*F2) * (k-F) +
244 (std::pow(x: F,y: -5 - beta)*(alpha2*alpha2*std::pow(x: -1 + beta,y: 3)*(-209 + 119*beta)*std::pow(x: F,y: 4*beta)*t + 30*alpha2*alpha*(-1 + beta)*(9 + beta*(-37 + 18*beta))*nu*std::pow(x: F,y: 1 + 3*beta)*rho*t -
245 30*alpha*nu*std::pow(x: F,y: 3 + beta)*rho*(24 + nu*nu*(-4*(1 + beta) + 3*(1 + 2*beta)*rho2)*t) +
246 10*alpha2*std::pow(x: F,y: 2 + 2*beta)*(24*(-4 + beta)*(-1 + beta) + nu*nu*(2*(-1 + beta)*(-7 + 4*beta) + 3*(-4 + beta*(-7 + 5*beta))*rho2)*t) +
247 nu*nu*F2*F2*(480 - 720*rho2 + nu*nu*(-64 + 75*rho2*(4 - 3*rho2))*t)))/(2880*alpha) * (k-F)*(k-F);
248 }
249 };
250 }
251
252 Real sabrFlochKennedyVolatility(Rate strike,
253 Rate forward,
254 Time expiryTime,
255 Real alpha,
256 Real beta,
257 Real nu,
258 Real rho) {
259 const SabrFlochKennedyVolatility v =
260 {.F: forward, .alpha: alpha, .beta: beta, .nu: nu, .rho: rho, .t: expiryTime};
261
262 return v(strike);
263 }
264
265}
266

source code of quantlib/ql/termstructures/volatility/sabr.cpp