QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
qdfpamericanengine.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) 2022 Klaus Spanderen
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/*! \file qrfpamericanengine.cpp
21*/
22
30#include <utility>
31#ifndef QL_BOOST_HAS_TANH_SINH
33#endif
34
35namespace QuantLib {
36
38 Size l, Size m, Size n, Size p):
39 m_(m), n_(n),
40 fpIntegrator_(ext::make_shared<GaussLegendreIntegrator>(l)),
41 exerciseBoundaryIntegrator_(
42 ext::make_shared<GaussLegendreIntegrator>(p)) {
43
44 QL_REQUIRE(m_ > 0, "at least one fixed point iteration step is needed");
45 QL_REQUIRE(n_ > 0, "at least one interpolation point is needed");
46 }
47
49 const {
50 return n_;
51 }
53 return m_-1;
54 }
56 const {
57 return Size(1);
58 }
59
60 ext::shared_ptr<Integrator>
62 return fpIntegrator_;
63 }
64 ext::shared_ptr<Integrator>
66 const {
68 }
69
71 Size m, Size n, Real eps)
72 : m_(m), n_(n),
73#ifdef QL_BOOST_HAS_TANH_SINH
74 integrator_(ext::make_shared<TanhSinhIntegral>(eps))
75#else
76 integrator_(ext::make_shared<GaussLobattoIntegral>(
77 100000, QL_MAX_REAL, 0.1*eps))
78#endif
79 {}
80
82 const {
83 return n_;
84 }
86 return m_-1;
87 }
89 const {
90 return Size(1);
91 }
92
93 ext::shared_ptr<Integrator>
95 return integrator_;
96 }
97 ext::shared_ptr<Integrator>
99 const {
100 return integrator_;
101 }
102
104 Size l, Size m, Size n, Real eps)
105 : QdFpLegendreScheme(l, m, n, 1),
106 eps_(eps) {}
107
108 ext::shared_ptr<Integrator>
110#ifdef QL_BOOST_HAS_TANH_SINH
111 return ext::make_shared<TanhSinhIntegral>(eps_);
112#else
113 return ext::make_shared<GaussLobattoIntegral>(
114 100000, QL_MAX_REAL, 0.1*eps_);
115#endif
116 }
117
118
119 class DqFpEquation {
120 public:
121 DqFpEquation(Rate _r,
122 Rate _q,
123 Volatility _vol,
124 std::function<Real(Real)> B,
125 ext::shared_ptr<Integrator> _integrator)
126 : r(_r), q(_q), vol(_vol), B(std::move(B)), integrator(std::move(_integrator)) {
127 const auto legendreIntegrator =
128 ext::dynamic_pointer_cast<GaussLegendreIntegrator>(integrator);
129
130 if (legendreIntegrator != nullptr) {
131 x_i = legendreIntegrator->getIntegration()->x();
132 w_i = legendreIntegrator->getIntegration()->weights();
133 }
134 }
135
136 virtual std::pair<Real, Real> NDd(Real tau, Real b) const = 0;
137 virtual std::tuple<Real, Real, Real> f(Real tau, Real b) const = 0;
138
139 virtual ~DqFpEquation() = default;
140
141 protected:
142 std::pair<Real, Real> d(Time t, Real z) const {
143 const Real v = vol * std::sqrt(t);
144 const Real m = (std::log(z) + (r-q)*t)/v + 0.5*v;
145
146 return std::make_pair(m, m-v);
147 }
148
149 Array x_i, w_i;
150 const Rate r, q;
151 const Volatility vol;
152
153 const std::function<Real(Real)> B;
154 const ext::shared_ptr<Integrator> integrator;
155
156 const NormalDistribution phi;
157 const CumulativeNormalDistribution Phi;
158 };
159
160 class DqFpEquation_B: public DqFpEquation {
161 public:
162 DqFpEquation_B(Real K,
163 Rate _r,
164 Rate _q,
165 Volatility _vol,
166 std::function<Real(Real)> B,
167 ext::shared_ptr<Integrator> _integrator);
168
169 std::pair<Real, Real> NDd(Real tau, Real b) const override;
170 std::tuple<Real, Real, Real> f(Real tau, Real b) const override;
171
172 private:
173 const Real K;
174 };
175
176 class DqFpEquation_A: public DqFpEquation {
177 public:
178 DqFpEquation_A(Real K,
179 Rate _r,
180 Rate _q,
181 Volatility _vol,
182 std::function<Real(Real)> B,
183 ext::shared_ptr<Integrator> _integrator);
184
185 std::pair<Real, Real> NDd(Real tau, Real b) const override;
186 std::tuple<Real, Real, Real> f(Real tau, Real b) const override;
187
188 private:
189 const Real K;
190 };
191
192 DqFpEquation_A::DqFpEquation_A(Real K,
193 Rate _r,
194 Rate _q,
195 Volatility _vol,
196 std::function<Real(Real)> B,
197 ext::shared_ptr<Integrator> _integrator)
198 : DqFpEquation(_r, _q, _vol, std::move(B), std::move(_integrator)), K(K) {}
199
200 std::tuple<Real, Real, Real> DqFpEquation_A::f(Real tau, Real b) const {
201 const Real v = vol * std::sqrt(tau);
202
203 Real N, D;
204 if (tau < squared(QL_EPSILON)) {
205 if (close_enough(b, K)) {
206 N = 1/(M_SQRT2*M_SQRTPI * v);
207 D = N + 0.5;
208 }
209 else {
210 N = 0.0;
211 D = (b > K)? 1.0: 0.0;
212 }
213 }
214 else {
215 const Real stv = std::sqrt(tau)/vol;
216
217 Real K12, K3;
218 if (!x_i.empty()) {
219 K12 = K3 = 0.0;
220
221 for (Integer i = x_i.size()-1; i >= 0; --i) {
222 const Real y = x_i[i];
223 const Real m = 0.25*tau*squared(1+y);
224 const std::pair<Real, Real> dpm = d(m, b/B(tau-m));
225
226 K12 += w_i[i] * std::exp(q*tau - q*m)
227 *(0.5*tau*(y+1)*Phi(dpm.first) + stv*phi(dpm.first));
228 K3 += w_i[i] * stv*std::exp(r*tau-r*m)*phi(dpm.second);
229 }
230 } else {
231 K12 = (*integrator)([&, this](Real y) -> Real {
232 const Real m = 0.25*tau*squared(1+y);
233 const Real df = std::exp(q*tau - q*m);
234
235 if (y <= 5*QL_EPSILON - 1) {
236 if (close_enough(b, B(tau-m)))
237 return df*stv/(M_SQRT2*M_SQRTPI);
238 else
239 return 0.0;
240 }
241 else {
242 const Real dp = d(m, b/B(tau-m)).first;
243 return df*(0.5*tau*(y+1)*Phi(dp) + stv*phi(dp));
244 }
245 }, -1, 1);
246
247 K3 = (*integrator)([&, this](Real y) -> Real {
248 const Real m = 0.25*tau*squared(1+y);
249 const Real df = std::exp(r*tau-r*m);
250
251 if (y <= 5*QL_EPSILON - 1) {
252 if (close_enough(b, B(tau-m)))
253 return df*stv/(M_SQRT2*M_SQRTPI);
254 else
255 return 0.0;
256 }
257 else
258 return df*stv*phi(d(m, b/B(tau-m)).second);
259 }, -1, 1);
260 }
261 const std::pair<Real, Real> dpm = d(tau, b/K);
262 N = phi(dpm.second)/v + r*K3;
263 D = phi(dpm.first)/v + Phi(dpm.first) + q*K12;
264 }
265
266 const Real alpha = K*std::exp(-(r-q)*tau);
267 Real fv;
268 if (tau < squared(QL_EPSILON)) {
269 if (close_enough(b, K))
270 fv = alpha;
271 else if (b > K)
272 fv = 0.0;
273 else {
274 if (close_enough(q, Real(0)))
275 fv = alpha*r*((q < 0)? -1.0 : 1.0)/QL_EPSILON;
276 else
277 fv = alpha*r/q;
278 }
279 }
280 else
281 fv = alpha*N/D;
282
283 return std::make_tuple(N, D, fv);
284 }
285
286 std::pair<Real, Real> DqFpEquation_A::NDd(Real tau, Real b) const {
287 Real Dd, Nd;
288
289 if (tau < squared(QL_EPSILON)) {
290 if (close_enough(b, K)) {
291 const Real sqTau = std::sqrt(tau);
292 const Real vol2 = vol*vol;
293 Dd = M_1_SQRTPI*M_SQRT_2*(
294 -(0.5*vol2 + r-q) / (b*vol*vol2*sqTau) + 1 / (b*vol*sqTau));
295 Nd = M_1_SQRTPI*M_SQRT_2 * (-0.5*vol2 + r-q) / (b*vol*vol2*sqTau);
296 }
297 else
298 Dd = Nd = 0.0;
299 }
300 else {
301 const std::pair<Real, Real> dpm = d(tau, b/K);
302
303 Dd = -phi(dpm.first) * dpm.first / (b*vol*vol*tau) +
304 phi(dpm.first) / (b*vol * std::sqrt(tau));
305 Nd = -phi(dpm.second) * dpm.second / (b*vol*vol*tau);
306 }
307
308 return std::make_pair(Nd, Dd);
309 }
310
311
312 DqFpEquation_B::DqFpEquation_B(Real K,
313 Rate _r,
314 Rate _q,
315 Volatility _vol,
316 std::function<Real(Real)> B,
317 ext::shared_ptr<Integrator> _integrator)
318 : DqFpEquation(_r, _q, _vol, std::move(B), std::move(_integrator)), K(K) {}
319
320
321 std::tuple<Real, Real, Real> DqFpEquation_B::f(Real tau, Real b) const {
322 Real N, D;
323 if (tau < squared(QL_EPSILON)) {
324 if (close_enough(b, K))
325 N = D = 0.5;
326 else if (b < K)
327 N = D = 0.0;
328 else
329 N = D = 1.0;
330 }
331 else {
332 Real ni, di;
333 if (!x_i.empty()) {
334 const Real c = 0.5*tau;
335
336 ni = di = 0.0;
337 for (Integer i = x_i.size()-1; i >= 0; --i) {
338 const Real u = c*x_i[i] + c;
339 const std::pair<Real, Real> dpm = d(tau - u, b/B(u));
340 ni += w_i[i] * std::exp(r*u)*Phi(dpm.second);
341 di += w_i[i] * std::exp(q*u)*Phi(dpm.first);
342 }
343 ni *= c;
344 di *= c;
345 } else {
346 ni = (*integrator)([&, this](Real u) -> Real {
347 const Real df = std::exp(r*u);
348 if (u >= tau*(1 - 5*QL_EPSILON)) {
349 if (close_enough(b, B(u)))
350 return 0.5*df;
351 else
352 return df*((b < B(u)? 0.0: 1.0));
353 }
354 else
355 return df*Phi(d(tau - u, b/B(u)).second);
356 }, 0, tau);
357 di = (*integrator)([&, this](Real u) -> Real {
358 const Real df = std::exp(q*u);
359 if (u >= tau*(1 - 5*QL_EPSILON)) {
360 if (close_enough(b, B(u)))
361 return 0.5*df;
362 else
363 return df*((b < B(u)? 0.0: 1.0));
364 }
365 else
366 return df*Phi(d(tau - u, b/B(u)).first);
367 }, 0, tau);
368 }
369
370 const std::pair<Real, Real> dpm = d(tau, b/K);
371
372 N = Phi(dpm.second) + r*ni;
373 D = Phi(dpm.first) + q*di;
374 }
375
376 Real fv;
377 const Real alpha = K*std::exp(-(r-q)*tau);
378 if (tau < squared(QL_EPSILON)) {
379 if (close_enough(b, K) || b > K)
380 fv = alpha;
381 else {
382 if (close_enough(q, Real(0)))
383 fv = alpha*r*((q < 0)? -1.0 : 1.0)/QL_EPSILON;
384 else
385 fv = alpha*r/q;
386 }
387 }
388 else
389 fv = alpha*N/D;
390
391 return std::make_tuple(N, D, fv);
392 }
393
394 std::pair<Real, Real> DqFpEquation_B::NDd(Real tau, Real b) const {
395 const std::pair<Real, Real> dpm = d(tau, b/K);
396 return std::make_pair(
397 phi(dpm.second) / (b*vol*std::sqrt(tau)),
398 phi(dpm.first) / (b*vol*std::sqrt(tau))
399 );
400 }
401
403 ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess,
404 ext::shared_ptr<QdFpIterationScheme> iterationScheme,
405 FixedPointEquation fpEquation)
406 : detail::QdPutCallParityEngine(std::move(bsProcess)),
407 iterationScheme_(std::move(iterationScheme)),
408 fpEquation_(fpEquation) {
409 }
410
411 ext::shared_ptr<QdFpIterationScheme>
413 static auto scheme = ext::make_shared<QdFpLegendreScheme>(7, 2, 7, 27);
414 return scheme;
415 }
416
417 ext::shared_ptr<QdFpIterationScheme>
419 static auto scheme = ext::make_shared<QdFpLegendreTanhSinhScheme>(25, 5, 13, 1e-8);
420 return scheme;
421 }
422
423 ext::shared_ptr<QdFpIterationScheme>
425 static auto scheme = ext::make_shared<QdFpTanhSinhIterationScheme>(10, 30, 1e-10);
426 return scheme;
427 }
428
430 Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
431
432 if (r < 0.0 && q < r)
433 QL_FAIL("double-boundary case q<r<0 for a put option is given");
434
435 const Real xmax = QdPlusAmericanEngine::xMax(K, r, q);
436 const Size n = iterationScheme_->getNumberOfChebyshevInterpolationNodes();
437
438 const ext::shared_ptr<ChebyshevInterpolation> interp =
441 .getPutExerciseBoundary(S, K, r, q, vol, T);
442
443 const Array z = interp->nodes();
444 const Array x = 0.5*std::sqrt(T)*(1.0+z);
445
446 const auto B = [xmax, T, &interp](Real tau) -> Real {
447 const Real z = 2*std::sqrt(std::abs(tau)/T)-1;
448 return xmax*std::exp(-std::sqrt(std::max(Real(0), (*interp)(z, true))));
449 };
450
451 const auto h = [=](Real fv) -> Real {
452 return squared(std::log(fv/xmax));
453 };
454
455 const ext::shared_ptr<DqFpEquation> eqn
456 = (fpEquation_ == FP_A
457 || (fpEquation_ == Auto && std::abs(r-q) < 0.001))?
458 ext::shared_ptr<DqFpEquation>(new DqFpEquation_A(
459 K, r, q, vol, B,
460 iterationScheme_->getFixedPointIntegrator()))
461 : ext::shared_ptr<DqFpEquation>(new DqFpEquation_B(
462 K, r, q, vol, B,
463 iterationScheme_->getFixedPointIntegrator()));
464
465 Array y(x.size());
466 y[0] = 0.0;
467
468 const Size n_newton
469 = iterationScheme_->getNumberOfJacobiNewtonFixedPointSteps();
470 for (Size k=0; k < n_newton; ++k) {
471 for (Size i=1; i < x.size(); ++i) {
472 const Real tau = squared(x[i]);
473 const Real b = B(tau);
474
475 const std::tuple<Real, Real, Real> results = eqn->f(tau, b);
476 const Real N = std::get<0>(results);
477 const Real D = std::get<1>(results);
478 const Real fv = std::get<2>(results);
479
480 if (tau < QL_EPSILON)
481 y[i] = h(fv);
482 else {
483 const std::pair<Real, Real> ndd = eqn->NDd(tau, b);
484 const Real Nd = std::get<0>(ndd);
485 const Real Dd = std::get<1>(ndd);
486
487 const Real fd = K*std::exp(-(r-q)*tau) * (Nd/D - Dd*N/(D*D));
488
489 y[i] = h(b - (fv - b)/ (fd-1));
490 }
491 }
492 interp->updateY(y);
493 }
494
495 const Size n_fp = iterationScheme_->getNumberOfNaiveFixedPointSteps();
496 for (Size k=0; k < n_fp; ++k) {
497 for (Size i=1; i < x.size(); ++i) {
498 const Real tau = squared(x[i]);
499 const Real fv = std::get<2>(eqn->f(tau, B(tau)));
500
501 y[i] = h(fv);
502 }
503 interp->updateY(y);
504 }
505
506 const detail::QdPlusAddOnValue aov(T, S, K, r, q, vol, xmax, interp);
507 const Real addOn =
508 (*iterationScheme_->getExerciseBoundaryToPriceIntegrator())(
509 aov, 0.0, std::sqrt(T));
510
511 const Real europeanValue = BlackCalculator(
512 Option::Put, K, S*std::exp((r-q)*T),
513 vol*std::sqrt(T), std::exp(-r*T)).value();
514
515 return std::max(europeanValue, 0.0) + std::max(0.0, addOn);
516 }
517
518}
519
520
521
522
Black-formula calculator class.
chebyshev interpolation between discrete Chebyshev nodes
1-D array used in linear algebra.
Definition: array.hpp:52
Size size() const
dimension of the array
Definition: array.hpp:495
Black 1976 calculator class.
Integral of a one-dimensional function.
static ext::shared_ptr< QdFpIterationScheme > highPrecisionScheme()
static ext::shared_ptr< QdFpIterationScheme > fastScheme()
QdFpAmericanEngine(ext::shared_ptr< GeneralizedBlackScholesProcess > bsProcess, ext::shared_ptr< QdFpIterationScheme > iterationScheme=accurateScheme(), FixedPointEquation fpEquation=Auto)
const FixedPointEquation fpEquation_
const ext::shared_ptr< QdFpIterationScheme > iterationScheme_
static ext::shared_ptr< QdFpIterationScheme > accurateScheme()
Real calculatePut(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const override
Gauss-Legendre (l,m,n)-p Scheme.
Size getNumberOfChebyshevInterpolationNodes() const override
QdFpLegendreScheme(Size l, Size m, Size n, Size p)
ext::shared_ptr< Integrator > getFixedPointIntegrator() const override
Size getNumberOfNaiveFixedPointSteps() const override
const ext::shared_ptr< Integrator > exerciseBoundaryIntegrator_
ext::shared_ptr< Integrator > getExerciseBoundaryToPriceIntegrator() const override
Size getNumberOfJacobiNewtonFixedPointSteps() const override
const ext::shared_ptr< Integrator > fpIntegrator_
QdFpLegendreTanhSinhScheme(Size l, Size m, Size n, Real eps)
ext::shared_ptr< Integrator > getExerciseBoundaryToPriceIntegrator() const override
Size getNumberOfChebyshevInterpolationNodes() const override
QdFpTanhSinhIterationScheme(Size m, Size n, Real eps)
ext::shared_ptr< Integrator > getFixedPointIntegrator() const override
Size getNumberOfNaiveFixedPointSteps() const override
ext::shared_ptr< Integrator > getExerciseBoundaryToPriceIntegrator() const override
Size getNumberOfJacobiNewtonFixedPointSteps() const override
const ext::shared_ptr< Integrator > integrator_
American engine based on the QD+ approximation to the exercise boundary.
ext::shared_ptr< ChebyshevInterpolation > getPutExerciseBoundary(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const
static Real xMax(Real K, Rate r, Rate q)
const ext::shared_ptr< GeneralizedBlackScholesProcess > process_
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
const Matrix m_
Definition: expm.cpp:49
ext::function< Real(Real)> b
Integral of a 1-dimensional function using the Gauss quadratures.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
#define QL_MAX_REAL
Definition: qldefines.hpp:176
#define QL_EPSILON
Definition: qldefines.hpp:178
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real Volatility
volatility
Definition: types.hpp:78
QL_INTEGER Integer
integer number
Definition: types.hpp:35
Real Rate
interest rates
Definition: types.hpp:70
std::size_t Size
size of a container
Definition: types.hpp:58
functionals and combinators not included in the STL
#define M_SQRTPI
#define M_SQRT_2
#define M_SQRT2
#define M_1_SQRTPI
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:182
STL namespace.
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
Real alpha
Definition: sabr.cpp:200