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 qrplusamericanengine.cpp
21*/
22
23#include <ql/exercise.hpp>
24#include <ql/utilities/null.hpp>
25#include <ql/math/functional.hpp>
26#include <ql/math/comparison.hpp>
27#include <ql/math/solvers1d/brent.hpp>
28#include <ql/math/solvers1d/ridder.hpp>
29#include <ql/math/solvers1d/newton.hpp>
30#include <ql/math/interpolations/chebyshevinterpolation.hpp>
31#include <ql/pricingengines/blackcalculator.hpp>
32#include <ql/pricingengines/vanilla/qdplusamericanengine.hpp>
33#include <ql/math/integrals/tanhsinhintegral.hpp>
34#ifndef QL_BOOST_HAS_TANH_SINH
35#include <ql/math/integrals/gausslobattointegral.hpp>
36#endif
37
38namespace QuantLib {
39
40 class QdPlusBoundaryEvaluator {
41 public:
42 QdPlusBoundaryEvaluator(
43 Real S, Real strike, Rate rf, Rate dy, Volatility vol, Time t, Time T)
44 : tau(t), K(strike), sigma(vol), sigma2(sigma * sigma), v(sigma * std::sqrt(x: tau)), r(rf),
45 q(dy), dr(std::exp(x: -r * tau)), dq(std::exp(x: -q * tau)),
46 ddr((std::abs(x: r*tau) > 1e-5)? Real(r/(1-dr)) : Real(1/(tau*(1-0.5*r*tau*(1-r*tau/3))))),
47 omega(2 * (r - q) / sigma2),
48 lambda(0.5 *
49 (-(omega - 1) - std::sqrt(x: squared(x: omega - 1) + 8 * ddr / sigma2))),
50 lambdaPrime(2 * ddr*ddr /
51 (sigma2 * std::sqrt(x: squared(x: omega - 1) + 8 * ddr / sigma2))),
52 alpha(2 * dr / (sigma2 * (2 * lambda + omega - 1))),
53 beta(alpha * (ddr + lambdaPrime / (2 * lambda + omega - 1)) - lambda),
54 xMax(QdPlusAmericanEngine::xMax(K: strike, r, q)),
55 xMin(QL_EPSILON * 1e4 * std::min(a: 0.5 * (strike + S), b: xMax)),
56
57 sc(Null<Real>()) {}
58
59 Real operator()(Real S) const {
60 ++nrEvaluations;
61
62 if (S != sc)
63 preCalculate(S);
64
65 if (close_enough(x: K-S, y: npv)) {
66 return (1-dq*Phi_dp)*S + alpha*theta/dr;
67 }
68 else {
69 const Real c0 = -beta - lambda + alpha*theta/(dr*(K-S-npv));
70 return (1-dq*Phi_dp)*S + (lambda+c0)*(K-S-npv);
71 }
72 }
73 Real derivative(Real S) const {
74 if (S != sc)
75 preCalculate(S);
76
77 return 1 - dq*Phi_dp + dq/v*phi_dp + beta*(1-dq*Phi_dp)
78 + alpha/dr*charm;
79 }
80 Real fprime2(Real S) const {
81 if (S != sc)
82 preCalculate(S);
83
84 const Real gamma = phi_dp*dq/(v*S);
85 const Real colour = gamma*(q + (r-q)*dp/v + (1-dp*dm)/(2*tau));
86
87 return dq*(phi_dp/(S*v) - phi_dp*dp/(S*v*v))
88 + beta*gamma + alpha/dr*colour;
89 }
90
91 Real xmin() const { return xMin; }
92 Real xmax() const { return xMax; }
93 Size evaluations() const { return nrEvaluations; }
94
95 private:
96 void preCalculate(Real S) const {
97 S = std::max(QL_EPSILON, b: S);
98 sc = S;
99 dp = std::log(x: S*dq/(K*dr))/v + 0.5*v;
100 dm = dp - v;
101 Phi_dp = Phi(-dp);
102 Phi_dm = Phi(-dm);
103 phi_dp = phi(dp);
104
105 npv = dr*K*Phi_dm - S*dq*Phi_dp;
106 theta = r*K*dr*Phi_dm - q*S*dq*Phi_dp - sigma2*S/(2*v)*dq*phi_dp;
107 charm = -dq*(phi_dp*((r-q)/v - dm/(2*tau)) +q*Phi_dp);
108 }
109
110 const CumulativeNormalDistribution Phi;
111 const NormalDistribution phi;
112 const Time tau;
113 const Real K;
114 const Volatility sigma, sigma2, v;
115 const Rate r, q;
116 const DiscountFactor dr, dq, ddr;
117 const Real omega, lambda, lambdaPrime, alpha, beta, xMax, xMin;
118 mutable Size nrEvaluations = 0;
119 mutable Real sc, dp, dm, Phi_dp, Phi_dm, phi_dp;
120 mutable Real npv, theta, charm;
121 };
122
123
124 detail::QdPlusAddOnValue::QdPlusAddOnValue(
125 Time T, Real S, Real K, Rate r, Rate q, Volatility vol,
126 const Real xmax, ext::shared_ptr<Interpolation> q_z)
127 : T_(T), S_(S), K_(K), xmax_(xmax),
128 r_(r), q_(q), vol_(vol), q_z_(std::move(q_z)) {}
129
130
131 Real detail::QdPlusAddOnValue::operator()(Real z) const {
132 const Real t = z*z;
133 const Real q = (*q_z_)(2*std::sqrt(x: std::max(a: 0.0, b: T_-t)/T_) - 1, true);
134 const Real b_t = xmax_*std::exp(x: -std::sqrt(x: std::max(a: 0.0, b: q)));
135
136 const Real dr = std::exp(x: -r_*t);
137 const Real dq = std::exp(x: -q_*t);
138 const Real v = vol_*std::sqrt(x: t);
139
140 Real r;
141 if (v >= QL_EPSILON) {
142 if (b_t > QL_EPSILON) {
143 const Real dp = std::log(x: S_*dq/(b_t*dr))/v + 0.5*v;
144 r = 2*z*(r_*K_*dr*Phi_(-dp+v) - q_*S_*dq*Phi_(-dp));
145 }
146 else
147 r = 0.0;
148 }
149 else if (close_enough(x: S_*dq, y: b_t*dr))
150 r = z*(r_*K_*dr - q_*S_*dq);
151 else if (b_t*dr > S_*dq)
152 r = 2*z*(r_*K_*dr - q_*S_*dq);
153 else
154 r = 0.0;
155
156 return r;
157 }
158
159
160 detail::QdPutCallParityEngine::QdPutCallParityEngine(
161 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
162 : process_(std::move(process)) {
163 registerWith(h: process_);
164 }
165
166 void detail::QdPutCallParityEngine::calculate() const {
167 QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
168 "not an American option");
169
170 const auto payoff =
171 ext::dynamic_pointer_cast<StrikedTypePayoff>(r: arguments_.payoff);
172 QL_REQUIRE(payoff, "non-striked payoff given");
173
174 const Real spot = process_->x0();
175 QL_REQUIRE(spot >= 0.0, "negative underlying given");
176
177 const auto maturity = arguments_.exercise->lastDate();
178 const Time T = process_->time(maturity);
179 const Real S = process_->x0();
180 const Real K = payoff->strike();
181 const Rate r = -std::log(x: process_->riskFreeRate()->discount(d: maturity))/T;
182 const Rate q = -std::log(x: process_->dividendYield()->discount(d: maturity))/T;
183 const Volatility vol = process_->blackVolatility()->blackVol(t: T, strike: K);
184
185 QL_REQUIRE(S >= 0, "zero or positive underlying value is required");
186 QL_REQUIRE(K >= 0, "zero or positive strike is required");
187 QL_REQUIRE(vol >= 0, "zero or positive volatility is required");
188
189 if (payoff->optionType() == Option::Put)
190 results_.value = calculatePutWithEdgeCases(S, K, r, q, vol, T);
191 else if (payoff->optionType() == Option::Call)
192 results_.value = calculatePutWithEdgeCases(S: K, K: S, r: q, q: r, vol, T);
193 else
194 QL_FAIL("unknown option type");
195 }
196
197 Real detail::QdPutCallParityEngine::calculatePutWithEdgeCases(
198 Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
199
200 if (close(x: K, y: 0.0))
201 return 0.0;
202
203 if (close(x: S, y: 0.0))
204 return std::max(a: K, b: K*std::exp(x: -r*T));
205
206 if (r <= 0.0 && r <= q)
207 return std::max(a: 0.0,
208 b: BlackCalculator(Option::Put, K, S*std::exp(x: (r-q)*T),
209 vol*std::sqrt(x: T), std::exp(x: -r*T)).value());
210
211 if (close(x: vol, y: 0.0)) {
212 const auto intrinsic = [&](Real t) -> Real {
213 return std::max(a: 0.0, b: K*std::exp(x: -r*t)-S*std::exp(x: -q*t));
214 };
215 const Real npv0 = intrinsic(0.0);
216 const Real npvT = intrinsic(T);
217 const Real extremT
218 = close_enough(x: r, y: q)? QL_MAX_REAL : Real(std::log(x: r*K/(q*S))/(r-q));
219
220 if (extremT > 0.0 && extremT < T)
221 return std::max(a: npv0, b: std::max(a: npvT, b: intrinsic(extremT)));
222 else
223 return std::max(a: npv0, b: npvT);
224 }
225
226 return calculatePut(S, K, r, q, vol, T);
227 }
228
229
230 Real QdPlusAmericanEngine::xMax(Real K, Rate r, Rate q) {
231 //Table 2 from Leif Andersen, Mark Lake (2021)
232 //"Fast American Option Pricing: The Double-Boundary Case"
233
234 if (r > 0.0 && q > 0.0)
235 return K*std::min(a: 1.0, b: r/q);
236 else if (r > 0.0 && q <= 0.0)
237 return K;
238 else if (r == 0.0 && q < 0.0)
239 return K;
240 else if (r == 0.0 && q >= 0.0)
241 return 0.0; // Eurpoean case
242 else if (r < 0.0 && q >= 0.0)
243 return 0.0; // European case
244 else if (r < 0.0 && q < r)
245 return K; // double boundary case
246 else if (r < 0.0 && r <= q && q < 0.0)
247 return 0; // European case
248 else
249 QL_FAIL("internal error");
250 }
251
252 QdPlusAmericanEngine::QdPlusAmericanEngine(
253 ext::shared_ptr<GeneralizedBlackScholesProcess> process,
254 Size interpolationPoints,
255 QdPlusAmericanEngine::SolverType solverType,
256 Real eps, Size maxIter)
257 : detail::QdPutCallParityEngine(std::move(process)),
258 interpolationPoints_(interpolationPoints),
259 solverType_(solverType),
260 eps_(eps),
261 maxIter_((maxIter == Null<Size>()) ? (
262 (solverType == Newton
263 || solverType == Brent || solverType== Ridder)? 100 : 10)
264 : maxIter ) { }
265
266 template <class Solver>
267 Real QdPlusAmericanEngine::buildInSolver(
268 const QdPlusBoundaryEvaluator& eval,
269 Solver solver, Real S, Real strike, Size maxIter, Real guess) const {
270
271 solver.setMaxEvaluations(maxIter);
272 solver.setLowerBound(eval.xmin());
273
274 const Real fxmin = eval(eval.xmin());
275 Real xmax = std::max(a: 0.5*(eval.xmax() + S), b: eval.xmax());
276 while (eval(xmax)*fxmin > 0.0 && eval.evaluations() < maxIter_)
277 xmax*=2;
278
279 if (guess == Null<Real>())
280 guess = 0.5*(xmax + S);
281
282 if (guess >= xmax)
283 guess = std::nextafter(x: xmax, y: Real(-1));
284 else if (guess <= eval.xmin())
285 guess = std::nextafter(x: eval.xmin(), QL_MAX_REAL);
286
287 return solver.solve(eval, eps_, guess, eval.xmin(), xmax);
288 }
289
290 std::pair<Size, Real> QdPlusAmericanEngine::putExerciseBoundaryAtTau(
291 Real S, Real K, Rate r, Rate q,
292 Volatility vol, Time T, Time tau) const {
293
294 if (tau < QL_EPSILON)
295 return std::pair<Size, Real>(
296 Size(0), xMax(K, r, q));
297
298 const QdPlusBoundaryEvaluator eval(S, K, r, q, vol, tau, T);
299
300 Real x;
301 switch (solverType_) {
302 case Brent:
303 x = buildInSolver(eval, solver: QuantLib::Brent(), S, strike: K, maxIter: maxIter_);
304 break;
305 case Newton:
306 x = buildInSolver(eval, solver: QuantLib::Newton(), S, strike: K, maxIter: maxIter_);
307 break;
308 case Ridder:
309 x = buildInSolver(eval, solver: QuantLib::Ridder(), S, strike: K, maxIter: maxIter_);
310 break;
311 case Halley:
312 case SuperHalley:
313 {
314 bool resultCloseEnough;
315 x = eval.xmax();
316 Real xOld, fx;
317 const Real xmin = eval.xmin();
318
319 do {
320 xOld = x;
321 fx = eval(x);
322 const Real fPrime = eval.derivative(S: x);
323 const Real lf = fx*eval.fprime2(S: x)/(fPrime*fPrime);
324 const Real step = (solverType_ == Halley)
325 ? Real(1/(1 - 0.5*lf)*fx/fPrime)
326 : Real((1 + 0.5*lf/(1-lf))*fx/fPrime);
327
328 x = std::max(a: xmin, b: x - step);
329 resultCloseEnough = std::fabs(x: x-xOld) < 0.5*eps_;
330 }
331 while (!resultCloseEnough && eval.evaluations() < maxIter_);
332
333 if (!resultCloseEnough && !close(x: std::fabs(x: fx), y: 0.0)) {
334 x = buildInSolver(eval, solver: QuantLib::Brent(), S, strike: K, maxIter: 10*maxIter_, guess: x);
335 }
336 }
337 break;
338 default:
339 QL_FAIL("unknown solver type");
340 }
341
342 return std::pair<Size, Real>(eval.evaluations(), x);
343 }
344
345 ext::shared_ptr<ChebyshevInterpolation>
346 QdPlusAmericanEngine::getPutExerciseBoundary(
347 Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
348
349 const Real xmax = xMax(K, r, q);
350
351 return ext::make_shared<ChebyshevInterpolation>(
352 args: interpolationPoints_,
353 args: [&, this](Real z) {
354 const Real x_sq = 0.25*T*squared(x: 1+z);
355 return squared(x: std::log(
356 x: this->putExerciseBoundaryAtTau(S, K, r, q, vol, T, tau: x_sq)
357 .second/xmax));
358 },
359 args: ChebyshevInterpolation::SecondKind
360 );
361 }
362
363 Real QdPlusAmericanEngine::calculatePut(
364 Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
365
366 if (r < 0.0 && q < r)
367 QL_FAIL("double-boundary case q<r<0 for a put option is given");
368
369 const ext::shared_ptr<Interpolation> q_z
370 = getPutExerciseBoundary(S, K, r, q, vol, T);
371
372 const Real xmax = xMax(K, r, q);
373 const detail::QdPlusAddOnValue aov(T, S, K, r, q, vol, xmax, q_z);
374
375#ifdef QL_BOOST_HAS_TANH_SINH
376 const Real addOn = TanhSinhIntegral(eps_)(aov, 0.0, std::sqrt(x: T));
377#else
378 const Real addOn = GaussLobattoIntegral(100000, QL_MAX_REAL, 0.1*eps_)
379 (aov, 0.0, std::sqrt(T));
380#endif
381
382 QL_REQUIRE(addOn > -10*eps_,
383 "negative early exercise value " << addOn);
384
385 const Real europeanValue = std::max(
386 a: 0.0,
387 b: BlackCalculator(
388 Option::Put, K,
389 S*std::exp(x: (r-q)*T),
390 vol*std::sqrt(x: T), std::exp(x: -r*T)).value()
391 );
392
393 return europeanValue + std::max(a: 0.0, b: addOn);
394 }
395}
396

source code of quantlib/ql/pricingengines/vanilla/qdplusamericanengine.cpp