26#include <ql/functional.hpp>
27#include <ql/instruments/payoffs.hpp>
28#include <ql/math/integrals/discreteintegrals.hpp>
29#include <ql/math/integrals/exponentialintegrals.hpp>
30#include <ql/math/integrals/gausslobattointegral.hpp>
31#include <ql/math/integrals/kronrodintegral.hpp>
32#include <ql/math/integrals/simpsonintegral.hpp>
33#include <ql/math/integrals/trapezoidintegral.hpp>
34#include <ql/math/solvers1d/brent.hpp>
35#include <ql/math/functional.hpp>
36#include <ql/pricingengines/blackcalculator.hpp>
37#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
40#if defined(QL_PATCH_MSVC)
41#pragma warning(disable: 4180)
53 integrand1(
Real c_inf, ext::function<
Real(
Real)> f) : c_inf_(c_inf), f_(
std::move(f)) {}
56 return f_(-std::log(0.5-0.5*x)/c_inf_)/((1.0-x)*c_inf_);
67 integrand2(
Real c_inf, ext::function<
Real(
Real)> f) : c_inf_(c_inf), f_(
std::move(f)) {}
70 return f_(-std::log(x)/c_inf_)/(x*c_inf_);
79 const integrand2 int_;
81 integrand3(
Real c_inf,
const ext::function<
Real(
Real)>& f)
84 Real operator()(
Real x)
const {
return int_(1.0-x); }
89 u_Max(
Real c_inf,
Real epsilon) : c_inf_(c_inf), logEpsilon_(
std::log(epsilon)) {}
93 return c_inf_*u + std::log(u) + logEpsilon_;
96 Size evaluations()
const {
return evaluations_; }
99 const Real c_inf_, logEpsilon_;
100 mutable Size evaluations_ = 0;
106 uHat_Max(
Real v0T2,
Real epsilon) : v0T2_(v0T2), logEpsilon_(
std::log(epsilon)) {}
110 return v0T2_*u*u + std::log(u) + logEpsilon_;
113 Size evaluations()
const {
return evaluations_; }
116 const Real v0T2_, logEpsilon_;
117 mutable Size evaluations_ = 0;
122 class AnalyticHestonEngine::Fj_Helper {
124 Fj_Helper(
const VanillaOption::arguments& arguments,
125 const ext::shared_ptr<HestonModel>& model,
126 const AnalyticHestonEngine* engine,
132 Fj_Helper(
Real kappa,
138 const AnalyticHestonEngine* engine,
145 Fj_Helper(
Real kappa,
162 const Real kappa_, theta_, sigma_, v0_;
167 const Real x_, sx_, dd_;
168 const Real sigma2_, rsigma_;
175 const AnalyticHestonEngine*
const engine_;
179 AnalyticHestonEngine::Fj_Helper::Fj_Helper(
180 const VanillaOption::arguments& arguments,
181 const ext::shared_ptr<HestonModel>& model,
182 const AnalyticHestonEngine*
const engine,
183 ComplexLogFormula cpxLog,
186 kappa_(model->kappa()), theta_(model->theta()),
187 sigma_(model->sigma()), v0_(model->v0()),
188 cpxLog_(cpxLog), term_(term),
189 x_(
std::log(model->process()->s0()->value())),
190 sx_(
std::log(ext::dynamic_pointer_cast<StrikedTypePayoff>
191 (arguments.payoff)->strike())),
192 dd_(x_-
std::log(ratio)),
193 sigma2_(sigma_*sigma_),
194 rsigma_(model->rho()*sigma_),
195 t0_(kappa_ - ((j_== 1)? model->rho()*sigma_ :
Real(0))),
201 AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa, Real theta,
202 Real sigma, Real v0, Real s0, Real rho,
203 const AnalyticHestonEngine*
const engine,
204 ComplexLogFormula cpxLog,
218 sx_(
std::log(strike)),
219 dd_(x_-
std::log(ratio)),
220 sigma2_(sigma_*sigma_),
222 t0_(kappa - ((j== 1)? rho*sigma :
Real(0))),
229 AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa,
235 ComplexLogFormula cpxLog,
240 : j_(j), kappa_(kappa), theta_(theta), sigma_(sigma), v0_(v0), cpxLog_(cpxLog), term_(term),
241 x_(
std::log(s0)), sx_(
std::log(strike)), dd_(x_ -
std::log(ratio)), sigma2_(sigma_ * sigma_),
242 rsigma_(rho * sigma_), t0_(kappa - ((j == 1) ? rho * sigma :
Real(0))), b_(0), g_km1_(0),
246 Real AnalyticHestonEngine::Fj_Helper::operator()(Real phi)
const
248 const Real rpsig(rsigma_*phi);
250 const std::complex<Real> t1 = t0_+std::complex<Real>(0, -rpsig);
251 const std::complex<Real> d =
252 std::sqrt(t1*t1 - sigma2_*phi
253 *std::complex<Real>(-phi, (j_== 1)? 1 : -1));
254 const std::complex<Real> ex = std::exp(-d*term_);
255 const std::complex<Real> addOnTerm =
256 engine_ !=
nullptr ? engine_->addOnTerm(phi, term_, j_) :
Real(0.0);
258 if (cpxLog_ == Gatheral) {
261 const std::complex<Real> p = (t1-d)/(t1+d);
262 const std::complex<Real>
g
263 = std::log((1.0 - p*ex)/(1.0 - p));
266 std::exp(v0_*(t1-d)*(1.0-ex)/(sigma2_*(1.0-ex*p))
267 + (kappa_*theta_)/sigma2_*((t1-d)*term_-2.0*g)
268 + std::complex<Real>(0.0, phi*(dd_-sx_))
273 const std::complex<Real> td = phi/(2.0*t1)
274 *std::complex<Real>(-phi, (j_== 1)? 1 : -1);
275 const std::complex<Real> p = td*sigma2_/(t1+d);
276 const std::complex<Real>
g = p*(1.0-ex);
279 std::exp(v0_*td*(1.0-ex)/(1.0-p*ex)
280 + (kappa_*theta_)*(td*term_-2.0*g/sigma2_)
281 + std::complex<Real>(0.0, phi*(dd_-sx_))
289 const Real kmr = rsigma_-kappa_;
290 if (std::fabs(kmr) > 1e-7) {
292 + (std::exp(kmr*term_)*kappa_*theta_
293 -kappa_*theta_*(kmr*term_+1.0) ) / (2*kmr*kmr)
294 - v0_*(1.0-std::exp(kmr*term_)) / (2.0*kmr);
298 return dd_-sx_ + 0.25*kappa_*theta_*term_*term_
303 - (std::exp(-kappa_*term_)*kappa_*theta_
304 +kappa_*theta_*(kappa_*term_-1.0))/(2*kappa_*kappa_)
305 - v0_*(1.0-std::exp(-kappa_*term_))/(2*kappa_);
309 else if (cpxLog_ == BranchCorrection) {
310 const std::complex<Real> p = (t1+d)/(t1-d);
313 std::complex<Real>
g;
316 const std::complex<Real> e = std::log(p)+d*term_;
320 g = std::log((1.0 - p/ex)/(1.0 - p));
323 g = d*term_ + std::log(p/(p - 1.0));
325 if (
g.imag() > M_PI ||
g.imag() <= -M_PI) {
327 Real im = std::fmod(
g.imag(), 2*M_PI);
330 else if (im <= -M_PI)
333 g = std::complex<Real>(
g.real(), im);
343 const Real tmp =
g.imag() - g_km1_;
350 g += std::complex<Real>(0, 2*b_*M_PI);
352 return std::exp(v0_*(t1+d)*(ex-1.0)/(sigma2_*(ex-p))
353 + (kappa_*theta_)/sigma2_*((t1+d)*term_-2.0*g)
354 + std::complex<Real>(0,phi*(dd_-sx_))
359 QL_FAIL(
"unknown complex logarithm formula");
364 AnalyticHestonEngine::AP_Helper::AP_Helper(
370 freq_(
std::log(fwd/strike)),
372 enginePtr_(enginePtr) {
373 QL_REQUIRE(enginePtr !=
nullptr,
"pricing engine required");
376 const Real kappa = enginePtr->
model_->kappa();
377 const Real theta = enginePtr->
model_->theta();
378 const Real sigma = enginePtr->
model_->sigma();
383 vAvg_ = (1-std::exp(-kappa*term))*(v0 - theta)
384 /(kappa*term) + theta;
387 vAvg_ = -8.0*std::log(enginePtr->
chF(
388 std::complex<Real>(0, -0.5), term).real())/term;
391 phi_ = -(v0+term*kappa*theta)/sigma
392 * std::complex<Real>(std::sqrt(1-rho*rho), rho);
394 psi_ = std::complex<Real>(
395 (kappa- 0.5*rho*sigma)*(v0 + term*kappa*theta)
396 + kappa*theta*std::log(4*(1-rho*rho)),
397 - ((0.5*rho*rho*sigma - kappa*rho)/std::sqrt(1-rho*rho)
398 *(v0 + kappa*theta*term)
399 - 2*kappa*theta*std::atan(rho/std::sqrt(1-rho*rho))))
403 QL_FAIL(
"unknown control variate");
408 QL_REQUIRE( enginePtr_->addOnTerm(u, term_, 1)
409 == std::complex<Real>(0.0)
410 && enginePtr_->addOnTerm(u, term_, 2)
411 == std::complex<Real>(0.0),
412 "only Heston model is supported");
414 const std::complex<Real> z(u, -0.5);
416 std::complex<Real> phiBS;
422 -0.5*vAvg_*term_*(z*z + std::complex<Real>(-z.imag(), z.real())));
425 phiBS = std::exp(u*phi_ + psi_);
428 QL_FAIL(
"unknown control variate");
431 return (std::exp(std::complex<Real>(0.0, u*freq_))
432 * (phiBS - enginePtr_->chF(z, term_)) / (u*u + 0.25)).real();
442 const std::complex<Real> phiFreq(phi_.real(), phi_.imag() + freq_);
444 using namespace ExponentialIntegral;
445 return fwd_ - std::sqrt(strike_*fwd_)/M_PI*
447 -2.0*Ci(-0.5*phiFreq)*std::sin(0.5*phiFreq)
448 +std::cos(0.5*phiFreq)*(M_PI+2.0*Si(0.5*phiFreq)))).real();
451 QL_FAIL(
"unknown control variate");
455 const std::complex<Real>& z,
Time t)
const {
463 const Real sigma2 = sigma*sigma;
466 const std::complex<Real> g
467 = kappa + rho*sigma*std::complex<Real>(z.imag(), -z.real());
469 const std::complex<Real> D = std::sqrt(
470 g*g + (z*z + std::complex<Real>(-z.imag(), z.real()))*sigma2);
472 const std::complex<Real> G = (g-D)/(g+D);
474 return std::exp(v0/sigma2*(1.0-std::exp(-D*t))/(1.0-G*std::exp(-D*t))
475 *(g-D) + kappa*theta/sigma2*((g-D)*t
476 -2.0*std::log((1.0-G*std::exp(-D*t))/(1.0-G))));
479 const Real kt = kappa*t;
480 const Real ekt = std::exp(kt);
481 const Real e2kt = std::exp(2*kt);
482 const Real rho2 = rho*rho;
483 const std::complex<Real> zpi = z + std::complex<Real>(0.0, 1.0);
485 return std::exp(-(((theta - v0 + ekt*((-1 + kt)*theta + v0))
486 *z*zpi)/ekt)/(2.*kappa))
487 + (std::exp(-(kt) - ((theta - v0 + ekt
488 *((-1 + kt)*theta + v0))*z*zpi)
489 /(2.*ekt*kappa))*rho*(2*theta + kt*theta -
490 v0 - kt*v0 + ekt*((-2 + kt)*theta + v0))
491 *(1.0 - std::complex<Real>(-z.imag(),z.real()))*z*z)
492 /(2.*kappa*kappa)*sigma
493 + (std::exp(-2*kt - ((theta - v0 + ekt
494 *((-1 + kt)*theta + v0))*z*zpi)/(2.*ekt*kappa))*z*z*zpi
495 *(-2*rho2*
squared(2*theta + kt*theta - v0 -
496 kt*v0 + ekt*((-2 + kt)*theta + v0))
497 *z*z*zpi + 2*kappa*v0*(-zpi
498 + e2kt*(zpi + 4*rho2*z) - 2*ekt*(2*rho2*z
499 + kt*(zpi + rho2*(2 + kt)*z))) + kappa*theta*(zpi + e2kt
500 *(-5.0*zpi - 24*rho2*z+ 2*kt*(zpi + 4*rho2*z)) +
501 4*ekt*(zpi + 6*rho2*z + kt*(zpi + rho2*(4 + kt)*z)))))
507 const std::complex<Real>& z,
Time T)
const {
508 return std::log(
chF(z, T));
512 const ext::shared_ptr<HestonModel>& model,
513 Size integrationOrder)
525 const ext::shared_ptr<HestonModel>& model,
526 Real relTolerance,
Size maxEvaluations)
533 relTolerance,
Null<
Real>(), maxEvaluations))),
534 andersenPiterbargEpsilon_(
Null<
Real>()) {
538 const ext::shared_ptr<HestonModel>& model,
541 const Real andersenPiterbargEpsilon)
548 andersenPiterbargEpsilon_(andersenPiterbargEpsilon) {
551 "Branch correction does not work in conjunction "
552 "with adaptive integration methods");
559 if (t > 0.1 && (v0+t*kappa*theta)/sigma*std::sqrt(1-rho*rho) < 0.055) {
573 Real dividendDiscount,
585 const Real ratio = riskFreeDiscount/dividendDiscount;
592 const Real c_inf = std::min(0.2, std::max(0.0001,
593 std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*term);
596 Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr,
597 cpxLog, term, strikePrice, ratio, 1))/M_PI;
601 Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr,
602 cpxLog, term, strikePrice, ratio, 2))/M_PI;
608 value = spotPrice*dividendDiscount*(p1+0.5)
609 - strikePrice*riskFreeDiscount*(p2+0.5);
612 value = spotPrice*dividendDiscount*(p1-0.5)
613 - strikePrice*riskFreeDiscount*(p2-0.5);
616 QL_FAIL(
"unknown option type");
625 std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*term)/sigma;
627 const Real fwdPrice = spotPrice / ratio;
630 *M_PI/(std::sqrt(strikePrice*fwdPrice)*riskFreeDiscount);
632 const ext::function<
Real()> uM = [&](){
636 AP_Helper cvHelper(term, fwdPrice, strikePrice,
646 * std::sqrt(strikePrice * fwdPrice)/M_PI;
652 value = (cvValue + h_cv)*riskFreeDiscount;
655 value = (cvValue + h_cv - (fwdPrice - strikePrice))*riskFreeDiscount;
658 QL_FAIL(
"unknown option type");
664 QL_FAIL(
"unknown complex log formula");
672 "not an European option");
675 ext::shared_ptr<PlainVanillaPayoff> payoff =
676 ext::dynamic_pointer_cast<PlainVanillaPayoff>(
arguments_.payoff);
677 QL_REQUIRE(payoff,
"non plain vanilla payoff given");
679 const ext::shared_ptr<HestonProcess>& process =
model_->process();
681 const Real riskFreeDiscount = process->riskFreeRate()->discount(
683 const Real dividendDiscount = process->dividendYield()->discount(
686 const Real spotPrice = process->s0()->value();
687 QL_REQUIRE(spotPrice > 0.0,
"negative or null underlying given");
689 const Real strikePrice = payoff->strike();
712 ext::shared_ptr<Integrator> integrator)
713 : intAlgo_(intAlgo), integrator_(
std::move(integrator)) {}
716 Algorithm intAlgo, ext::shared_ptr<GaussianQuadrature> gaussianQuadrature)
717 : intAlgo_(intAlgo), gaussianQuadrature_(
std::move(gaussianQuadrature)) {}
722 Size maxEvaluations) {
724 ext::shared_ptr<Integrator>(
733 Size maxEvaluations) {
735 ext::shared_ptr<Integrator>(
742 Size maxEvaluations) {
744 ext::shared_ptr<Integrator>(
751 Size maxEvaluations) {
753 ext::shared_ptr<Integrator>(
760 QL_REQUIRE(intOrder <= 192,
"maximum integraton order (192) exceeded");
762 ext::shared_ptr<GaussianQuadrature>(
769 ext::shared_ptr<GaussianQuadrature>(
776 ext::shared_ptr<GaussianQuadrature>(
783 ext::shared_ptr<GaussianQuadrature>(
790 DiscreteSimpson, ext::shared_ptr<Integrator>(
797 DiscreteTrapezoid, ext::shared_ptr<Integrator>(
802 if (integrator_ !=
nullptr) {
803 return integrator_->numberOfEvaluations();
804 }
else if (gaussianQuadrature_ !=
nullptr) {
805 return gaussianQuadrature_->order();
807 QL_FAIL(
"neither Integrator nor GaussianQuadrature given");
812 return intAlgo_ == GaussLobatto
813 || intAlgo_ == GaussKronrod
814 || intAlgo_ == Simpson
815 || intAlgo_ == Trapezoid;
821 const ext::function<
Real()>& maxBound)
const {
827 retVal = (*gaussianQuadrature_)(f);
831 case GaussChebyshev2nd:
832 retVal = (*gaussianQuadrature_)(integrand1(c_inf, f));
839 retVal = (*integrator_)(f, 0.0, maxBound());
841 retVal = (*integrator_)(integrand2(c_inf, f), 0.0, 1.0);
843 case DiscreteTrapezoid:
844 case DiscreteSimpson:
846 retVal = (*integrator_)(f, 0.0, maxBound());
848 retVal = (*integrator_)(integrand3(c_inf, f), 0.0, 1.0);
851 QL_FAIL(
"unknwon integration algorithm");
860 Real maxBound)
const {
863 c_inf, f, [=](){
return maxBound; });
869 const Real uMaxGuess = -std::log(epsilon)/c_inf;
870 const Real uMaxStep = 0.1*uMaxGuess;
877 std::sqrt(uMaxGuess), 0.1*std::sqrt(uMaxGuess));
879 return std::max(uMax, uHatMax);
std::complex< Real > psi_
Real controlVariateValue() const
std::complex< Real > phi_
Real operator()(Real u) const
const ComplexLogFormula cpxLog_
static Real andersenPiterbargIntegrationLimit(Real c_inf, Real epsilon, Real v0, Real t)
Size numberOfEvaluations() const
bool isAdaptiveIntegration() const
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 gaussChebyshev(Size integrationOrder=128)
static Integration gaussLobatto(Real relTolerance, Real absTolerance, Size maxEvaluations=1000)
Real calculate(Real c_inf, const ext::function< Real(Real)> &f, const ext::function< Real()> &maxBound={}) const
static Integration trapezoid(Real absTolerance, Size maxEvaluations=1000)
static Integration gaussLaguerre(Size integrationOrder=128)
static Integration gaussKronrod(Real absTolerance, Size maxEvaluations=1000)
static Integration discreteTrapezoid(Size evaluation=1000)
analytic Heston-model engine based on Fourier transform
const Real andersenPiterbargEpsilon_
Size numberOfEvaluations() const
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)
void calculate() const override
AnalyticHestonEngine(const ext::shared_ptr< HestonModel > &model, Real relTolerance, Size maxEvaluations)
const ComplexLogFormula cpxLog_
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.
Gauss-Chebyshev integration (second kind)
Gauss-Chebyshev integration.
Integral of a 1-dimensional function using the Gauss-Kronrod methods.
generalized Gauss-Laguerre integration
Gauss-Legendre integration.
Integral of a one-dimensional function.
Base class for some pricing engine on a particular model.
Handle< ModelType > model_
Heston model for the stochastic volatility of an asset.
template class providing a null value for a given type.
Integral of a one-dimensional function.
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Integral of a one-dimensional function.
Intermediate class for put/call payoffs.
Option::Type optionType() const
Vanilla option (no discrete dividends, no barriers) on a single asset.
Real Time
continuous quantity with 1-year units
std::size_t Size
size of a container