text: /* -*- 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.
  */
  
  /*! \file analytichestonengine.hpp
      \brief analytic Heston-model engine
  */
  
  #ifndef quantlib_analytic_heston_engine_hpp
  #define quantlib_analytic_heston_engine_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 <boost/function.hpp>
  #include <complex>
  
  namespace QuantLib {
  
      //! analytic Heston-model engine based on Fourier transform
  
      /*! Integration detail:
          Two algebraically equivalent formulations of the complex
          logarithm of the Heston model exist. Gatherals [2005]
          (also Duffie, Pan and Singleton [2000], and Schoutens,
          Simons and Tistaert[2004]) version does not cause
          discoutinuities whereas the original version (e.g. Heston [1993])
          needs some sort of "branch correction" to work properly.
          Gatheral's version does also work with adaptive integration
          routines and should be preferred over the original Heston version.
      */
  
      /*! References:
  
          Heston, Steven L., 1993. A Closed-Form Solution for Options
          with Stochastic Volatility with Applications to Bond and
          Currency Options.  The review of Financial Studies, Volume 6,
          Issue 2, 327-343.
  
          A. Sepp, Pricing European-Style Options under Jump Diffusion
          Processes with Stochastic Volatility: Applications of Fourier
          Transform (<http://math.ut.ee/~spartak/papers/stochjumpvols.pdf>)
  
          R. Lord and C. Kahl, Why the rotation count algorithm works,
          http://papers.ssrn.com/sol3/papers.cfm?abstract_id=921335
  
          H. Albrecher, P. Mayer, W.Schoutens and J. Tistaert,
          The Little Heston Trap, http://www.schoutens.be/HestonTrap.pdf
  
          J. Gatheral, The Volatility Surface: A Practitioner's Guide,
          Wiley Finance
  
          \ingroup vanillaengines
  
          \test the correctness of the returned value is tested by
                reproducing results available in web/literature
                and comparison with Black pricing.
      */
      class AnalyticHestonEngine
          : public GenericModelEngine<HestonModel,
                                      VanillaOption::arguments,
                                      VanillaOption::results> {
        public:
          class Integration;
          enum ComplexLogFormula { Gatheral, BranchCorrection };
  
          // 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 boost::shared_ptr<HestonModel>& model,
                               Real relTolerance, Size maxEvaluations);
  
          // Constructor using Laguerre integration
          // and Gatheral's version of complex log.
          AnalyticHestonEngine(const boost::shared_ptr<HestonModel>& model,
                               Size integrationOrder = 144);
  
          // Constructor giving full control
          // over the Fourier integration algorithm
          AnalyticHestonEngine(const boost::shared_ptr<HestonModel>& model,
                               ComplexLogFormula cpxLog, const Integration& itg);
  
  
          void calculate() const;
          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,
                                               const ComplexLogFormula cpxLog,
                                               const AnalyticHestonEngine* const enginePtr,
                                               Real& value,
                                               Size& evaluations);
  
        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 boost::shared_ptr<Integration> integration_;
      };
  
  
      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);
  
          Real calculate(Real c_inf,
                         const boost::function1<Real, Real>& f) const;
  
          Size numberOfEvaluations() const;
          bool isAdaptiveIntegration() const;
  
        private:
          enum Algorithm
              { GaussLobatto, GaussKronrod, Simpson, Trapezoid,
                GaussLaguerre, GaussLegendre,
                GaussChebyshev, GaussChebyshev2nd };
  
          Integration(Algorithm intAlgo,
                      const boost::shared_ptr<GaussianQuadrature>& quadrature);
  
          Integration(Algorithm intAlgo,
                      const boost::shared_ptr<Integrator>& integrator);
  
          const Algorithm intAlgo_;
          const boost::shared_ptr<Integrator> integrator_;
          const boost::shared_ptr<GaussianQuadrature> gaussianQuadrature_;
      };
  
      // inline
  
      inline 
      std::complex<Real> AnalyticHestonEngine::addOnTerm(Real,
                                                         Time,
                                                         Size) const {
          return std::complex<Real>(0,0);
      }
  }
  
  #endif
  