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
37namespace 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 override;
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 = {}) const;
223
224 Real calculate(Real c_inf,
225 const ext::function<Real(Real)>& f,
226 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, ext::shared_ptr<GaussianQuadrature> quadrature);
239
240 Integration(Algorithm intAlgo, ext::shared_ptr<Integrator> integrator);
241
242 const Algorithm intAlgo_;
243 const ext::shared_ptr<Integrator> integrator_;
244 const ext::shared_ptr<GaussianQuadrature> gaussianQuadrature_;
245 };
246
247 // inline
248
249 inline
250 std::complex<Real> AnalyticHestonEngine::addOnTerm(Real,
251 Time,
252 Size) const {
253 return std::complex<Real>(0,0);
254 }
255}
256
257#endif
258

source code of quantlib/ql/pricingengines/vanilla/analytichestonengine.hpp