analytichestonengine.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2004, 2005, 2008 Klaus Spanderen
 Copyright (C) 2007 StatPro Italia srl
 
 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/
 
 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.
 
 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/
 
#ifndef quantlib_analytic_heston_engine_hpp
#define quantlib_analytic_heston_engine_hpp
 
#include <ql/utilities/null.hpp>
#include <ql/math/integrals/integral.hpp>
#include <ql/math/integrals/gaussianquadratures.hpp>
#include <ql/pricingengines/genericmodelengine.hpp>
#include <ql/models/equity/hestonmodel.hpp>
#include <ql/instruments/vanillaoption.hpp>
#include <ql/functional.hpp>
#include <complex>
 
namespace QuantLib {
 
 
    class AnalyticHestonEngine
        : public GenericModelEngine<HestonModel,
                                    VanillaOption::arguments,
                                    VanillaOption::results> {
      public:
        class Integration;
        enum ComplexLogFormula {
            // Gatheral form of characteristic function w/o control variate
            Gatheral,
            // old branch correction form of the characteristic function w/o control variate
            BranchCorrection,
            // Gatheral form with Andersen-Piterbarg control variate
            AndersenPiterbarg,
            // same as AndersenPiterbarg, but a slightly better control variate
            AndersenPiterbargOptCV,
            // Gatheral form with asymptotic expansion of the characteristic function as control variate
            // https://hpcquantlib.wordpress.com/2020/08/30/a-novel-control-variate-for-the-heston-model
            AsymptoticChF,
            // auto selection of best control variate algorithm from above
            OptimalCV
        };
 
        // Simple to use constructor: Using adaptive
        // Gauss-Lobatto integration and Gatheral's version of complex log.
        // Be aware: using a too large number for maxEvaluations might result
        // in a stack overflow as the Lobatto integration is a recursive
        // algorithm.
        AnalyticHestonEngine(const ext::shared_ptr<HestonModel>& model,
                             Real relTolerance, Size maxEvaluations);
 
        // Constructor using Laguerre integration
        // and Gatheral's version of complex log.
        AnalyticHestonEngine(const ext::shared_ptr<HestonModel>& model,
                             Size integrationOrder = 144);
 
        // Constructor giving full control
        // over the Fourier integration algorithm
        AnalyticHestonEngine(const ext::shared_ptr<HestonModel>& model,
                             ComplexLogFormula cpxLog, const Integration& itg,
                             Real andersenPiterbargEpsilon = 1e-8);
 
        // normalized characteristic function
        std::complex<Real> chF(const std::complex<Real>& z, Time t) const;
        std::complex<Real> lnChF(const std::complex<Real>& z, Time t) const;
 
        void calculate() const override;
        Size numberOfEvaluations() 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);
 
        static ComplexLogFormula optimalControlVariate(
             Time t, Real v0, Real kappa, Real theta, Real sigma, Real rho);
 
        class AP_Helper {
          public:
            AP_Helper(Time term,
                      Real fwd,
                      Real strike,
                      ComplexLogFormula cpxLog,
                      const AnalyticHestonEngine* enginePtr);
 
            Real operator()(Real u) const;
            Real controlVariateValue() const;
 
          private:
            const Time term_;
            const Real fwd_, strike_, freq_;
            const ComplexLogFormula cpxLog_;
            const AnalyticHestonEngine* const enginePtr_;
            Real vAvg_;
            std::complex<Real> phi_, psi_;
        };
 
      protected:
        // call back for extended stochastic volatility
        // plus jump diffusion engines like bates model
        virtual std::complex<Real> addOnTerm(Real phi,
                                             Time t,
                                             Size j) const;
 
      private:
        class Fj_Helper;
 
        mutable Size evaluations_;
        const ComplexLogFormula cpxLog_;
        const ext::shared_ptr<Integration> integration_;
        const Real andersenPiterbargEpsilon_;
    };
 
 
    class AnalyticHestonEngine::Integration {
      public:
        // non adaptive integration algorithms based on Gaussian quadrature
        static Integration gaussLaguerre    (Size integrationOrder = 128);
        static Integration gaussLegendre    (Size integrationOrder = 128);
        static Integration gaussChebyshev   (Size integrationOrder = 128);
        static Integration gaussChebyshev2nd(Size integrationOrder = 128);
 
        // for an adaptive integration algorithm Gatheral's version has to
        // be used.Be aware: using a too large number for maxEvaluations might
        // result in a stack overflow as the these integrations are based on
        // recursive algorithms.
        static Integration gaussLobatto(Real relTolerance, Real absTolerance,
                                        Size maxEvaluations = 1000);
 
        // usually these routines have a poor convergence behavior.
        static Integration gaussKronrod(Real absTolerance,
                                        Size maxEvaluations = 1000);
        static Integration simpson(Real absTolerance,
                                   Size maxEvaluations = 1000);
        static Integration trapezoid(Real absTolerance,
                                     Size maxEvaluations = 1000);
        static Integration discreteSimpson(Size evaluation = 1000);
        static Integration discreteTrapezoid(Size evaluation = 1000);
 
        static Real andersenPiterbargIntegrationLimit(
            Real c_inf, Real epsilon, Real v0, Real t);
 
        Real calculate(Real c_inf,
                       const ext::function<Real(Real)>& f,
                       const ext::function<Real()>& maxBound = {}) const;
 
        Real calculate(Real c_inf,
                       const ext::function<Real(Real)>& f,
                       Real maxBound) const;
 
        Size numberOfEvaluations() const;
        bool isAdaptiveIntegration() const;
 
      private:
        enum Algorithm
            { GaussLobatto, GaussKronrod, Simpson, Trapezoid,
              DiscreteTrapezoid, DiscreteSimpson,
              GaussLaguerre, GaussLegendre,
              GaussChebyshev, GaussChebyshev2nd };
 
        Integration(Algorithm intAlgo, ext::shared_ptr<GaussianQuadrature> quadrature);
 
        Integration(Algorithm intAlgo, ext::shared_ptr<Integrator> integrator);
 
        const Algorithm intAlgo_;
        const ext::shared_ptr<Integrator> integrator_;
        const ext::shared_ptr<GaussianQuadrature> gaussianQuadrature_;
    };
 
    // inline
 
    inline 
    std::complex<Real> AnalyticHestonEngine::addOnTerm(Real,
                                                       Time,
                                                       Size) const {
        return std::complex<Real>(0,0);
    }
}
 
#endif