QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
qdplusamericanengine.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 qrplusamericanengine.cpp
21*/
22
23#include <ql/exercise.hpp>
24#include <ql/utilities/null.hpp>
34#ifndef QL_BOOST_HAS_TANH_SINH
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(tau)), r(rf),
45 q(dy), dr(std::exp(-r * tau)), dq(std::exp(-q * tau)),
46 ddr((std::abs(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(squared(omega - 1) + 8 * ddr / sigma2))),
50 lambdaPrime(2 * ddr*ddr /
51 (sigma2 * std::sqrt(squared(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(strike, r, q)),
55 xMin(QL_EPSILON * 1e4 * std::min(0.5 * (strike + S), 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(K-S, 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, S);
98 sc = S;
99 dp = std::log(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
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
132 const Real t = z*z;
133 const Real q = (*q_z_)(2*std::sqrt(std::max(0.0, T_-t)/T_) - 1, true);
134 const Real b_t = xmax_*std::exp(-std::sqrt(std::max(0.0, q)));
135
136 const Real dr = std::exp(-r_*t);
137 const Real dq = std::exp(-q_*t);
138 const Real v = vol_*std::sqrt(t);
139
140 Real r;
141 if (v >= QL_EPSILON) {
142 if (b_t > QL_EPSILON) {
143 const Real dp = std::log(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(S_*dq, 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
161 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
162 : process_(std::move(process)) {
164 }
165
167 QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
168 "not an American option");
169
170 const auto payoff =
171 ext::dynamic_pointer_cast<StrikedTypePayoff>(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(process_->riskFreeRate()->discount(maturity))/T;
182 const Rate q = -std::log(process_->dividendYield()->discount(maturity))/T;
183 const Volatility vol = process_->blackVolatility()->blackVol(T, 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(K, S, q, r, vol, T);
193 else
194 QL_FAIL("unknown option type");
195 }
196
198 Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const {
199
200 if (close(K, 0.0))
201 return 0.0;
202
203 if (close(S, 0.0))
204 return std::max(K, K*std::exp(-r*T));
205
206 if (r <= 0.0 && r <= q)
207 return std::max(0.0,
208 BlackCalculator(Option::Put, K, S*std::exp((r-q)*T),
209 vol*std::sqrt(T), std::exp(-r*T)).value());
210
211 if (close(vol, 0.0)) {
212 const auto intrinsic = [&](Real t) -> Real {
213 return std::max(0.0, K*std::exp(-r*t)-S*std::exp(-q*t));
214 };
215 const Real npv0 = intrinsic(0.0);
216 const Real npvT = intrinsic(T);
217 const Real extremT
218 = close_enough(r, q)? QL_MAX_REAL : Real(std::log(r*K/(q*S))/(r-q));
219
220 if (extremT > 0.0 && extremT < T)
221 return std::max(npv0, std::max(npvT, intrinsic(extremT)));
222 else
223 return std::max(npv0, npvT);
224 }
225
226 return calculatePut(S, K, r, q, vol, T);
227 }
228
229
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(1.0, 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
253 ext::shared_ptr<GeneralizedBlackScholesProcess> process,
254 Size interpolationPoints,
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>
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(0.5*(eval.xmax() + S), 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(xmax, Real(-1));
284 else if (guess <= eval.xmin())
285 guess = std::nextafter(eval.xmin(), QL_MAX_REAL);
286
287 return solver.solve(eval, eps_, guess, eval.xmin(), xmax);
288 }
289
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, QuantLib::Brent(), S, K, maxIter_);
304 break;
305 case Newton:
306 x = buildInSolver(eval, QuantLib::Newton(), S, K, maxIter_);
307 break;
308 case Ridder:
309 x = buildInSolver(eval, QuantLib::Ridder(), S, K, 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(x);
323 const Real lf = fx*eval.fprime2(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(xmin, x - step);
329 resultCloseEnough = std::fabs(x-xOld) < 0.5*eps_;
330 }
331 while (!resultCloseEnough && eval.evaluations() < maxIter_);
332
333 if (!resultCloseEnough && !close(std::fabs(fx), 0.0)) {
334 x = buildInSolver(eval, QuantLib::Brent(), S, K, 10*maxIter_, 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>
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>(
353 [&, this](Real z) {
354 const Real x_sq = 0.25*T*squared(1+z);
355 return squared(std::log(
356 this->putExerciseBoundaryAtTau(S, K, r, q, vol, T, x_sq)
357 .second/xmax));
358 },
360 );
361 }
362
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(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 0.0,
388 Option::Put, K,
389 S*std::exp((r-q)*T),
390 vol*std::sqrt(T), std::exp(-r*T)).value()
391 );
392
393 return europeanValue + std::max(0.0, addOn);
394 }
395}
Black-formula calculator class.
Brent 1-D solver.
ext::shared_ptr< SimpleQuote > vol_
Definition: cdsoption.cpp:62
const Instrument::results * results_
Definition: cdsoption.cpp:63
chebyshev interpolation between discrete Chebyshev nodes
Black 1976 calculator class.
Brent 1-D solver
Definition: brent.hpp:37
Integral of a one-dimensional function.
Newton 1-D solver
Definition: newton.hpp:40
template class providing a null value for a given type.
Definition: null.hpp:76
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
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)
std::pair< Size, Real > putExerciseBoundaryAtTau(Real S, Real K, Rate r, Rate q, Volatility vol, Time T, Time tau) const
Real calculatePut(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const override
QdPlusAmericanEngine(ext::shared_ptr< GeneralizedBlackScholesProcess >, Size interpolationPoints=8, SolverType solverType=Halley, Real eps=1e-6, Size maxIter=Null< Size >())
Real buildInSolver(const QdPlusBoundaryEvaluator &eval, Solver solver, Real S, Real strike, Size maxIter, Real guess=Null< Real >()) const
Ridder 1-D solver
Definition: ridder.hpp:37
QdPlusAddOnValue(Time T, Real S, Real K, Rate r, Rate q, Volatility vol, Real xmax, ext::shared_ptr< Interpolation > q_z)
QdPutCallParityEngine(ext::shared_ptr< GeneralizedBlackScholesProcess > process)
const ext::shared_ptr< GeneralizedBlackScholesProcess > process_
Real calculatePutWithEdgeCases(Real S, Real K, Rate r, Rate q, Volatility vol, Time T) const
floating-point comparisons
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
Option exercise classes and payoff function.
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 DiscountFactor
discount factor between dates
Definition: types.hpp:66
Real Volatility
volatility
Definition: types.hpp:78
Real Rate
interest rates
Definition: types.hpp:70
std::size_t Size
size of a container
Definition: types.hpp:58
Real theta
Real sigma
ext::shared_ptr< QuantLib::Payoff > payoff
functionals and combinators not included in the STL
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
bool close(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:163
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:182
STL namespace.
Newton 1-D solver.
null values
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
Ridder 1-D solver.
Real beta
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200