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
20#include <ql/experimental/volatility/noarbsabr.hpp>
21
22#include <ql/math/solvers1d/brent.hpp>
23#include <ql/math/modifiedbessel.hpp>
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(a: f - strike, b: 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
57 QL_REQUIRE(expiryTime > 0.0 && expiryTime <= detail::NoArbSabrModel::expiryTime_max,
58 "expiryTime (" << expiryTime << ") out of bounds");
59 QL_REQUIRE(forward > 0.0, "forward (" << forward << ") must be positive");
60 QL_REQUIRE(beta >= detail::NoArbSabrModel::beta_min && beta <= detail::NoArbSabrModel::beta_max,
61 "beta (" << beta << ") out of bounds");
62 Real sigmaI = alpha * std::pow(x: forward, y: beta - 1.0);
63 QL_REQUIRE(sigmaI >= detail::NoArbSabrModel::sigmaI_min &&
64 sigmaI <= detail::NoArbSabrModel::sigmaI_max,
65 "sigmaI = alpha*forward^(beta-1.0) ("
66 << sigmaI << ") out of bounds, alpha=" << alpha
67 << " beta=" << beta << " forward=" << forward);
68 QL_REQUIRE(nu >= detail::NoArbSabrModel::nu_min && nu <= detail::NoArbSabrModel::nu_max,
69 "nu (" << nu << ") out of bounds");
70 QL_REQUIRE(rho >= detail::NoArbSabrModel::rho_min && rho <= detail::NoArbSabrModel::rho_max,
71 "rho (" << rho << ") out of bounds");
72
73 // determine a region sufficient for integration in the normal case
74
75 fmin_ = fmax_ = forward_;
76 for (Real tmp = p(f: fmax_);
77 tmp > std::max(a: detail::NoArbSabrModel::i_accuracy / std::max(a: 1.0, b: fmax_ - fmin_),
78 b: detail::NoArbSabrModel::density_threshold);
79 tmp = p(f: fmax_)) {
80 fmax_ *= 2.0;
81 }
82 for (Real tmp = p(f: fmin_);
83 tmp > std::max(a: detail::NoArbSabrModel::i_accuracy / std::max(a: 1.0, b: fmax_ - fmin_),
84 b: detail::NoArbSabrModel::density_threshold);
85 tmp = p(f: fmin_)) {
86 fmin_ *= 0.5;
87 }
88 fmin_ = std::max(a: detail::NoArbSabrModel::strike_min, b: fmin_);
89
90 QL_REQUIRE(fmax_ > fmin_, "could not find a reasonable integration domain");
91
92 integrator_ =
93 ext::make_shared<GaussLobattoIntegral>(
94 args: detail::NoArbSabrModel::i_max_iterations, args: detail::NoArbSabrModel::i_accuracy);
95
96 detail::D0Interpolator d0(forward_, expiryTime_, alpha_, beta_, nu_, rho_);
97 absProb_ = d0();
98
99 try {
100 Brent b;
101 Real start = std::sqrt(x: externalForward_ - detail::NoArbSabrModel::strike_min);
102 Real tmp =
103 b.solve(f: [&](Real x){ return forwardError(forward: x); },
104 accuracy: detail::NoArbSabrModel::forward_accuracy, guess: start,
105 step: std::min(a: detail::NoArbSabrModel::forward_search_step, b: start / 2.0));
106 forward_ = tmp * tmp + detail::NoArbSabrModel::strike_min;
107 } catch (Error&) {
108 // fall back to unadjusted forward
109 forward_ = externalForward_;
110 }
111
112 Real d = forwardError(forward: std::sqrt(x: forward_ - detail::NoArbSabrModel::strike_min));
113 numericalForward_ = d + externalForward_;
114}
115
116Real NoArbSabrModel::optionPrice(const Real strike) const {
117 if (p(f: std::max(a: forward_, b: strike)) < detail::NoArbSabrModel::density_threshold)
118 return 0.0;
119 return (1.0 - absProb_) *
120 ((*integrator_)(integrand(this, strike),
121 strike, std::max(a: fmax_, b: 2.0 * strike)) /
122 numericalIntegralOverP_);
123}
124
125Real NoArbSabrModel::digitalOptionPrice(const Real strike) const {
126 if (strike < QL_MIN_POSITIVE_REAL)
127 return 1.0;
128 if (p(f: std::max(a: forward_, b: strike)) < detail::NoArbSabrModel::density_threshold)
129 return 0.0;
130 return (1.0 - absProb_)
131 * ((*integrator_)(p_integrand(this),
132 strike, std::max(a: fmax_, b: 2.0 * strike)) /
133 numericalIntegralOverP_);
134}
135
136Real NoArbSabrModel::forwardError(const Real forward) const {
137 forward_ = forward * forward + detail::NoArbSabrModel::strike_min;
138 numericalIntegralOverP_ = (*integrator_)(p_integrand(this),
139 fmin_, fmax_);
140 return optionPrice(strike: 0.0) - externalForward_;
141}
142
143Real NoArbSabrModel::p(const Real f) const {
144
145 if (f < detail::NoArbSabrModel::density_lower_bound ||
146 forward_ < detail::NoArbSabrModel::density_lower_bound)
147 return 0.0;
148
149 Real fOmB = std::pow(x: f, y: 1.0 - beta_);
150 Real FOmB = std::pow(x: forward_, y: 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(x: 1.0 + 2.0 * rho_ * nu_ * zf + nu_ * nu_ * zf * zf);
158 Real Jz = std::sqrt(x: 1.0 - 2.0 * rho_ * nu_ * z + nu_ * nu_ * z * z);
159
160 Real xz = std::log(x: (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(x: 1.0 - rho_ * rho_);
168 Real h = 0.5 * beta_ * rho_ / ((1.0 - beta_) * Jmzf * Jmzf) *
169 (nu_ * zf * std::log(x: zf * Jz / zF) +
170 (1 + rho_ * nu_ * zf) / sqrtOmR *
171 (std::atan(x: (nu_ * z - rho_) / sqrtOmR) +
172 std::atan(x: rho_ / sqrtOmR)));
173
174 Real res =
175 std::pow(x: Jz, y: -1.5) / (alpha_ * std::pow(x: f, y: beta_) * expiryTime_) *
176 std::pow(x: zf, y: 1.0 - gamma) * std::pow(x: zF, y: gamma) *
177 std::exp(x: -(xz * xz) / (2.0 * expiryTime_) +
178 (h + kappa1 * expiryTime_)) *
179 modifiedBesselFunction_i_exponentiallyWeighted(nu: gamma,
180 x: 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(x: forward_, y: 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
220Real D0Interpolator::operator()() const {
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(first: tauG_.begin(), last: tauG_.end(), val: 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(first: sigmaIG_.rbegin(), last: sigmaIG_.rend(), val: 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(first: rhoG_.rbegin(), last: rhoG_.rend(), val: 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(first: nuG_.begin(), last: nuG_.end(), val: 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(first: betaG_.begin(), last: betaG_.end(), val: 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 =
291 phi(d0: detail::NoArbSabrModel::tiny_prob);
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(d0: (Real)sabrabsprob[ind] /
305 detail::NoArbSabrModel::nsim);
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(phi: phiRes);
319}
320
321Real D0Interpolator::phi(const Real d0) const {
322 if (d0 < 1e-14)
323 return detail::NoArbSabrModel::phiByTau_cutoff * expiryTime_;
324 return boost::math::gamma_q_inv(a: gamma_, p: d0) * expiryTime_;
325}
326
327Real D0Interpolator::d0(const Real phi) const {
328 return boost::math::gamma_q(a: gamma_, z: std::max(a: 0.0, b: phi / expiryTime_));
329}
330
331} // namespace detail
332
333} // namespace QuantLib
334

source code of quantlib/ql/experimental/volatility/noarbsabr.cpp