QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
noarbsabr.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) 2014 Peter Caspers
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
21
24#include <boost/math/special_functions/gamma.hpp>
25#include <boost/functional/hash.hpp>
26
27namespace QuantLib {
28
29class NoArbSabrModel::integrand {
30 const NoArbSabrModel* model;
31 Real strike;
32 public:
33 integrand(const NoArbSabrModel* model, Real strike)
34 : model(model), strike(strike) {}
35 Real operator()(Real f) const {
36 return std::max(f - strike, 0.0) * model->p(f);
37 }
38};
39
40class NoArbSabrModel::p_integrand {
41 const NoArbSabrModel* model;
42 public:
43 explicit p_integrand(const NoArbSabrModel* model)
44 : model(model) {}
45 Real operator()(Real f) const {
46 return model->p(f);
47 }
48};
49
50NoArbSabrModel::NoArbSabrModel(const Real expiryTime, const Real forward,
51 const Real alpha, const Real beta, const Real nu,
52 const Real rho)
53 : expiryTime_(expiryTime), externalForward_(forward), alpha_(alpha),
54 beta_(beta), nu_(nu), rho_(rho), forward_(forward),
55 numericalForward_(forward) {
56
58 "expiryTime (" << expiryTime << ") out of bounds");
59 QL_REQUIRE(forward > 0.0, "forward (" << forward << ") must be positive");
61 "beta (" << beta << ") out of bounds");
62 Real sigmaI = alpha * std::pow(forward, beta - 1.0);
65 "sigmaI = alpha*forward^(beta-1.0) ("
66 << sigmaI << ") out of bounds, alpha=" << alpha
67 << " beta=" << beta << " forward=" << forward);
69 "nu (" << nu << ") out of bounds");
71 "rho (" << rho << ") out of bounds");
72
73 // determine a region sufficient for integration in the normal case
74
76 for (Real tmp = p(fmax_);
77 tmp > std::max(detail::NoArbSabrModel::i_accuracy / std::max(1.0, fmax_ - fmin_),
79 tmp = p(fmax_)) {
80 fmax_ *= 2.0;
81 }
82 for (Real tmp = p(fmin_);
83 tmp > std::max(detail::NoArbSabrModel::i_accuracy / std::max(1.0, fmax_ - fmin_),
85 tmp = p(fmin_)) {
86 fmin_ *= 0.5;
87 }
89
90 QL_REQUIRE(fmax_ > fmin_, "could not find a reasonable integration domain");
91
93 ext::make_shared<GaussLobattoIntegral>(
95
97 absProb_ = d0();
98
99 try {
100 Brent b;
102 Real tmp =
103 b.solve([&](Real x){ return forwardError(x); },
105 std::min(detail::NoArbSabrModel::forward_search_step, start / 2.0));
107 } catch (Error&) {
108 // fall back to unadjusted forward
110 }
111
114}
115
118 return 0.0;
119 return (1.0 - absProb_) *
120 ((*integrator_)(integrand(this, strike),
121 strike, std::max(fmax_, 2.0 * strike)) /
123}
124
126 if (strike < QL_MIN_POSITIVE_REAL)
127 return 1.0;
129 return 0.0;
130 return (1.0 - absProb_)
131 * ((*integrator_)(p_integrand(this),
132 strike, std::max(fmax_, 2.0 * strike)) /
134}
135
138 numericalIntegralOverP_ = (*integrator_)(p_integrand(this),
139 fmin_, fmax_);
140 return optionPrice(0.0) - externalForward_;
141}
142
144
147 return 0.0;
148
149 Real fOmB = std::pow(f, 1.0 - beta_);
150 Real FOmB = std::pow(forward_, 1.0 - beta_);
151
152 Real zf = fOmB / (alpha_ * (1.0 - beta_));
153 Real zF = FOmB / (alpha_ * (1.0 - beta_));
154 Real z = zF - zf;
155
156 // Real JzF = std::sqrt(1.0 - 2.0 * rho_ * nu_ * zF + nu_ * nu_ * zF * zF);
157 Real Jmzf = std::sqrt(1.0 + 2.0 * rho_ * nu_ * zf + nu_ * nu_ * zf * zf);
158 Real Jz = std::sqrt(1.0 - 2.0 * rho_ * nu_ * z + nu_ * nu_ * z * z);
159
160 Real xz = std::log((Jz - rho_ + nu_ * z) / (1.0 - rho_)) / nu_;
161 Real Bp_B = beta_ / FOmB;
162 // Real Bpp_B = beta_ * (2.0 * beta_ - 1.0) / (FOmB * FOmB);
163 Real kappa1 = 0.125 * nu_ * nu_ * (2.0 - 3.0 * rho_ * rho_) -
164 0.25 * rho_ * nu_ * alpha_ * Bp_B;
165 // Real kappa2 = alpha_ * alpha_ * (0.25 * Bpp_B - 0.375 * Bp_B * Bp_B);
166 Real gamma = 1.0 / (2.0 * (1.0 - beta_));
167 Real sqrtOmR = std::sqrt(1.0 - rho_ * rho_);
168 Real h = 0.5 * beta_ * rho_ / ((1.0 - beta_) * Jmzf * Jmzf) *
169 (nu_ * zf * std::log(zf * Jz / zF) +
170 (1 + rho_ * nu_ * zf) / sqrtOmR *
171 (std::atan((nu_ * z - rho_) / sqrtOmR) +
172 std::atan(rho_ / sqrtOmR)));
173
174 Real res =
175 std::pow(Jz, -1.5) / (alpha_ * std::pow(f, beta_) * expiryTime_) *
176 std::pow(zf, 1.0 - gamma) * std::pow(zF, gamma) *
177 std::exp(-(xz * xz) / (2.0 * expiryTime_) +
178 (h + kappa1 * expiryTime_)) *
180 Real(zF * zf / expiryTime_));
181 return res;
182}
183
184namespace detail {
185
186D0Interpolator::D0Interpolator(const Real forward, const Real expiryTime,
187 const Real alpha, const Real beta, const Real nu,
188 const Real rho)
189 : forward_(forward), expiryTime_(expiryTime), alpha_(alpha), beta_(beta),
190 nu_(nu), rho_(rho), gamma_(1.0 / (2.0 * (1.0 - beta_))) {
191
192 sigmaI_ = alpha_ * std::pow(forward_, beta_ - 1.0);
193
194 tauG_ = {
195 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0,
196 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0, 5.25, 5.5, 5.75, 6.0, 6.25,
197 6.5, 6.75, 7.0, 7.25, 7.5, 7.75, 8.0, 8.25, 8.5, 8.75, 9.0, 9.25, 9.5,
198 9.75, 10.0, 10.25, 10.5, 10.75, 11.0, 11.25, 11.5, 11.75, 12.0, 12.25,
199 12.5, 12.75, 13.0, 13.25, 13.5, 13.75, 14.0, 14.25, 14.5, 14.75, 15.0,
200 15.25, 15.5, 15.75, 16.0, 16.25, 16.5, 16.75, 17.0, 17.25, 17.5, 17.75,
201 18.0, 18.25, 18.5, 18.75, 19.0, 19.25, 19.5, 19.75, 20.0, 20.25, 20.5,
202 20.75, 21.0, 21.25, 21.5, 21.75, 22.0, 22.25, 22.5, 22.75, 23.0, 23.25,
203 23.5, 23.75, 24.0, 24.25, 24.5, 24.75, 25.0, 25.25, 25.5, 25.75, 26.0,
204 26.25, 26.5, 26.75, 27.0, 27.25, 27.5, 27.75, 28.0, 28.25, 28.5, 28.75,
205 29.0, 29.25, 29.5, 29.75, 30.0
206 };
207
208 sigmaIG_ = {
209 1.0, 0.8, 0.7, 0.6, 0.5, 0.45, 0.4, 0.35, 0.3, 0.27, 0.24, 0.21,
210 0.18, 0.15, 0.125, 0.1, 0.075, 0.05
211 };
212
213 rhoG_ = { 0.75, 0.50, 0.25, 0.00, -0.25, -0.50, -0.75 };
214
215 nuG_ = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 };
216
217 betaG_ = { 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
218}
219
221
222 // we do not need to check the indices here, because this is already
223 // done in the NoArbSabr constructor
224
225 Size tauInd = std::upper_bound(tauG_.begin(), tauG_.end(), expiryTime_) -
226 tauG_.begin();
227 if (tauInd == tauG_.size())
228 --tauInd; // tau at upper bound
229 Real expiryTimeTmp = expiryTime_;
230 if (tauInd == 0) {
231 ++tauInd;
232 expiryTimeTmp = tauG_.front();
233 }
234 Real tauL = (expiryTimeTmp - tauG_[tauInd - 1]) /
235 (tauG_[tauInd] - tauG_[tauInd - 1]);
236
237 Size sigmaIInd =
238 sigmaIG_.size() -
239 (std::upper_bound(sigmaIG_.rbegin(), sigmaIG_.rend(), sigmaI_) -
240 sigmaIG_.rbegin());
241 if (sigmaIInd == 0)
242 ++sigmaIInd; // sigmaI at upper bound
243 Real sigmaIL = (sigmaI_ - sigmaIG_[sigmaIInd - 1]) /
244 (sigmaIG_[sigmaIInd] - sigmaIG_[sigmaIInd - 1]);
245
246 Size rhoInd =
247 rhoG_.size() -
248 (std::upper_bound(rhoG_.rbegin(), rhoG_.rend(), rho_) - rhoG_.rbegin());
249 if (rhoInd == 0) {
250 rhoInd++;
251 }
252 if (rhoInd == rhoG_.size()) {
253 rhoInd--;
254 }
255 Real rhoL =
256 (rho_ - rhoG_[rhoInd - 1]) / (rhoG_[rhoInd] - rhoG_[rhoInd - 1]);
257
258 // for nu = 0 we know phi = 0.5*z_F^2
259 Size nuInd = std::upper_bound(nuG_.begin(), nuG_.end(), nu_) - nuG_.begin();
260 if (nuInd == nuG_.size())
261 --nuInd; // nu at upper bound
262 Real tmpNuG = nuInd > 0 ? nuG_[nuInd - 1] : 0.0;
263 Real nuL = (nu_ - tmpNuG) / (nuG_[nuInd] - tmpNuG);
264
265 // for beta = 1 we know phi = 0.0
266 Size betaInd =
267 std::upper_bound(betaG_.begin(), betaG_.end(), beta_) - betaG_.begin();
268 Real tmpBetaG;
269 if (betaInd == betaG_.size())
270 tmpBetaG = 1.0;
271 else
272 tmpBetaG = betaG_[betaInd];
273 Real betaL =
274 (beta_ - betaG_[betaInd - 1]) / (tmpBetaG - betaG_[betaInd - 1]);
275
276 Real phiRes = 0.0;
277 for (int iTau = -1; iTau <= 0; ++iTau) {
278 for (int iSigma = -1; iSigma <= 0; ++iSigma) {
279 for (int iRho = -1; iRho <= 0; ++iRho) {
280 for (int iNu = -1; iNu <= 0; ++iNu) {
281 for (int iBeta = -1; iBeta <= 0; ++iBeta) {
282 Real phiTmp;
283 if (iNu == -1 && nuInd == 0) {
284 phiTmp =
285 0.5 /
286 (sigmaI_ * sigmaI_ * (1.0 - beta_) *
287 (1.0 - beta_)); // this is 0.5*z_F^2, see above
288 } else {
289 if (iBeta == 0 && betaInd == betaG_.size()) {
290 phiTmp =
292 } else {
293 int ind = (tauInd + iTau +
294 (sigmaIInd + iSigma +
295 (rhoInd + iRho +
296 (nuInd + iNu + ((betaInd + iBeta) *
297 nuG_.size())) *
298 rhoG_.size()) *
299 sigmaIG_.size()) *
300 tauG_.size());
301 QL_REQUIRE(ind >= 0 && ind < 1209600,
302 "absorption matrix index ("
303 << ind << ") invalid");
304 phiTmp = phi((Real)sabrabsprob[ind] /
306 }
307 }
308 phiRes += phiTmp * (iTau == -1 ? (1.0 - tauL) : tauL) *
309 (iSigma == -1 ? (1.0 - sigmaIL) : sigmaIL) *
310 (iRho == -1 ? (1.0 - rhoL) : rhoL) *
311 (iNu == -1 ? (1.0 - nuL) : nuL) *
312 (iBeta == -1 ? (1.0 - betaL) : betaL);
313 }
314 }
315 }
316 }
317 }
318 return d0(phiRes);
319}
320
322 if (d0 < 1e-14)
324 return boost::math::gamma_q_inv(gamma_, d0) * expiryTime_;
325}
326
327Real D0Interpolator::d0(const Real phi) const {
328 return boost::math::gamma_q(gamma_, std::max(0.0, phi / expiryTime_));
329}
330
331} // namespace detail
332
333} // namespace QuantLib
Brent 1-D solver.
Brent 1-D solver
Definition: brent.hpp:37
Base error class.
Definition: errors.hpp:39
Real p(Real f) const
Definition: noarbsabr.cpp:143
NoArbSabrModel(Real expiryTime, Real forward, Real alpha, Real beta, Real nu, Real rho)
Definition: noarbsabr.cpp:50
Real digitalOptionPrice(Real strike) const
Definition: noarbsabr.cpp:125
Real optionPrice(Real strike) const
Definition: noarbsabr.cpp:116
ext::shared_ptr< GaussLobattoIntegral > integrator_
Definition: noarbsabr.hpp:131
Real forwardError(Real forward) const
Definition: noarbsabr.cpp:136
std::vector< Real > sigmaIG_
Definition: noarbsabr.hpp:150
D0Interpolator(Real forward, Real expiryTime, Real alpha, Real beta, Real nu, Real rho)
Definition: noarbsabr.cpp:186
std::vector< Real > betaG_
Definition: noarbsabr.hpp:150
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Date d
ext::function< Real(Real)> b
#define QL_MIN_POSITIVE_REAL
Definition: qldefines.hpp:177
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Real rho
modified Bessel functions of first and second kind
const unsigned long sabrabsprob[1209600]
Definition: noarbsabr.hpp:138
Definition: any.hpp:35
Real modifiedBesselFunction_i_exponentiallyWeighted(Real nu, Real x)
No-arbitrage SABR.
Real beta
Definition: sabr.cpp:200
Real nu
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200