QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
sabr.cpp
Go to the documentation of this file.
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
28#include <ql/errors.hpp>
30
31namespace QuantLib {
32
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(forward*strike, oneMinusBeta);
43 const Real sqrtA= std::sqrt(A);
44 Real logM;
45 if (!close(forward, strike))
46 logM = std::log(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(B)+z-rho)/(1.0-rho);
55 const Real xx = std::log(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(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
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 + shift, forward + shift, expiryTime, alpha, beta, nu, rho);
85 } else {
86 return unsafeSabrLogNormalVolatility(strike + shift, forward + shift, expiryTime, alpha, beta, nu, rho);
87 }
88 }
89
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(forward * strike, oneMinusBeta);
95 const Real sqrtA = std::sqrt(A);
96 Real logM;
97 if (!close(forward, strike))
98 logM = std::log(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(B) + z - rho) / (1.0 - rho);
108 const Real xx = std::log(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(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(forward * strike, beta / 2.0);
126
127 return F * E * multiplier * d;
128 }
129
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
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
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");
174 return unsafeSabrVolatility(strike, forward, expiryTime, alpha, beta, nu, rho,
175 volatilityType);
176 }
177
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");
194 return unsafeShiftedSabrVolatility(strike, forward, expiryTime,
195 alpha, beta, nu, rho,shift, volatilityType);
196 }
197
198 namespace {
199 struct SabrFlochKennedyVolatility {
201
202 Real y(Real k) const {
203 return -1.0/(1.0-beta)*(std::pow(F,1-beta)-std::pow(k,1-beta));
204 }
205
206 Real Dint(Real k) const {
207 return 1/nu*std::log( ( std::sqrt(1+2*rho*nu/alpha*y(k)
208 + squared(nu/alpha*y(k)) )
209 - rho - nu/alpha*y(k) ) / (1-rho) );
210 }
211
212 Real D(Real k) const {
213 return std::sqrt(alpha*alpha+2*alpha*rho*nu*y(k)
214 + squared(nu*y(k)))*std::pow(k,beta);
215 }
216
217 Real omega0(Real k) const {
218 return std::log(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(k,beta)-std::pow(F,beta))/(k-F)*t)
226 -omega0(k)/squared(Dint(k))*(std::log(
227 omega0(k)) + 0.5*std::log((F*k/(D(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(F,-3 + beta)*(alpha2*squared(-1 + beta)*std::pow(F,2*beta)*t + 6*alpha*beta*nu*std::pow(F,1 + beta)*rho*t +
240 F2*(24 + nu*nu*(2 - 3*rho2)*t)))/24.0 +
241 (3*alpha2*alpha*std::pow(-1 + beta,3)*std::pow(F,3*beta)*t +
242 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) +
243 alpha*std::pow(F,2 + beta)*(24*(-1 + beta) + nu*nu*(2*(-1 + beta) + 3*(1 + beta)*rho2)*t))/(48.*F2*F2) * (k-F) +
244 (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 -
245 30*alpha*nu*std::pow(F,3 + beta)*rho*(24 + nu*nu*(-4*(1 + beta) + 3*(1 + 2*beta)*rho2)*t) +
246 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) +
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
253 Rate forward,
254 Time expiryTime,
255 Real alpha,
256 Real beta,
257 Real nu,
258 Real rho) {
259 const SabrFlochKennedyVolatility v =
260 {forward, alpha, beta, nu, rho, expiryTime};
261
262 return v(strike);
263 }
264
265}
floating-point comparisons
output manipulators
const DefaultType & t
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Date d
#define QL_EPSILON
Definition: qldefines.hpp:178
detail::percent_holder rate(Rate)
output rates and spreads as percentages
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real Rate
interest rates
Definition: types.hpp:70
Real rho
functionals and combinators not included in the STL
Definition: any.hpp:35
Real sabrVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho, VolatilityType volatilityType)
Definition: sabr.cpp:159
Real unsafeShiftedSabrVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho, Real shift, VolatilityType volatilityType)
Definition: sabr.cpp:74
Real unsafeSabrVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho, VolatilityType volatilityType)
Definition: sabr.cpp:130
Real sabrFlochKennedyVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho)
Definition: sabr.cpp:252
Real unsafeSabrNormalVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho)
Definition: sabr.cpp:90
Real unsafeSabrLogNormalVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho)
Definition: sabr.cpp:33
T squared(T x)
Definition: functional.hpp:37
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
void validateSabrParameters(Real alpha, Real beta, Real nu, Real rho)
Definition: sabr.cpp:145
Real shiftedSabrVolatility(Rate strike, Rate forward, Time expiryTime, Real alpha, Real beta, Real nu, Real rho, Real shift, VolatilityType volatilityType)
Definition: sabr.cpp:178
ext::shared_ptr< BlackVolTermStructure > v
Real beta
Definition: sabr.cpp:200
Real F
Definition: sabr.cpp:200
Real nu
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
SABR functions.
volatility types