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 analytichestonengine.hpp 22 \brief analytic Heston-model engine 23 */ 24 25 #ifndef quantlib_analytic_heston_engine_hpp 26 #define quantlib_analytic_heston_engine_hpp 27 28 #include <ql/utilities/null.hpp> 29 #include <ql/math/integrals/integral.hpp> 30 #include <ql/math/integrals/gaussianquadratures.hpp> 31 #include <ql/pricingengines/genericmodelengine.hpp> 32 #include <ql/models/equity/hestonmodel.hpp> 33 #include <ql/instruments/vanillaoption.hpp> 34 #include <ql/functional.hpp> 35 #include <complex> 36 37 namespace QuantLib { 38 39 //! analytic Heston-model engine based on Fourier transform 40 41 /*! Integration detail: 42 Two algebraically equivalent formulations of the complex 43 logarithm of the Heston model exist. Gatherals [2005] 44 (also Duffie, Pan and Singleton [2000], and Schoutens, 45 Simons and Tistaert[2004]) version does not cause 46 discoutinuities whereas the original version (e.g. Heston [1993]) 47 needs some sort of "branch correction" to work properly. 48 Gatheral's version does also work with adaptive integration 49 routines and should be preferred over the original Heston version. 50 */ 51 52 /*! References: 53 54 Heston, Steven L., 1993. A Closed-Form Solution for Options 55 with Stochastic Volatility with Applications to Bond and 56 Currency Options. The review of Financial Studies, Volume 6, 57 Issue 2, 327-343. 58 59 A. Sepp, Pricing European-Style Options under Jump Diffusion 60 Processes with Stochastic Volatility: Applications of Fourier 61 Transform (<http://math.ut.ee/~spartak/papers/stochjumpvols.pdf>) 62 63 R. Lord and C. Kahl, Why the rotation count algorithm works, 64 http://papers.ssrn.com/sol3/papers.cfm?abstract_id=921335 65 66 H. Albrecher, P. Mayer, W.Schoutens and J. Tistaert, 67 The Little Heston Trap, http://www.schoutens.be/HestonTrap.pdf 68 69 J. Gatheral, The Volatility Surface: A Practitioner's Guide, 70 Wiley Finance 71 72 F. Le Floc'h, Fourier Integration and Stochastic Volatility 73 Calibration, 74 https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2362968 75 76 L. Andersen, and V. Piterbarg, 2010, 77 Interest Rate Modeling, Volume I: Foundations and Vanilla Models, 78 Atlantic Financial Press London. 79 80 81 \ingroup vanillaengines 82 83 \test the correctness of the returned value is tested by 84 reproducing results available in web/literature 85 and comparison with Black pricing. 86 */ 87 class AnalyticHestonEngine 88 : public GenericModelEngine<HestonModel, 89 VanillaOption::arguments, 90 VanillaOption::results> { 91 public: 92 class Integration; 93 enum ComplexLogFormula { 94 // Gatheral form of characteristic function w/o control variate 95 Gatheral, 96 // old branch correction form of the characteristic function w/o control variate 97 BranchCorrection, 98 // Gatheral form with Andersen-Piterbarg control variate 99 AndersenPiterbarg, 100 // same as AndersenPiterbarg, but a slightly better control variate 101 AndersenPiterbargOptCV, 102 // Gatheral form with asymptotic expansion of the characteristic function as control variate 103 // https://hpcquantlib.wordpress.com/2020/08/30/a-novel-control-variate-for-the-heston-model 104 AsymptoticChF, 105 // auto selection of best control variate algorithm from above 106 OptimalCV 107 }; 108 109 // Simple to use constructor: Using adaptive 110 // Gauss-Lobatto integration and Gatheral's version of complex log. 111 // Be aware: using a too large number for maxEvaluations might result 112 // in a stack overflow as the Lobatto integration is a recursive 113 // algorithm. 114 AnalyticHestonEngine(const ext::shared_ptr<HestonModel>& model, 115 Real relTolerance, Size maxEvaluations); 116 117 // Constructor using Laguerre integration 118 // and Gatheral's version of complex log. 119 AnalyticHestonEngine(const ext::shared_ptr<HestonModel>& model, 120 Size integrationOrder = 144); 121 122 // Constructor giving full control 123 // over the Fourier integration algorithm 124 AnalyticHestonEngine(const ext::shared_ptr<HestonModel>& model, 125 ComplexLogFormula cpxLog, const Integration& itg, 126 Real andersenPiterbargEpsilon = 1e-8); 127 128 // normalized characteristic function 129 std::complex<Real> chF(const std::complex<Real>& z, Time t) const; 130 std::complex<Real> lnChF(const std::complex<Real>& z, Time t) const; 131 132 void calculate() const; 133 Size numberOfEvaluations() const; 134 135 static void doCalculation(Real riskFreeDiscount, 136 Real dividendDiscount, 137 Real spotPrice, 138 Real strikePrice, 139 Real term, 140 Real kappa, 141 Real theta, 142 Real sigma, 143 Real v0, 144 Real rho, 145 const TypePayoff& type, 146 const Integration& integration, 147 ComplexLogFormula cpxLog, 148 const AnalyticHestonEngine* enginePtr, 149 Real& value, 150 Size& evaluations); 151 152 static ComplexLogFormula optimalControlVariate( 153 Time t, Real v0, Real kappa, Real theta, Real sigma, Real rho); 154 155 class AP_Helper { 156 public: 157 AP_Helper(Time term, 158 Real fwd, 159 Real strike, 160 ComplexLogFormula cpxLog, 161 const AnalyticHestonEngine* enginePtr); 162 163 Real operator()(Real u) const; 164 Real controlVariateValue() const; 165 166 private: 167 const Time term_; 168 const Real fwd_, strike_, freq_; 169 const ComplexLogFormula cpxLog_; 170 const AnalyticHestonEngine* const enginePtr_; 171 Real vAvg_; 172 std::complex<Real> phi_, psi_; 173 }; 174 175 protected: 176 // call back for extended stochastic volatility 177 // plus jump diffusion engines like bates model 178 virtual std::complex<Real> addOnTerm(Real phi, 179 Time t, 180 Size j) const; 181 182 private: 183 class Fj_Helper; 184 185 mutable Size evaluations_; 186 const ComplexLogFormula cpxLog_; 187 const ext::shared_ptr<Integration> integration_; 188 const Real andersenPiterbargEpsilon_; 189 }; 190 191 192 class AnalyticHestonEngine::Integration { 193 public: 194 // non adaptive integration algorithms based on Gaussian quadrature 195 static Integration gaussLaguerre (Size integrationOrder = 128); 196 static Integration gaussLegendre (Size integrationOrder = 128); 197 static Integration gaussChebyshev (Size integrationOrder = 128); 198 static Integration gaussChebyshev2nd(Size integrationOrder = 128); 199 200 // for an adaptive integration algorithm Gatheral's version has to 201 // be used.Be aware: using a too large number for maxEvaluations might 202 // result in a stack overflow as the these integrations are based on 203 // recursive algorithms. 204 static Integration gaussLobatto(Real relTolerance, Real absTolerance, 205 Size maxEvaluations = 1000); 206 207 // usually these routines have a poor convergence behavior. 208 static Integration gaussKronrod(Real absTolerance, 209 Size maxEvaluations = 1000); 210 static Integration simpson(Real absTolerance, 211 Size maxEvaluations = 1000); 212 static Integration trapezoid(Real absTolerance, 213 Size maxEvaluations = 1000); 214 static Integration discreteSimpson(Size evaluation = 1000); 215 static Integration discreteTrapezoid(Size evaluation = 1000); 216 217 static Real andersenPiterbargIntegrationLimit( 218 Real c_inf, Real epsilon, Real v0, Real t); 219 220 Real calculate(Real c_inf, 221 const ext::function<Real(Real)>& f, 222 const ext::function<Real()>& maxBound = 223 ext::function<Real()>()) const; 224 225 Real calculate(Real c_inf, 226 const ext::function<Real(Real)>& f, Real maxBound) const; 227 228 Size numberOfEvaluations() const; 229 bool isAdaptiveIntegration() const; 230 231 private: 232 enum Algorithm 233 { GaussLobatto, GaussKronrod, Simpson, Trapezoid, 234 DiscreteTrapezoid, DiscreteSimpson, 235 GaussLaguerre, GaussLegendre, 236 GaussChebyshev, GaussChebyshev2nd }; 237 238 Integration(Algorithm intAlgo, 239 const ext::shared_ptr<GaussianQuadrature>& quadrature); 240 241 Integration(Algorithm intAlgo, 242 const ext::shared_ptr<Integrator>& integrator); 243 244 const Algorithm intAlgo_; 245 const ext::shared_ptr<Integrator> integrator_; 246 const ext::shared_ptr<GaussianQuadrature> gaussianQuadrature_; 247 }; 248 249 // inline 250 251 inline addOnTerm(Real,Time,Size) const252 std::complex<Real> AnalyticHestonEngine::addOnTerm(Real, 253 Time, 254 Size) const { 255 return std::complex<Real>(0,0); 256 } 257 } 258 259 #endif 260