QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
analytichestonengine.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) 2004, 2005, 2008 Klaus Spanderen
5 Copyright (C) 2007 StatPro Italia srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21/*! \file hestonmodel.hpp
22 \brief analytic pricing engine for a heston option
23 based on fourier transformation
24*/
25
26#include <ql/functional.hpp>
36#include <ql/math/expm1.hpp>
40
41#include <boost/math/tools/minima.hpp>
42#include <boost/math/special_functions/sign.hpp>
43
44#include <cmath>
45#include <limits>
46#include <utility>
47
48#if defined(QL_PATCH_MSVC)
49#pragma warning(disable: 4180)
50#endif
51
52namespace QuantLib {
53
54 namespace {
55
56 class integrand1 {
57 private:
58 const Real c_inf_;
59 const ext::function<Real(Real)> f_;
60 public:
61 integrand1(Real c_inf, ext::function<Real(Real)> f) : c_inf_(c_inf), f_(std::move(f)) {}
62 Real operator()(Real x) const {
63 if ((1.0-x)*c_inf_ > QL_EPSILON)
64 return f_(-std::log(0.5-0.5*x)/c_inf_)/((1.0-x)*c_inf_);
65 else
66 return 0.0;
67 }
68 };
69
70 class integrand2 {
71 private:
72 const Real c_inf_;
73 const ext::function<Real(Real)> f_;
74 public:
75 integrand2(Real c_inf, ext::function<Real(Real)> f) : c_inf_(c_inf), f_(std::move(f)) {}
76 Real operator()(Real x) const {
77 if (x*c_inf_ > QL_EPSILON) {
78 return f_(-std::log(x)/c_inf_)/(x*c_inf_);
79 } else {
80 return 0.0;
81 }
82 }
83 };
84
85 class integrand3 {
86 private:
87 const integrand2 int_;
88 public:
89 integrand3(Real c_inf, const ext::function<Real(Real)>& f)
90 : int_(c_inf, f) {}
91
92 Real operator()(Real x) const { return int_(1.0-x); }
93 };
94
95 class u_Max {
96 public:
97 u_Max(Real c_inf, Real epsilon) : c_inf_(c_inf), logEpsilon_(std::log(epsilon)) {}
98
99 Real operator()(Real u) const {
100 ++evaluations_;
101 return c_inf_*u + std::log(u) + logEpsilon_;
102 }
103
104 Size evaluations() const { return evaluations_; }
105
106 private:
108 mutable Size evaluations_ = 0;
109 };
110
111
112 class uHat_Max {
113 public:
114 uHat_Max(Real v0T2, Real epsilon) : v0T2_(v0T2), logEpsilon_(std::log(epsilon)) {}
115
116 Real operator()(Real u) const {
117 ++evaluations_;
118 return v0T2_*u*u + std::log(u) + logEpsilon_;
119 }
120
121 Size evaluations() const { return evaluations_; }
122
123 private:
125 mutable Size evaluations_ = 0;
126 };
127 }
128
129 // helper class for integration
130 class AnalyticHestonEngine::Fj_Helper {
131 public:
132 Fj_Helper(Real kappa, Real theta, Real sigma, Real v0,
133 Real s0, Real rho,
134 const AnalyticHestonEngine* engine,
135 ComplexLogFormula cpxLog,
136 Time term, Real strike, Real ratio, Size j);
137
138 Real operator()(Real phi) const;
139
140 private:
141 const Size j_;
142 // const VanillaOption::arguments& arg_;
143 const Real kappa_, theta_, sigma_, v0_;
144 const ComplexLogFormula cpxLog_;
145
146 // helper variables
147 const Time term_;
148 const Real x_, sx_, dd_;
149 const Real sigma2_, rsigma_;
150 const Real t0_;
151
152 // log branch counter
153 mutable int b_ = 0; // log branch counter
154 mutable Real g_km1_ = 0; // imag part of last log value
155
156 const AnalyticHestonEngine* const engine_;
157 };
158
159 AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa, Real theta,
160 Real sigma, Real v0, Real s0, Real rho,
161 const AnalyticHestonEngine* const engine,
162 ComplexLogFormula cpxLog,
163 Time term, Real strike, Real ratio, Size j)
164 :
165 j_(j),
166 kappa_(kappa),
167 theta_(theta),
168 sigma_(sigma),
169 v0_(v0),
170 cpxLog_(cpxLog),
171 term_(term),
172 x_(std::log(s0)),
173 sx_(std::log(strike)),
174 dd_(x_-std::log(ratio)),
175 sigma2_(sigma_*sigma_),
176 rsigma_(rho*sigma_),
177 t0_(kappa - ((j== 1)? rho*sigma : Real(0))),
178
179 engine_(engine)
180 {
181 }
182
183 Real AnalyticHestonEngine::Fj_Helper::operator()(Real phi) const
184 {
185 const Real rpsig(rsigma_*phi);
186
187 const std::complex<Real> t1 = t0_+std::complex<Real>(0, -rpsig);
188 const std::complex<Real> d =
189 std::sqrt(t1*t1 - sigma2_*phi
190 *std::complex<Real>(-phi, (j_== 1)? 1 : -1));
191 const std::complex<Real> ex = std::exp(-d*term_);
192 const std::complex<Real> addOnTerm =
193 engine_ != nullptr ? engine_->addOnTerm(phi, term_, j_) : Real(0.0);
194
195 if (cpxLog_ == Gatheral) {
196 if (phi != 0.0) {
197 if (sigma_ > 1e-5) {
198 const std::complex<Real> p = (t1-d)/(t1+d);
199 const std::complex<Real> g
200 = std::log((1.0 - p*ex)/(1.0 - p));
201
202 return
203 std::exp(v0_*(t1-d)*(1.0-ex)/(sigma2_*(1.0-ex*p))
204 + (kappa_*theta_)/sigma2_*((t1-d)*term_-2.0*g)
205 + std::complex<Real>(0.0, phi*(dd_-sx_))
206 + addOnTerm
207 ).imag()/phi;
208 }
209 else {
210 const std::complex<Real> td = phi/(2.0*t1)
211 *std::complex<Real>(-phi, (j_== 1)? 1 : -1);
212 const std::complex<Real> p = td*sigma2_/(t1+d);
213 const std::complex<Real> g = p*(1.0-ex);
214
215 return
216 std::exp(v0_*td*(1.0-ex)/(1.0-p*ex)
217 + (kappa_*theta_)*(td*term_-2.0*g/sigma2_)
218 + std::complex<Real>(0.0, phi*(dd_-sx_))
219 + addOnTerm
220 ).imag()/phi;
221 }
222 }
223 else {
224 // use l'Hospital's rule to get lim_{phi->0}
225 if (j_ == 1) {
226 const Real kmr = rsigma_-kappa_;
227 if (std::fabs(kmr) > 1e-7) {
228 return dd_-sx_
229 + (std::exp(kmr*term_)*kappa_*theta_
230 -kappa_*theta_*(kmr*term_+1.0) ) / (2*kmr*kmr)
231 - v0_*(1.0-std::exp(kmr*term_)) / (2.0*kmr);
232 }
233 else
234 // \kappa = \rho * \sigma
235 return dd_-sx_ + 0.25*kappa_*theta_*term_*term_
236 + 0.5*v0_*term_;
237 }
238 else {
239 return dd_-sx_
240 - (std::exp(-kappa_*term_)*kappa_*theta_
241 +kappa_*theta_*(kappa_*term_-1.0))/(2*kappa_*kappa_)
242 - v0_*(1.0-std::exp(-kappa_*term_))/(2*kappa_);
243 }
244 }
245 }
246 else if (cpxLog_ == BranchCorrection) {
247 const std::complex<Real> p = (t1+d)/(t1-d);
248
249 // next term: g = std::log((1.0 - p*std::exp(d*term_))/(1.0 - p))
250 std::complex<Real> g;
251
252 // the exp of the following expression is needed.
253 const std::complex<Real> e = std::log(p)+d*term_;
254
255 // does it fit to the machine precision?
256 if (std::exp(-e.real()) > QL_EPSILON) {
257 g = std::log((1.0 - p/ex)/(1.0 - p));
258 } else {
259 // use a "big phi" approximation
260 g = d*term_ + std::log(p/(p - 1.0));
261
262 if (g.imag() > M_PI || g.imag() <= -M_PI) {
263 // get back to principal branch of the complex logarithm
264 Real im = std::fmod(g.imag(), 2*M_PI);
265 if (im > M_PI)
266 im -= 2*M_PI;
267 else if (im <= -M_PI)
268 im += 2*M_PI;
269
270 g = std::complex<Real>(g.real(), im);
271 }
272 }
273
274 // be careful here as we have to use a log branch correction
275 // to deal with the discontinuities of the complex logarithm.
276 // the principal branch is not always the correct one.
277 // (s. A. Sepp, chapter 4)
278 // remark: there is still the change that we miss a branch
279 // if the order of the integration is not high enough.
280 const Real tmp = g.imag() - g_km1_;
281 if (tmp <= -M_PI)
282 ++b_;
283 else if (tmp > M_PI)
284 --b_;
285
286 g_km1_ = g.imag();
287 g += std::complex<Real>(0, 2*b_*M_PI);
288
289 return std::exp(v0_*(t1+d)*(ex-1.0)/(sigma2_*(ex-p))
290 + (kappa_*theta_)/sigma2_*((t1+d)*term_-2.0*g)
291 + std::complex<Real>(0,phi*(dd_-sx_))
292 + addOnTerm
293 ).imag()/phi;
294 }
295 else {
296 QL_FAIL("unknown complex logarithm formula");
297 }
298 }
299
300 AnalyticHestonEngine::OptimalAlpha::OptimalAlpha(
301 const Time t,
302 const AnalyticHestonEngine* const enginePtr)
303 : t_(t),
304 fwd_(enginePtr->model_->process()->s0()->value()
305 * enginePtr->model_->process()->dividendYield()->discount(t)
306 / enginePtr->model_->process()->riskFreeRate()->discount(t)),
307 kappa_(enginePtr->model_->kappa()),
308 theta_(enginePtr->model_->theta()),
309 sigma_(enginePtr->model_->sigma()),
310 rho_(enginePtr->model_->rho()),
311 eps_(std::pow(2, -int(0.5*std::numeric_limits<Real>::digits))),
312 enginePtr_(enginePtr)
313 {
314 km_ = k(0.0, -1);
315 kp_ = k(0.0, 1);
316 }
317
319 const Real eps = 1e-8;
320 const auto cm = [this](Real k) -> Real { return M(k) - t_; };
321
322 Real alpha_max;
323 const Real adx = kappa_ - sigma_*rho_;
324 if (adx > 0.0) {
325 const Real kp_2pi = k(2*M_PI, 1);
326
327 alpha_max = Brent().solve(
328 cm, eps_, 0.5*(kp_ + kp_2pi), (1+eps)*kp_, (1-eps)*kp_2pi
329 ) - 1.0;
330 }
331 else if (adx < 0.0) {
332 const Time tCut = -2/(kappa_ - sigma_*rho_*kp_);
333 if (t_ < tCut) {
334 const Real kp_pi = k(M_PI, 1);
335 alpha_max = Brent().solve(
336 cm, eps_, 0.5*(kp_ + kp_pi), (1+eps)*kp_, (1-eps)*kp_pi
337 ) - 1.0;
338 }
339 else {
340 alpha_max = Brent().solve(
341 cm, eps_, 0.5*(1.0 + kp_),
342 1 + eps, (1-eps)*kp_
343 ) - 1.0;
344 }
345 }
346 else { // adx == 0.0
347 const Real kp_pi = k(M_PI, 1);
348 alpha_max = Brent().solve(
349 cm, eps_, 0.5*(kp_ + kp_pi), (1+eps)*kp_, (1-eps)*kp_pi
350 ) - 1.0;
351 }
352
353 QL_REQUIRE(alpha_max >= 0.0, "alpha max must be larger than zero");
354
355 return alpha_max;
356 }
357
358 std::pair<Real, Real>
360 const Real alpha_max = alphaMax(strike);
361
362 return findMinima(eps_, std::max(2*eps_, (1.0-1e-6)*alpha_max), strike);
363 }
364
366 const auto cm = [this](Real k) -> Real { return M(k) - t_; };
367
368 const Real km_2pi = k(2*M_PI, -1);
369
370 const Real alpha_min = Brent().solve(
371 cm, eps_, 0.5*(km_2pi + km_), (1-1e-8)*km_2pi, (1+1e-8)*km_
372 ) - 1.0;
373
374 QL_REQUIRE(alpha_min <= -1.0, "alpha min must be smaller than minus one");
375
376 return alpha_min;
377 }
378
379 std::pair<Real, Real>
381 const Real alpha_min = alphaMin(strike);
382
383 return findMinima(
384 std::min(-1.0-1e-6, -1.0 + (1.0-1e-6)*(alpha_min + 1.0)),
385 -1.0 - eps_, strike
386 );
387 }
388
390 try {
391 const std::pair<Real, Real> minusOne = alphaSmallerMinusOne(strike);
392 const std::pair<Real, Real> greaterZero = alphaGreaterZero(strike);
393
394 if (minusOne.second < greaterZero.second) {
395 return minusOne.first;
396 }
397 else {
398 return greaterZero.first;
399 }
400 }
401 catch (const Error&) {
402 return -0.5;
403 }
404
405 }
406
408 return evaluations_;
409 }
410
412 Real lower, Real upper, Real strike) const {
413 const Real freq = std::log(fwd_/strike);
414
415 return boost::math::tools::brent_find_minima(
416 [freq, this](Real alpha) -> Real {
417 ++evaluations_;
418 const std::complex<Real> z(0, -(alpha+1));
419 return enginePtr_->lnChF(z, t_).real()
420 - std::log(alpha*(alpha+1)) + alpha*freq;
421 },
422 lower, upper, int(0.5*std::numeric_limits<Real>::digits)
423 );
424 }
425
427 const Real beta = kappa_ - sigma_*rho_*k;
428
429 if (k >= km_ && k <= kp_) {
430 const Real D = std::sqrt(beta*beta - sigma_*sigma_*k*(k-1));
431 return std::log((beta-D)/(beta+D)) / D;
432 }
433 else {
434 const Real D_imag =
435 std::sqrt(-(beta*beta - sigma_*sigma_*k*(k-1)));
436
437 return 2/D_imag
438 * ( ((beta>0.0)? M_PI : 0.0) - std::atan(D_imag/beta) );
439 }
440 }
441
443 return ( (sigma_ - 2*rho_*kappa_)
444 + sgn*std::sqrt(
445 squared(sigma_-2*rho_*kappa_)
446 + 4*(kappa_*kappa_ + x*x/(t_*t_))*(1-rho_*rho_)))
447 /(2*sigma_*(1-rho_*rho_));
448 }
449
451 Time term, Real fwd, Real strike, ComplexLogFormula cpxLog,
452 const AnalyticHestonEngine* const enginePtr,
453 const Real alpha)
454 : term_(term),
455 fwd_(fwd),
456 strike_(strike),
457 freq_(std::log(fwd/strike)),
458 cpxLog_(cpxLog),
459 enginePtr_(enginePtr),
460 alpha_(alpha),
461 s_alpha_(std::exp(alpha*freq_))
462 {
463 QL_REQUIRE(enginePtr != nullptr, "pricing engine required");
464
465 const Real v0 = enginePtr->model_->v0();
466 const Real kappa = enginePtr->model_->kappa();
467 const Real theta = enginePtr->model_->theta();
468 const Real sigma = enginePtr->model_->sigma();
469 const Real rho = enginePtr->model_->rho();
470
471 switch(cpxLog_) {
473 vAvg_ = (1-std::exp(-kappa*term))*(v0 - theta)
474 /(kappa*term) + theta;
475 break;
477 vAvg_ = -8.0*std::log(enginePtr->chF(
478 std::complex<Real>(0, alpha_), term).real())/term;
479 break;
480 case AsymptoticChF:
481 phi_ = -(v0+term*kappa*theta)/sigma
482 * std::complex<Real>(std::sqrt(1-rho*rho), rho);
483
484 psi_ = std::complex<Real>(
485 (kappa- 0.5*rho*sigma)*(v0 + term*kappa*theta)
486 + kappa*theta*std::log(4*(1-rho*rho)),
487 - ((0.5*rho*rho*sigma - kappa*rho)/std::sqrt(1-rho*rho)
488 *(v0 + kappa*theta*term)
489 - 2*kappa*theta*std::atan(rho/std::sqrt(1-rho*rho))))
490 /(sigma*sigma);
491 case AngledContour:
492 vAvg_ = (1-std::exp(-kappa*term))*(v0 - theta)
493 /(kappa*term) + theta;
495 {
496 const Real r = rho - sigma*freq_ / (v0 + kappa*theta*term);
497 tanPhi_ = std::tan(
498 (r*freq_ < 0.0)? M_PI/12*boost::math::sign(freq_) : 0.0
499 );
500 }
501 break;
502 default:
503 QL_FAIL("unknown control variate");
504 }
505 }
506
508 QL_REQUIRE( enginePtr_->addOnTerm(u, term_, 1)
509 == std::complex<Real>(0.0)
510 && enginePtr_->addOnTerm(u, term_, 2)
511 == std::complex<Real>(0.0),
512 "only Heston model is supported");
513
514 constexpr std::complex<double> i(0, 1);
515
517 const std::complex<Real> h_u(u, u*tanPhi_ - alpha_);
518 const std::complex<Real> hPrime(h_u-i);
519
520 std::complex<Real> phiBS(0.0);
521 if (cpxLog_ == AngledContour)
522 phiBS = std::exp(
523 -0.5*vAvg_*term_*(hPrime*hPrime +
524 std::complex<Real>(-hPrime.imag(), hPrime.real())));
525 else if (cpxLog_ == AsymptoticChF)
526 phiBS = std::exp(u*std::complex<Real>(1, tanPhi_)*phi_ + psi_);
527
528 return std::exp(-u*tanPhi_*freq_)
529 *(std::exp(std::complex<Real>(0.0, u*freq_))
530 *std::complex<Real>(1, tanPhi_)
531 *(phiBS - enginePtr_->chF(hPrime, term_))/(h_u*hPrime)
532 ).real()*s_alpha_;
533 }
535 const std::complex<Real> z(u, -alpha_);
536 const std::complex<Real> zPrime(u, -alpha_-1);
537 const std::complex<Real> phiBS = std::exp(
538 -0.5*vAvg_*term_*(zPrime*zPrime +
539 std::complex<Real>(-zPrime.imag(), zPrime.real()))
540 );
541
542 return (std::exp(std::complex<Real> (0.0, u*freq_))
543 * (phiBS - enginePtr_->chF(zPrime, term_)) / (z*zPrime)
544 ).real()*s_alpha_;
545 }
546 else
547 QL_FAIL("unknown control variate");
548 }
549
551 if ( cpxLog_ == AngledContour
553 return BlackCalculator(
554 Option::Call, strike_, fwd_, std::sqrt(vAvg_*term_))
555 .value();
556 }
557 else if (cpxLog_ == AsymptoticChF) {
558 QL_REQUIRE(alpha_ == -0.5, "alpha must be equal to -0.5");
559
560 const std::complex<Real> phiFreq(phi_.real(), phi_.imag() + freq_);
561
562 using namespace ExponentialIntegral;
563 return fwd_ - std::sqrt(strike_*fwd_)/M_PI*
564 (std::exp(psi_)*(
565 -2.0*Ci(-0.5*phiFreq)*std::sin(0.5*phiFreq)
566 +std::cos(0.5*phiFreq)*(M_PI+2.0*Si(0.5*phiFreq)))).real();
567 }
568 else if (cpxLog_ == AngledContourNoCV) {
569 return ((alpha_ <= 0.0)? fwd_ : 0.0)
570 - ((alpha_ <= -1.0)? strike_ : 0.0)
571 -0.5*((alpha_ == 0.0)? fwd_ :0.0)
572 +0.5*((alpha_ == -1.0)? strike_: 0.0);
573 }
574 else
575 QL_FAIL("unknown control variate");
576 }
577
578 std::complex<Real> AnalyticHestonEngine::chF(
579 const std::complex<Real>& z, Time t) const {
580 if (model_->sigma() > 1e-6 || model_->kappa() < 1e-8) {
581 return std::exp(lnChF(z, t));
582 }
583 else {
584 const Real kappa = model_->kappa();
585 const Real sigma = model_->sigma();
586 const Real theta = model_->theta();
587 const Real rho = model_->rho();
588 const Real v0 = model_->v0();
589
590 const Real sigma2 = sigma*sigma;
591
592 const Real kt = kappa*t;
593 const Real ekt = std::exp(kt);
594 const Real e2kt = std::exp(2*kt);
595 const Real rho2 = rho*rho;
596 const std::complex<Real> zpi = z + std::complex<Real>(0.0, 1.0);
597
598 return std::exp(-(((theta - v0 + ekt*((-1 + kt)*theta + v0))
599 *z*zpi)/ekt)/(2.*kappa))
600
601 + (std::exp(-(kt) - ((theta - v0 + ekt
602 *((-1 + kt)*theta + v0))*z*zpi)
603 /(2.*ekt*kappa))*rho*(2*theta + kt*theta -
604 v0 - kt*v0 + ekt*((-2 + kt)*theta + v0))
605 *(1.0 - std::complex<Real>(-z.imag(),z.real()))*z*z)
606 /(2.*kappa*kappa)*sigma
607
608 + (std::exp(-2*kt - ((theta - v0 + ekt
609 *((-1 + kt)*theta + v0))*z*zpi)/(2.*ekt*kappa))*z*z*zpi
610 *(-2*rho2*squared(2*theta + kt*theta - v0 -
611 kt*v0 + ekt*((-2 + kt)*theta + v0))
612 *z*z*zpi + 2*kappa*v0*(-zpi
613 + e2kt*(zpi + 4*rho2*z) - 2*ekt*(2*rho2*z
614 + kt*(zpi + rho2*(2 + kt)*z))) + kappa*theta*(zpi + e2kt
615 *(-5.0*zpi - 24*rho2*z+ 2*kt*(zpi + 4*rho2*z)) +
616 4*ekt*(zpi + 6*rho2*z + kt*(zpi + rho2*(4 + kt)*z)))))
617 /(16.*squared(squared(kappa)))*sigma2;
618 }
619 }
620
621 std::complex<Real> AnalyticHestonEngine::lnChF(
622 const std::complex<Real>& z, Time t) const {
623
624 const Real kappa = model_->kappa();
625 const Real sigma = model_->sigma();
626 const Real theta = model_->theta();
627 const Real rho = model_->rho();
628 const Real v0 = model_->v0();
629
630 const Real sigma2 = sigma*sigma;
631
632 const std::complex<Real> g
633 = kappa + rho*sigma*std::complex<Real>(z.imag(), -z.real());
634
635 const std::complex<Real> D = std::sqrt(
636 g*g + (z*z + std::complex<Real>(-z.imag(), z.real()))*sigma2);
637
638 // reduce cancelation errors, see. L. Andersen and M. Lake
639 std::complex<Real> r(g-D);
640 if (g.real()*D.real() + g.imag()*D.imag() > 0.0) {
641 r = -sigma2*z*std::complex<Real>(z.real(), z.imag()+1)/(g+D);
642 }
643
644 std::complex<Real> y;
645 if (D.real() != 0.0 || D.imag() != 0.0) {
646 y = expm1(-D*t)/(2.0*D);
647 }
648 else
649 y = -0.5*t;
650
651 const std::complex<Real> A
652 = kappa*theta/sigma2*(r*t - 2.0*log1p(-r*y));
653 const std::complex<Real> B
654 = z*std::complex<Real>(z.real(), z.imag()+1)*y/(1.0-r*y);
655
656 return A+v0*B;
657 }
658
660 const ext::shared_ptr<HestonModel>& model,
661 Size integrationOrder)
664 VanillaOption::results>(model),
665 evaluations_(0),
668 Integration::gaussLaguerre(integrationOrder))),
670 alpha_(-0.5) {
671 }
672
674 const ext::shared_ptr<HestonModel>& model,
675 Real relTolerance, Size maxEvaluations)
678 VanillaOption::results>(model),
679 evaluations_(0),
680 cpxLog_(OptimalCV),
681 integration_(new Integration(Integration::gaussLobatto(
682 relTolerance, Null<Real>(), maxEvaluations))),
683 andersenPiterbargEpsilon_(1e-40),
684 alpha_(-0.5) {
685 }
686
688 const ext::shared_ptr<HestonModel>& model,
689 ComplexLogFormula cpxLog,
690 const Integration& integration,
691 Real andersenPiterbargEpsilon,
692 Real alpha)
695 VanillaOption::results>(model),
696 evaluations_(0),
697 cpxLog_(cpxLog),
698 integration_(new Integration(integration)),
699 andersenPiterbargEpsilon_(andersenPiterbargEpsilon),
700 alpha_(alpha) {
702 || !integration.isAdaptiveIntegration(),
703 "Branch correction does not work in conjunction "
704 "with adaptive integration methods");
705 }
706
710
711 if (t > 0.15 && (v0+t*kappa*theta)/sigma*std::sqrt(1-rho*rho) < 0.15
712 && ((kappa- 0.5*rho*sigma)*(v0 + t*kappa*theta)
713 + kappa*theta*std::log(4*(1-rho*rho)))/(sigma*sigma) < 0.1) {
714 return AsymptoticChF;
715 }
716 else {
717 return AngledContour;
718 }
719 }
720
722 return evaluations_;
723 }
724
726 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
727 const Date& maturity) const {
728
729 const ext::shared_ptr<HestonProcess>& process = model_->process();
730 const Real fwd = process->s0()->value()
731 * process->dividendYield()->discount(maturity)
732 / process->riskFreeRate()->discount(maturity);
733
734 return priceVanillaPayoff(payoff, process->time(maturity), fwd);
735 }
736
738 const ext::shared_ptr<PlainVanillaPayoff>& payoff, Time maturity) const {
739
740 const ext::shared_ptr<HestonProcess>& process = model_->process();
741 const Real fwd = process->s0()->value()
742 * process->dividendYield()->discount(maturity)
743 / process->riskFreeRate()->discount(maturity);
744
745 return priceVanillaPayoff(payoff, maturity, fwd);
746 }
747
749 const ext::shared_ptr<PlainVanillaPayoff>& payoff,
750 Time maturity, Real fwd) const {
751
752 Real value;
753
754 const ext::shared_ptr<HestonProcess>& process = model_->process();
755 const DiscountFactor dr = process->riskFreeRate()->discount(maturity);
756
757 const Real strike = payoff->strike();
758 const Real spot = process->s0()->value();
759 QL_REQUIRE(spot > 0.0, "negative or null underlying given");
760
761 const DiscountFactor df = spot/fwd;
762 const DiscountFactor dd = dr/df;
763
764 const Real kappa = model_->kappa();
765 const Real sigma = model_->sigma();
766 const Real theta = model_->theta();
767 const Real rho = model_->rho();
768 const Real v0 = model_->v0();
769
770 evaluations_ = 0;
771
772 switch(cpxLog_) {
773 case Gatheral:
774 case BranchCorrection: {
775 const Real c_inf = std::min(0.2, std::max(0.0001,
776 std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*maturity);
777
778 const Real p1 = integration_->calculate(c_inf,
779 Fj_Helper(kappa, theta, sigma, v0, spot, rho, this,
780 cpxLog_, maturity, strike, df, 1))/M_PI;
781 evaluations_ += integration_->numberOfEvaluations();
782
783 const Real p2 = integration_->calculate(c_inf,
784 Fj_Helper(kappa, theta, sigma, v0, spot, rho, this,
785 cpxLog_, maturity, strike, df, 2))/M_PI;
786 evaluations_ += integration_->numberOfEvaluations();
787
788 switch (payoff->optionType())
789 {
790 case Option::Call:
791 value = spot*dd*(p1+0.5)
792 - strike*dr*(p2+0.5);
793 break;
794 case Option::Put:
795 value = spot*dd*(p1-0.5)
796 - strike*dr*(p2-0.5);
797 break;
798 default:
799 QL_FAIL("unknown option type");
800 }
801 }
802 break;
805 case AsymptoticChF:
806 case AngledContour:
808 case OptimalCV: {
809 const Real c_inf =
810 std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*maturity)/sigma;
811
812 const Real epsilon = andersenPiterbargEpsilon_
813 *M_PI/(std::sqrt(strike*fwd)*dr);
814
815 const ext::function<Real()> uM = [&](){
817 c_inf, epsilon, v0, maturity);
818 };
819
820 const ComplexLogFormula finalLog = (cpxLog_ == OptimalCV)
822 : cpxLog_;
823
824 const AP_Helper cvHelper(
825 maturity, fwd, strike, finalLog, this, alpha_
826 );
827
828 const Real cvValue = cvHelper.controlVariateValue();
829
830 const Real vAvg = (1-std::exp(-kappa*maturity))*(v0-theta)/(kappa*maturity) + theta;
831
832 const Real scalingFactor = (cpxLog_ != OptimalCV && cpxLog_ != AsymptoticChF)
833 ? std::max(0.25, std::min(1000.0, 0.25/std::sqrt(0.5*vAvg*maturity)))
834 : Real(1.0);
835
836 const Real h_cv =
837 fwd/M_PI*integration_->calculate(c_inf, cvHelper, uM, scalingFactor);
838
839 evaluations_ += integration_->numberOfEvaluations();
840
841 switch (payoff->optionType())
842 {
843 case Option::Call:
844 value = (cvValue + h_cv)*dr;
845 break;
846 case Option::Put:
847 value = (cvValue + h_cv - (fwd - strike))*dr;
848 break;
849 default:
850 QL_FAIL("unknown option type");
851 }
852 }
853 break;
854 default:
855 QL_FAIL("unknown complex log formula");
856 }
857
858 return value;
859 }
860
862 {
863 // this is a european option pricer
864 QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
865 "not an European option");
866
867 // plain vanilla
868 ext::shared_ptr<PlainVanillaPayoff> payoff =
869 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
870 QL_REQUIRE(payoff, "non plain vanilla payoff given");
871
872 const Date exerciseDate = arguments_.exercise->lastDate();
873
874 results_.value = priceVanillaPayoff(payoff, exerciseDate);
875 }
876
877
879 ext::shared_ptr<Integrator> integrator)
880 : intAlgo_(intAlgo), integrator_(std::move(integrator)) {}
881
883 Algorithm intAlgo, ext::shared_ptr<GaussianQuadrature> gaussianQuadrature)
884 : intAlgo_(intAlgo), gaussianQuadrature_(std::move(gaussianQuadrature)) {}
885
888 Real relTolerance, Real absTolerance, Size maxEvaluations, bool useConvergenceEstimate) {
889 return Integration(GaussLobatto,
890 ext::shared_ptr<Integrator>(
891 new GaussLobattoIntegral(maxEvaluations,
892 absTolerance,
893 relTolerance,
894 useConvergenceEstimate)));
895 }
896
899 Size maxEvaluations) {
900 return Integration(GaussKronrod,
901 ext::shared_ptr<Integrator>(
902 new GaussKronrodAdaptive(absTolerance,
903 maxEvaluations)));
904 }
905
908 Size maxEvaluations) {
909 return Integration(Simpson,
910 ext::shared_ptr<Integrator>(
911 new SimpsonIntegral(absTolerance,
912 maxEvaluations)));
913 }
914
917 Size maxEvaluations) {
918 return Integration(Trapezoid,
919 ext::shared_ptr<Integrator>(
920 new TrapezoidIntegral<Default>(absTolerance,
921 maxEvaluations)));
922 }
923
926 QL_REQUIRE(intOrder <= 192, "maximum integraton order (192) exceeded");
927 return Integration(GaussLaguerre,
928 ext::shared_ptr<GaussianQuadrature>(
929 new GaussLaguerreIntegration(intOrder)));
930 }
931
934 return Integration(GaussLegendre,
935 ext::shared_ptr<GaussianQuadrature>(
936 new GaussLegendreIntegration(intOrder)));
937 }
938
941 return Integration(GaussChebyshev,
942 ext::shared_ptr<GaussianQuadrature>(
943 new GaussChebyshevIntegration(intOrder)));
944 }
945
948 return Integration(GaussChebyshev2nd,
949 ext::shared_ptr<GaussianQuadrature>(
950 new GaussChebyshev2ndIntegration(intOrder)));
951 }
952
955 return Integration(
956 DiscreteSimpson, ext::shared_ptr<Integrator>(
957 new DiscreteSimpsonIntegrator(evaluations)));
958 }
959
962 return Integration(
963 DiscreteTrapezoid, ext::shared_ptr<Integrator>(
964 new DiscreteTrapezoidIntegrator(evaluations)));
965 }
966
969 return Integration(
970 ExpSinh, ext::shared_ptr<Integrator>(
971 new ExpSinhIntegral(relTolerance)));
972 }
973
975 if (integrator_ != nullptr) {
976 return integrator_->numberOfEvaluations();
977 } else if (gaussianQuadrature_ != nullptr) {
978 return gaussianQuadrature_->order();
979 } else {
980 QL_FAIL("neither Integrator nor GaussianQuadrature given");
981 }
982 }
983
985 return intAlgo_ == GaussLobatto
986 || intAlgo_ == GaussKronrod
987 || intAlgo_ == Simpson
988 || intAlgo_ == Trapezoid
989 || intAlgo_ == ExpSinh;
990 }
991
993 Real c_inf,
994 const ext::function<Real(Real)>& f,
995 const ext::function<Real()>& maxBound,
996 const Real scaling) const {
997
998 Real retVal;
999
1000 switch(intAlgo_) {
1001 case GaussLaguerre:
1002 retVal = (*gaussianQuadrature_)(f);
1003 break;
1004 case GaussLegendre:
1005 case GaussChebyshev:
1006 case GaussChebyshev2nd:
1007 retVal = (*gaussianQuadrature_)(integrand1(c_inf, f));
1008 break;
1009 case ExpSinh:
1010 retVal = scaling*(*integrator_)(
1011 [scaling, f](Real x) -> Real { return f(scaling*x);},
1012 0.0, std::numeric_limits<Real>::max());
1013 break;
1014 case Simpson:
1015 case Trapezoid:
1016 case GaussLobatto:
1017 case GaussKronrod:
1018 if (maxBound && maxBound() != Null<Real>())
1019 retVal = (*integrator_)(f, 0.0, maxBound());
1020 else
1021 retVal = (*integrator_)(integrand2(c_inf, f), 0.0, 1.0);
1022 break;
1023 case DiscreteTrapezoid:
1024 case DiscreteSimpson:
1025 if (maxBound && maxBound() != Null<Real>())
1026 retVal = (*integrator_)(f, 0.0, maxBound());
1027 else
1028 retVal = (*integrator_)(integrand3(c_inf, f), 0.0, 1.0);
1029 break;
1030 default:
1031 QL_FAIL("unknwon integration algorithm");
1032 }
1033
1034 return retVal;
1035 }
1036
1038 Real c_inf,
1039 const ext::function<Real(Real)>& f,
1040 Real maxBound) const {
1041
1043 c_inf, f, [=](){ return maxBound; });
1044 }
1045
1047 Real c_inf, Real epsilon, Real v0, Real t) {
1048
1049 const Real uMaxGuess = -std::log(epsilon)/c_inf;
1050 const Real uMaxStep = 0.1*uMaxGuess;
1051
1052 const Real uMax = Brent().solve(u_Max(c_inf, epsilon),
1053 QL_EPSILON*uMaxGuess, uMaxGuess, uMaxStep);
1054
1055 try {
1056 const Real uHatMaxGuess = std::sqrt(-std::log(epsilon)/(0.5*v0*t));
1057 const Real uHatMax = Brent().solve(uHat_Max(0.5*v0*t, epsilon),
1058 QL_EPSILON*uHatMaxGuess, uHatMaxGuess, 0.001*uHatMaxGuess);
1059
1060 return std::max(uMax, uHatMax);
1061 }
1062 catch (const Error&) {
1063 return uMax;
1064 }
1065 }
1066
1067
1069 Real dividendDiscount,
1070 Real spotPrice,
1071 Real strikePrice,
1072 Real term,
1074 const TypePayoff& type,
1075 const Integration& integration,
1076 const ComplexLogFormula cpxLog,
1077 const AnalyticHestonEngine* const enginePtr,
1078 Real& value,
1079 Size& evaluations) {
1080
1081 const Real ratio = riskFreeDiscount/dividendDiscount;
1082
1083 evaluations = 0;
1084
1085 switch(cpxLog) {
1086 case Gatheral:
1087 case BranchCorrection: {
1088 const Real c_inf = std::min(0.2, std::max(0.0001,
1089 std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*term);
1090
1091 const Real p1 = integration.calculate(c_inf,
1092 Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr,
1093 cpxLog, term, strikePrice, ratio, 1))/M_PI;
1094 evaluations += integration.numberOfEvaluations();
1095
1096 const Real p2 = integration.calculate(c_inf,
1097 Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr,
1098 cpxLog, term, strikePrice, ratio, 2))/M_PI;
1099 evaluations += integration.numberOfEvaluations();
1100
1101 switch (type.optionType())
1102 {
1103 case Option::Call:
1104 value = spotPrice*dividendDiscount*(p1+0.5)
1105 - strikePrice*riskFreeDiscount*(p2+0.5);
1106 break;
1107 case Option::Put:
1108 value = spotPrice*dividendDiscount*(p1-0.5)
1109 - strikePrice*riskFreeDiscount*(p2-0.5);
1110 break;
1111 default:
1112 QL_FAIL("unknown option type");
1113 }
1114 }
1115 break;
1116 case AndersenPiterbarg:
1118 case AsymptoticChF:
1119 case OptimalCV: {
1120 const Real c_inf =
1121 std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*term)/sigma;
1122
1123 const Real fwdPrice = spotPrice / ratio;
1124
1125 const Real epsilon = enginePtr->andersenPiterbargEpsilon_
1126 *M_PI/(std::sqrt(strikePrice*fwdPrice)*riskFreeDiscount);
1127
1128 const ext::function<Real()> uM = [&](){
1129 return Integration::andersenPiterbargIntegrationLimit(c_inf, epsilon, v0, term);
1130 };
1131
1132 AP_Helper cvHelper(term, fwdPrice, strikePrice,
1133 (cpxLog == OptimalCV)
1135 : cpxLog,
1136 enginePtr
1137 );
1138
1139 const Real cvValue = cvHelper.controlVariateValue();
1140
1141 const Real h_cv = integration.calculate(c_inf, cvHelper, uM)
1142 * std::sqrt(strikePrice * fwdPrice)/M_PI;
1143 evaluations += integration.numberOfEvaluations();
1144
1145 switch (type.optionType())
1146 {
1147 case Option::Call:
1148 value = (cvValue + h_cv)*riskFreeDiscount;
1149 break;
1150 case Option::Put:
1151 value = (cvValue + h_cv - (fwdPrice - strikePrice))*riskFreeDiscount;
1152 break;
1153 default:
1154 QL_FAIL("unknown option type");
1155 }
1156 }
1157 break;
1158
1159 default:
1160 QL_FAIL("unknown complex log formula");
1161 }
1162 }
1163}
const Real kappa_
const integrand2 int_
Size evaluations_
const Real v0T2_
const Real logEpsilon_
analytic Heston-model engine
Black-formula calculator class.
Brent 1-D solver.
ext::shared_ptr< PricingEngine > engine_
Definition: cdsoption.cpp:60
AP_Helper(Time term, Real fwd, Real strike, ComplexLogFormula cpxLog, const AnalyticHestonEngine *enginePtr, Real alpha=-0.5)
static Real andersenPiterbargIntegrationLimit(Real c_inf, Real epsilon, Real v0, Real t)
Integration(Algorithm intAlgo, ext::shared_ptr< GaussianQuadrature > quadrature)
static Integration discreteSimpson(Size evaluation=1000)
static Integration gaussChebyshev2nd(Size integrationOrder=128)
static Integration simpson(Real absTolerance, Size maxEvaluations=1000)
static Integration gaussLegendre(Size integrationOrder=128)
static Integration gaussLobatto(Real relTolerance, Real absTolerance, Size maxEvaluations=1000, bool useConvergenceEstimate=false)
static Integration gaussChebyshev(Size integrationOrder=128)
static Integration expSinh(Real relTolerance=1e-8)
static Integration trapezoid(Real absTolerance, Size maxEvaluations=1000)
static Integration gaussLaguerre(Size integrationOrder=128)
static Integration gaussKronrod(Real absTolerance, Size maxEvaluations=1000)
Real calculate(Real c_inf, const ext::function< Real(Real)> &f, const ext::function< Real()> &maxBound={}, Real scaling=1.0) const
static Integration discreteTrapezoid(Size evaluation=1000)
std::pair< Real, Real > alphaGreaterZero(Real strike) const
std::pair< Real, Real > alphaSmallerMinusOne(Real strike) const
std::pair< Real, Real > findMinima(Real lower, Real upper, Real strike) const
analytic Heston-model engine based on Fourier transform
std::complex< Real > lnChF(const std::complex< Real > &z, Time t) const
static void doCalculation(Real riskFreeDiscount, Real dividendDiscount, Real spotPrice, Real strikePrice, Real term, Real kappa, Real theta, Real sigma, Real v0, Real rho, const TypePayoff &type, const Integration &integration, ComplexLogFormula cpxLog, const AnalyticHestonEngine *enginePtr, Real &value, Size &evaluations)
Real priceVanillaPayoff(const ext::shared_ptr< PlainVanillaPayoff > &payoff, const Date &maturity) const
AnalyticHestonEngine(const ext::shared_ptr< HestonModel > &model, Real relTolerance, Size maxEvaluations)
const ext::shared_ptr< Integration > integration_
std::complex< Real > chF(const std::complex< Real > &z, Time t) const
static ComplexLogFormula optimalControlVariate(Time t, Real v0, Real kappa, Real theta, Real sigma, Real rho)
Black 1976 calculator class.
Brent 1-D solver
Definition: brent.hpp:37
Concrete date class.
Definition: date.hpp:125
Base error class.
Definition: errors.hpp:39
Gauss-Chebyshev integration (second kind)
Integral of a 1-dimensional function using the Gauss-Kronrod methods.
generalized Gauss-Laguerre integration
Integral of a one-dimensional function.
Base class for some pricing engine on a particular model.
Heston model for the stochastic volatility of an asset.
Definition: hestonmodel.hpp:42
template class providing a null value for a given type.
Definition: null.hpp:76
Integral of a one-dimensional function.
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
Integral of a one-dimensional function.
Intermediate class for put/call payoffs.
Definition: payoffs.hpp:49
Option::Type optionType() const
Definition: payoffs.hpp:51
Vanilla option (no discrete dividends, no barriers) on a single asset.
ext::function< Real(Real)> f_
const DefaultType & t
integrals on non uniform grids
#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
complex versions of expm1 and logp1
Maps function, bind and cref to either the boost or std implementation.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
#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
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
Real kappa
Real theta
const Real c_inf_
Real v0
Real rho
Real sigma
ext::shared_ptr< QuantLib::Payoff > payoff
Integral of a 1-dimensional function using the Gauss-Kronrod method.
const VF_R b_
functionals and combinators not included in the STL
#define M_PI
Definition: any.hpp:35
std::complex< Real > log1p(const std::complex< Real > &z)
Definition: expm1.cpp:42
T squared(T x)
Definition: functional.hpp:37
std::complex< Real > expm1(const std::complex< Real > &z)
Definition: expm1.cpp:26
STL namespace.
Payoffs for various options.
ext::shared_ptr< YieldTermStructure > r
Real beta
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
integral of a one-dimensional function using Simpson formula
integral of a one-dimensional function using the trapezoid formula
std::uint64_t x_