41#include <boost/math/tools/minima.hpp>
42#include <boost/math/special_functions/sign.hpp>
48#if defined(QL_PATCH_MSVC)
49#pragma warning(disable: 4180)
64 return f_(-std::log(0.5-0.5*x)/c_inf_)/((1.0-x)*c_inf_);
78 return f_(-std::log(x)/c_inf_)/(x*
c_inf_);
130 class AnalyticHestonEngine::Fj_Helper {
134 const AnalyticHestonEngine* engine,
149 const Real sigma2_, rsigma_;
154 mutable Real g_km1_ = 0;
156 const AnalyticHestonEngine*
const engine_;
161 const AnalyticHestonEngine*
const engine,
162 ComplexLogFormula cpxLog,
173 sx_(
std::log(strike)),
174 dd_(
x_-
std::log(ratio)),
183 Real AnalyticHestonEngine::Fj_Helper::operator()(Real phi)
const
185 const Real rpsig(rsigma_*phi);
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 =
195 if (cpxLog_ == Gatheral) {
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));
203 std::exp(v0_*(t1-
d)*(1.0-ex)/(sigma2_*(1.0-ex*p))
205 + std::complex<Real>(0.0, phi*(dd_-sx_))
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);
216 std::exp(v0_*td*(1.0-ex)/(1.0-p*ex)
218 + std::complex<Real>(0.0, phi*(dd_-sx_))
227 if (std::fabs(kmr) > 1e-7) {
231 - v0_*(1.0-std::exp(kmr*term_)) / (2.0*kmr);
246 else if (cpxLog_ == BranchCorrection) {
247 const std::complex<Real> p = (t1+
d)/(t1-
d);
250 std::complex<Real>
g;
253 const std::complex<Real> e = std::log(p)+
d*term_;
257 g = std::log((1.0 - p/ex)/(1.0 - p));
260 g =
d*term_ + std::log(p/(p - 1.0));
267 else if (im <= -
M_PI)
270 g = std::complex<Real>(
g.real(), im);
280 const Real tmp =
g.imag() - g_km1_;
287 g += std::complex<Real>(0, 2*
b_*
M_PI);
289 return std::exp(v0_*(t1+
d)*(ex-1.0)/(sigma2_*(ex-p))
291 + std::complex<Real>(0,phi*(dd_-sx_))
296 QL_FAIL(
"unknown complex logarithm formula");
300 AnalyticHestonEngine::OptimalAlpha::OptimalAlpha(
304 fwd_(enginePtr->model_->process()->s0()->value()
305 * enginePtr->model_->process()->dividendYield()->discount(
t)
306 / enginePtr->model_->process()->riskFreeRate()->discount(
t)),
310 rho_(enginePtr->model_->
rho()),
311 eps_(
std::pow(2, -int(0.5*
std::numeric_limits<
Real>::digits))),
312 enginePtr_(enginePtr)
319 const Real eps = 1e-8;
320 const auto cm = [
this](
Real k) ->
Real {
return M(k) -
t_; };
328 cm, eps_, 0.5*(kp_ + kp_2pi), (1+eps)*kp_, (1-eps)*kp_2pi
331 else if (adx < 0.0) {
336 cm, eps_, 0.5*(kp_ + kp_pi), (1+eps)*kp_, (1-eps)*kp_pi
341 cm, eps_, 0.5*(1.0 + kp_),
349 cm, eps_, 0.5*(kp_ + kp_pi), (1+eps)*kp_, (1-eps)*kp_pi
353 QL_REQUIRE(alpha_max >= 0.0,
"alpha max must be larger than zero");
358 std::pair<Real, Real>
360 const Real alpha_max = alphaMax(strike);
362 return findMinima(eps_, std::max(2*eps_, (1.0-1e-6)*alpha_max), strike);
366 const auto cm = [
this](
Real k) ->
Real {
return M(k) -
t_; };
371 cm, eps_, 0.5*(km_2pi + km_), (1-1e-8)*km_2pi, (1+1e-8)*km_
374 QL_REQUIRE(alpha_min <= -1.0,
"alpha min must be smaller than minus one");
379 std::pair<Real, Real>
381 const Real alpha_min = alphaMin(strike);
384 std::min(-1.0-1e-6, -1.0 + (1.0-1e-6)*(alpha_min + 1.0)),
391 const std::pair<Real, Real> minusOne = alphaSmallerMinusOne(strike);
392 const std::pair<Real, Real> greaterZero = alphaGreaterZero(strike);
394 if (minusOne.second < greaterZero.second) {
395 return minusOne.first;
398 return greaterZero.first;
401 catch (
const Error&) {
413 const Real freq = std::log(fwd_/strike);
415 return boost::math::tools::brent_find_minima(
418 const std::complex<Real> z(0, -(
alpha+1));
419 return enginePtr_->lnChF(z,
t_).real()
422 lower, upper, int(0.5*std::numeric_limits<Real>::digits)
429 if (k >= km_ && k <= kp_) {
431 return std::log((
beta-D)/(
beta+D)) / D;
438 * ( ((
beta>0.0)?
M_PI : 0.0) - std::atan(D_imag/
beta) );
447 /(2*
sigma_*(1-rho_*rho_));
457 freq_(
std::log(fwd/strike)),
459 enginePtr_(enginePtr),
463 QL_REQUIRE(enginePtr !=
nullptr,
"pricing engine required");
477 vAvg_ = -8.0*std::log(enginePtr->
chF(
478 std::complex<Real>(0,
alpha_), term).real())/term;
482 * std::complex<Real>(std::sqrt(1-
rho*
rho),
rho);
484 psi_ = std::complex<Real>(
503 QL_FAIL(
"unknown control variate");
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");
514 constexpr std::complex<double> i(0, 1);
517 const std::complex<Real> h_u(u, u*tanPhi_ -
alpha_);
518 const std::complex<Real> hPrime(h_u-i);
520 std::complex<Real> phiBS(0.0);
523 -0.5*vAvg_*term_*(hPrime*hPrime +
524 std::complex<Real>(-hPrime.imag(), hPrime.real())));
526 phiBS = std::exp(u*std::complex<Real>(1, tanPhi_)*phi_ + psi_);
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)
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()))
542 return (std::exp(std::complex<Real> (0.0, u*freq_))
543 * (phiBS - enginePtr_->chF(zPrime, term_)) / (z*zPrime)
547 QL_FAIL(
"unknown control variate");
560 const std::complex<Real> phiFreq(phi_.real(), phi_.imag() + freq_);
562 using namespace ExponentialIntegral;
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();
569 return ((
alpha_ <= 0.0)? fwd_ : 0.0)
571 -0.5*((
alpha_ == 0.0)? fwd_ :0.0)
575 QL_FAIL(
"unknown control variate");
579 const std::complex<Real>& z,
Time t)
const {
581 return std::exp(
lnChF(z,
t));
593 const Real ekt = std::exp(kt);
594 const Real e2kt = std::exp(2*kt);
596 const std::complex<Real> zpi = z + std::complex<Real>(0.0, 1.0);
599 *z*zpi)/ekt)/(2.*
kappa))
601 + (std::exp(-(kt) - ((
theta -
v0 + ekt
602 *((-1 + kt)*
theta +
v0))*z*zpi)
605 *(1.0 - std::complex<Real>(-z.imag(),z.real()))*z*z)
608 + (std::exp(-2*kt - ((
theta -
v0 + ekt
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)))))
622 const std::complex<Real>& z,
Time t)
const {
632 const std::complex<Real> g
635 const std::complex<Real> D = std::sqrt(
636 g*g + (z*z + std::complex<Real>(-z.imag(), z.real()))*sigma2);
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);
644 std::complex<Real>
y;
645 if (D.real() != 0.0 || D.imag() != 0.0) {
651 const std::complex<Real> A
653 const std::complex<Real> B
654 = z*std::complex<Real>(z.real(), z.imag()+1)*
y/(1.0-
r*
y);
660 const ext::shared_ptr<HestonModel>& model,
661 Size integrationOrder)
674 const ext::shared_ptr<HestonModel>& model,
675 Real relTolerance,
Size maxEvaluations)
682 relTolerance,
Null<
Real>(), maxEvaluations))),
683 andersenPiterbargEpsilon_(1e-40),
688 const ext::shared_ptr<HestonModel>& model,
691 Real andersenPiterbargEpsilon,
699 andersenPiterbargEpsilon_(andersenPiterbargEpsilon),
703 "Branch correction does not work in conjunction "
704 "with adaptive integration methods");
726 const ext::shared_ptr<PlainVanillaPayoff>&
payoff,
727 const Date& maturity)
const {
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);
738 const ext::shared_ptr<PlainVanillaPayoff>&
payoff,
Time maturity)
const {
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);
749 const ext::shared_ptr<PlainVanillaPayoff>&
payoff,
754 const ext::shared_ptr<HestonProcess>& process =
model_->process();
755 const DiscountFactor dr = process->riskFreeRate()->discount(maturity);
758 const Real spot = process->s0()->value();
759 QL_REQUIRE(spot > 0.0,
"negative or null underlying given");
775 const Real c_inf = std::min(0.2, std::max(0.0001,
788 switch (
payoff->optionType())
791 value = spot*dd*(p1+0.5)
792 - strike*dr*(p2+0.5);
795 value = spot*dd*(p1-0.5)
796 - strike*dr*(p2-0.5);
799 QL_FAIL(
"unknown option type");
813 *
M_PI/(std::sqrt(strike*fwd)*dr);
815 const ext::function<
Real()> uM = [&](){
817 c_inf, epsilon,
v0, maturity);
825 maturity, fwd, strike, finalLog,
this,
alpha_
833 ? std::max(0.25, std::min(1000.0, 0.25/std::sqrt(0.5*vAvg*maturity)))
841 switch (
payoff->optionType())
844 value = (cvValue + h_cv)*dr;
847 value = (cvValue + h_cv - (fwd - strike))*dr;
850 QL_FAIL(
"unknown option type");
855 QL_FAIL(
"unknown complex log formula");
865 "not an European option");
868 ext::shared_ptr<PlainVanillaPayoff>
payoff =
869 ext::dynamic_pointer_cast<PlainVanillaPayoff>(
arguments_.payoff);
879 ext::shared_ptr<Integrator> integrator)
880 : intAlgo_(intAlgo), integrator_(
std::move(integrator)) {}
883 Algorithm intAlgo, ext::shared_ptr<GaussianQuadrature> gaussianQuadrature)
884 : intAlgo_(intAlgo), gaussianQuadrature_(
std::move(gaussianQuadrature)) {}
888 Real relTolerance,
Real absTolerance,
Size maxEvaluations,
bool useConvergenceEstimate) {
890 ext::shared_ptr<Integrator>(
894 useConvergenceEstimate)));
899 Size maxEvaluations) {
901 ext::shared_ptr<Integrator>(
908 Size maxEvaluations) {
910 ext::shared_ptr<Integrator>(
917 Size maxEvaluations) {
919 ext::shared_ptr<Integrator>(
926 QL_REQUIRE(intOrder <= 192,
"maximum integraton order (192) exceeded");
928 ext::shared_ptr<GaussianQuadrature>(
935 ext::shared_ptr<GaussianQuadrature>(
942 ext::shared_ptr<GaussianQuadrature>(
949 ext::shared_ptr<GaussianQuadrature>(
956 DiscreteSimpson, ext::shared_ptr<Integrator>(
963 DiscreteTrapezoid, ext::shared_ptr<Integrator>(
970 ExpSinh, ext::shared_ptr<Integrator>(
975 if (integrator_ !=
nullptr) {
976 return integrator_->numberOfEvaluations();
977 }
else if (gaussianQuadrature_ !=
nullptr) {
978 return gaussianQuadrature_->order();
980 QL_FAIL(
"neither Integrator nor GaussianQuadrature given");
985 return intAlgo_ == GaussLobatto
986 || intAlgo_ == GaussKronrod
987 || intAlgo_ == Simpson
988 || intAlgo_ == Trapezoid
989 || intAlgo_ == ExpSinh;
995 const ext::function<
Real()>& maxBound,
996 const Real scaling)
const {
1002 retVal = (*gaussianQuadrature_)(
f);
1005 case GaussChebyshev:
1006 case GaussChebyshev2nd:
1007 retVal = (*gaussianQuadrature_)(integrand1(c_inf,
f));
1010 retVal = scaling*(*integrator_)(
1011 [scaling,
f](
Real x) ->
Real {
return f(scaling*x);},
1012 0.0, std::numeric_limits<Real>::max());
1019 retVal = (*integrator_)(
f, 0.0, maxBound());
1021 retVal = (*integrator_)(integrand2(c_inf,
f), 0.0, 1.0);
1023 case DiscreteTrapezoid:
1024 case DiscreteSimpson:
1026 retVal = (*integrator_)(
f, 0.0, maxBound());
1028 retVal = (*integrator_)(integrand3(c_inf,
f), 0.0, 1.0);
1031 QL_FAIL(
"unknwon integration algorithm");
1040 Real maxBound)
const {
1043 c_inf,
f, [=](){
return maxBound; });
1049 const Real uMaxGuess = -std::log(epsilon)/c_inf;
1050 const Real uMaxStep = 0.1*uMaxGuess;
1056 const Real uHatMaxGuess = std::sqrt(-std::log(epsilon)/(0.5*
v0*
t));
1058 QL_EPSILON*uHatMaxGuess, uHatMaxGuess, 0.001*uHatMaxGuess);
1060 return std::max(uMax, uHatMax);
1062 catch (
const Error&) {
1069 Real dividendDiscount,
1079 Size& evaluations) {
1081 const Real ratio = riskFreeDiscount/dividendDiscount;
1088 const Real c_inf = std::min(0.2, std::max(0.0001,
1093 cpxLog, term, strikePrice, ratio, 1))/
M_PI;
1098 cpxLog, term, strikePrice, ratio, 2))/
M_PI;
1104 value = spotPrice*dividendDiscount*(p1+0.5)
1105 - strikePrice*riskFreeDiscount*(p2+0.5);
1108 value = spotPrice*dividendDiscount*(p1-0.5)
1109 - strikePrice*riskFreeDiscount*(p2-0.5);
1112 QL_FAIL(
"unknown option type");
1123 const Real fwdPrice = spotPrice / ratio;
1126 *
M_PI/(std::sqrt(strikePrice*fwdPrice)*riskFreeDiscount);
1128 const ext::function<
Real()> uM = [&](){
1132 AP_Helper cvHelper(term, fwdPrice, strikePrice,
1141 const Real h_cv = integration.
calculate(c_inf, cvHelper, uM)
1142 * std::sqrt(strikePrice * fwdPrice)/
M_PI;
1148 value = (cvValue + h_cv)*riskFreeDiscount;
1151 value = (cvValue + h_cv - (fwdPrice - strikePrice))*riskFreeDiscount;
1154 QL_FAIL(
"unknown option type");
1160 QL_FAIL(
"unknown complex log formula");
analytic Heston-model engine
Black-formula calculator class.
ext::shared_ptr< PricingEngine > engine_
std::complex< Real > psi_
Real controlVariateValue() const
std::complex< Real > phi_
Real operator()(Real u) const
const ComplexLogFormula cpxLog_
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)
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 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
Size numberOfEvaluations() const
std::pair< Real, Real > alphaSmallerMinusOne(Real strike) const
Real alphaMax(Real strike) const
Real k(Real x, Integer sgn) const
Real alphaMin(Real strike) const
std::pair< Real, Real > findMinima(Real lower, Real upper, Real strike) const
Real operator()(Real strike) const
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
Real priceVanillaPayoff(const ext::shared_ptr< PlainVanillaPayoff > &payoff, const Date &maturity) const
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.
ext::function< Real(Real)> f_
integrals on non uniform grids
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
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
Real Time
continuous quantity with 1-year units
Real DiscountFactor
discount factor between dates
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
ext::shared_ptr< QuantLib::Payoff > payoff
Integral of a 1-dimensional function using the Gauss-Kronrod method.
functionals and combinators not included in the STL
std::complex< Real > log1p(const std::complex< Real > &z)
std::complex< Real > expm1(const std::complex< Real > &z)
Payoffs for various options.
ext::shared_ptr< YieldTermStructure > r
integral of a one-dimensional function using Simpson formula
integral of a one-dimensional function using the trapezoid formula