30110271200

Committed 24 Jul 2026 04:44PM UTC coverage: 74.943% (+0.06%) from 74.885%

Build # 30110271200

Build Type

push

github

Committed by

web-flow

Commit Message

Rough Heston Model with Semi-Analytic Fourier Pricing Engine and Fractional ODE Solver (#2648)

Coverage Stats

283 of 299 new or added lines in 5 files covered. (94.65%)

1 existing line in 1 file now uncovered.

59962 of 80010 relevant lines covered (74.94%)

8659768.11 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

87.88

/ql/math/integrals/fourierintegration.cpp

/* -*- 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
 <https://www.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.
*/

#include <ql/math/integrals/discreteintegrals.hpp>
#include <ql/math/integrals/expsinhintegral.hpp>
#include <ql/math/integrals/fourierintegration.hpp>
#include <ql/math/integrals/gausslobattointegral.hpp>
#include <ql/math/integrals/kronrodintegral.hpp>
#include <ql/math/integrals/simpsonintegral.hpp>
#include <ql/math/integrals/trapezoidintegral.hpp>
#include <ql/math/solvers1d/brent.hpp>
#include <ql/utilities/null.hpp>
#include <cmath>
#include <limits>
#include <utility>

namespace QuantLib {

    namespace {

        class integrand1 {
          private:
            const Real c_inf_;
            const std::function<Real(Real)> f_;
          public:
            integrand1(Real c_inf, std::function<Real(Real)> f) : c_inf_(c_inf), f_(std::move(f)) {}
            Real operator()(Real x) const {
                if ((1.0-x)*c_inf_ > QL_EPSILON)
                    return f_(-std::log(0.5-0.5*x)/c_inf_)/((1.0-x)*c_inf_);
                else
                    return 0.0;
            }
        };

        class integrand2 {
          private:
            const Real c_inf_;
            const std::function<Real(Real)> f_;
          public:
            integrand2(Real c_inf, std::function<Real(Real)> f) : c_inf_(c_inf), f_(std::move(f)) {}
            Real operator()(Real x) const {
                if (x*c_inf_ > QL_EPSILON) {
                    return f_(-std::log(x)/c_inf_)/(x*c_inf_);
                } else {
                    return 0.0;
                }
            }
        };

        class integrand3 {
          private:
            const integrand2 int_;
          public:
            integrand3(Real c_inf, const std::function<Real(Real)>& f)
            : int_(c_inf, f) {}

            Real operator()(Real x) const { return int_(1.0-x); }
        };

        class u_Max {
          public:
            u_Max(Real c_inf, Real epsilon) : c_inf_(c_inf), logEpsilon_(std::log(epsilon)) {}

            Real operator()(Real u) const {
                ++evaluations_;
                return c_inf_*u + std::log(u) + logEpsilon_;
            }

            Size evaluations() const { return evaluations_; }

          private:
            const Real c_inf_, logEpsilon_;
            mutable Size evaluations_ = 0;
        };


        class uHat_Max {
          public:
            uHat_Max(Real v0T2, Real epsilon) : v0T2_(v0T2), logEpsilon_(std::log(epsilon)) {}

            Real operator()(Real u) const {
                ++evaluations_;
                return v0T2_*u*u + std::log(u) + logEpsilon_;
            }

            Size evaluations() const { return evaluations_; }

          private:
            const Real v0T2_, logEpsilon_;
            mutable Size evaluations_ = 0;
        };
    }


    FourierIntegration::FourierIntegration(Algorithm intAlgo,
                                           ext::shared_ptr<Integrator> integrator)
    : intAlgo_(intAlgo), integrator_(std::move(integrator)) {}

    FourierIntegration::FourierIntegration(
        Algorithm intAlgo, ext::shared_ptr<GaussianQuadrature> gaussianQuadrature)
    : intAlgo_(intAlgo), gaussianQuadrature_(std::move(gaussianQuadrature)) {}

    FourierIntegration FourierIntegration::gaussLobatto(
       Real relTolerance, Real absTolerance, Size maxEvaluations, bool useConvergenceEstimate) {
       return FourierIntegration(GaussLobatto,
                           ext::make_shared<GaussLobattoIntegral>(maxEvaluations,
                                                        absTolerance,
                                                        relTolerance,
                                                        useConvergenceEstimate));
    }

    FourierIntegration FourierIntegration::gaussKronrod(Real absTolerance,
                                                        Size maxEvaluations) {
        return FourierIntegration(GaussKronrod,
                           ext::make_shared<GaussKronrodAdaptive>(absTolerance,
                                                        maxEvaluations));
    }

    FourierIntegration FourierIntegration::simpson(Real absTolerance,
                                                   Size maxEvaluations) {
        return FourierIntegration(Simpson,
                           ext::make_shared<SimpsonIntegral>(absTolerance,
                                                   maxEvaluations));
    }

    FourierIntegration FourierIntegration::trapezoid(Real absTolerance,
                                                     Size maxEvaluations) {
        return FourierIntegration(Trapezoid,
                           ext::make_shared<TrapezoidIntegral<Default>>(absTolerance,
                                                             maxEvaluations));
    }

    FourierIntegration FourierIntegration::gaussLaguerre(Size intOrder) {
        QL_REQUIRE(intOrder <= 192, "maximum integraton order (192) exceeded");
        return FourierIntegration(GaussLaguerre,
                           ext::make_shared<GaussLaguerreIntegration>(intOrder));
    }

    FourierIntegration FourierIntegration::gaussLegendre(Size intOrder) {
        return FourierIntegration(GaussLegendre,
                           ext::make_shared<GaussLegendreIntegration>(intOrder));
    }

    FourierIntegration FourierIntegration::gaussChebyshev(Size intOrder) {
        return FourierIntegration(GaussChebyshev,
                           ext::make_shared<GaussChebyshevIntegration>(intOrder));
    }

    FourierIntegration FourierIntegration::gaussChebyshev2nd(Size intOrder) {
        return FourierIntegration(GaussChebyshev2nd,
                           ext::make_shared<GaussChebyshev2ndIntegration>(intOrder));
    }

    FourierIntegration FourierIntegration::discreteSimpson(Size evaluations) {
        return FourierIntegration(
            DiscreteSimpson, ext::make_shared<DiscreteSimpsonIntegrator>(evaluations));
    }

    FourierIntegration FourierIntegration::discreteTrapezoid(Size evaluations) {
        return FourierIntegration(
            DiscreteTrapezoid, ext::make_shared<DiscreteTrapezoidIntegrator>(evaluations));
    }

    FourierIntegration FourierIntegration::expSinh(Real relTolerance) {
        return FourierIntegration(
            ExpSinh, ext::make_shared<ExpSinhIntegral>(relTolerance));
    }

    Size FourierIntegration::numberOfEvaluations() const {
        if (integrator_ != nullptr) {
            return integrator_->numberOfEvaluations();
        } else if (gaussianQuadrature_ != nullptr) {
            return gaussianQuadrature_->order();
        } else {
            QL_FAIL("neither Integrator nor GaussianQuadrature given");
        }
    }

    bool FourierIntegration::isAdaptiveIntegration() const {
        return intAlgo_ == GaussLobatto
            || intAlgo_ == GaussKronrod
            || intAlgo_ == Simpson
            || intAlgo_ == Trapezoid
            || intAlgo_ == ExpSinh;
    }

    Real FourierIntegration::calculate(
        Real c_inf,
        const std::function<Real(Real)>& f,
        const std::function<Real()>& maxBound,
        const Real scaling) const {

        Real retVal;

        switch(intAlgo_) {
          case GaussLaguerre:
            retVal = (*gaussianQuadrature_)(f);
            break;
          case GaussLegendre:
          case GaussChebyshev:
          case GaussChebyshev2nd:
            retVal = (*gaussianQuadrature_)(integrand1(c_inf, f));
            break;
          case ExpSinh:
            retVal = scaling*(*integrator_)(
                [scaling, f](Real x) -> Real { return f(scaling*x);},
                0.0, std::numeric_limits<Real>::max());
            break;
          case Simpson:
          case Trapezoid:
          case GaussLobatto:
          case GaussKronrod:
              if (maxBound && maxBound() != Null<Real>())
                  retVal = (*integrator_)(f, 0.0, maxBound());
              else
                  retVal = (*integrator_)(integrand2(c_inf, f), 0.0, 1.0);
              break;
          case DiscreteTrapezoid:
          case DiscreteSimpson:
              if (maxBound && maxBound() != Null<Real>())
                  retVal = (*integrator_)(f, 0.0, maxBound());
              else
                  retVal = (*integrator_)(integrand3(c_inf, f), 0.0, 1.0);
              break;
          default:
            QL_FAIL("unknwon integration algorithm");
        }

        return retVal;
     }

    Real FourierIntegration::calculate(
        Real c_inf,
        const std::function<Real(Real)>& f,
        Real maxBound) const {

        return FourierIntegration::calculate(
            c_inf, f, [=](){ return maxBound; });
    }

    Real FourierIntegration::andersenPiterbargIntegrationLimit(
        Real c_inf, Real epsilon, Real v0, Real t) {

        const Real uMaxGuess = -std::log(epsilon)/c_inf;
        const Real uMaxStep = 0.1*uMaxGuess;

        const Real uMax = Brent().solve(u_Max(c_inf, epsilon),
            QL_EPSILON*uMaxGuess, uMaxGuess, uMaxStep);

        try {
            const Real uHatMaxGuess = std::sqrt(-std::log(epsilon)/(0.5*v0*t));
            const Real uHatMax = Brent().solve(uHat_Max(0.5*v0*t, epsilon),
                QL_EPSILON*uHatMaxGuess, uHatMaxGuess, 0.001*uHatMaxGuess);

            return std::max(uMax, uHatMax);
        }
        catch (const Error&) {
            return uMax;
        }
    }
}

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	<https://www.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	#include <ql/math/integrals/discreteintegrals.hpp>
22	#include <ql/math/integrals/expsinhintegral.hpp>
23	#include <ql/math/integrals/fourierintegration.hpp>
24	#include <ql/math/integrals/gausslobattointegral.hpp>
25	#include <ql/math/integrals/kronrodintegral.hpp>
26	#include <ql/math/integrals/simpsonintegral.hpp>
27	#include <ql/math/integrals/trapezoidintegral.hpp>
28	#include <ql/math/solvers1d/brent.hpp>
29	#include <ql/utilities/null.hpp>
30	#include <cmath>
31	#include <limits>
32	#include <utility>
33
34	namespace QuantLib {
35
36	namespace {
37
NEW 38	class integrand1 {	×
39	private:
40	const Real c_inf_;
41	const std::function<Real(Real)> f_;
42	public:
43	integrand1(Real c_inf, std::function<Real(Real)> f) : c_inf_(c_inf), f_(std::move(f)) {}	2,535✔
44	Real operator()(Real x) const {	1,296,128✔
45	if ((1.0-x)*c_inf_ > QL_EPSILON)	1,296,128✔
46	return f_(-std::log(0.5-0.5x)/c_inf_)/((1.0-x)c_inf_);	1,296,128✔
47	else
48	return 0.0;
49	}
50	};
51
52	class integrand2 {	24✔
53	private:
54	const Real c_inf_;
55	const std::function<Real(Real)> f_;
56	public:
57	integrand2(Real c_inf, std::function<Real(Real)> f) : c_inf_(c_inf), f_(std::move(f)) {}	14✔
58	Real operator()(Real x) const {	94,784✔
59	if (x*c_inf_ > QL_EPSILON) {	94,784✔
60	return f_(-std::log(x)/c_inf_)/(x*c_inf_);	94,768✔
61	} else {
62	return 0.0;
63	}
64	}
65	};
66
67	class integrand3 {	8✔
68	private:
69	const integrand2 int_;
70	public:
71	integrand3(Real c_inf, const std::function<Real(Real)>& f)	4✔
72	: int_(c_inf, f) {}	4✔
73
74	Real operator()(Real x) const { return int_(1.0-x); }	2,048✔
75	};
76
77	class u_Max {
78	public:
79	u_Max(Real c_inf, Real epsilon) : c_inf_(c_inf), logEpsilon_(std::log(epsilon)) {}	117,845✔
80
81	Real operator()(Real u) const {	1,005,050✔
82	++evaluations_;	1,005,050✔
83	return c_inf_*u + std::log(u) + logEpsilon_;	1,005,050✔
84	}
85
86	Size evaluations() const { return evaluations_; }
87
88	private:
89	const Real c_inf_, logEpsilon_;
90	mutable Size evaluations_ = 0;
91	};
92
93
94	class uHat_Max {
95	public:
96	uHat_Max(Real v0T2, Real epsilon) : v0T2_(v0T2), logEpsilon_(std::log(epsilon)) {}	117,845✔
97
98	Real operator()(Real u) const {	1,511,427✔
99	++evaluations_;	1,511,427✔
100	return v0T2_uu + std::log(u) + logEpsilon_;	1,511,427✔
101	}
102
103	Size evaluations() const { return evaluations_; }
104
105	private:
106	const Real v0T2_, logEpsilon_;
107	mutable Size evaluations_ = 0;
108	};
109	}
110
111
112	FourierIntegration::FourierIntegration(Algorithm intAlgo,	42✔
113	ext::shared_ptr<Integrator> integrator)	42✔
114	: intAlgo_(intAlgo), integrator_(std::move(integrator)) {}	42✔
115
116	FourierIntegration::FourierIntegration(	2,016✔
117	Algorithm intAlgo, ext::shared_ptr<GaussianQuadrature> gaussianQuadrature)	2,016✔
118	: intAlgo_(intAlgo), gaussianQuadrature_(std::move(gaussianQuadrature)) {}	2,016✔
119
120	FourierIntegration FourierIntegration::gaussLobatto(	19✔
121	Real relTolerance, Real absTolerance, Size maxEvaluations, bool useConvergenceEstimate) {
122	return FourierIntegration(GaussLobatto,
123	ext::make_shared<GaussLobattoIntegral>(maxEvaluations,	19✔
124	absTolerance,
125	relTolerance,
126	useConvergenceEstimate));	38✔
127	}
128
129	FourierIntegration FourierIntegration::gaussKronrod(Real absTolerance,	2✔
130	Size maxEvaluations) {
131	return FourierIntegration(GaussKronrod,
132	ext::make_shared<GaussKronrodAdaptive>(absTolerance,	2✔
133	maxEvaluations));	4✔
134	}
135
136	FourierIntegration FourierIntegration::simpson(Real absTolerance,	2✔
137	Size maxEvaluations) {
138	return FourierIntegration(Simpson,
139	ext::make_shared<SimpsonIntegral>(absTolerance,	2✔
140	maxEvaluations));	4✔
141	}
142
143	FourierIntegration FourierIntegration::trapezoid(Real absTolerance,	3✔
144	Size maxEvaluations) {
145	return FourierIntegration(Trapezoid,
146	ext::make_shared<TrapezoidIntegral<Default>>(absTolerance,	3✔
147	maxEvaluations));	6✔
148	}
149
150	FourierIntegration FourierIntegration::gaussLaguerre(Size intOrder) {	1,992✔
151	QL_REQUIRE(intOrder <= 192, "maximum integraton order (192) exceeded");	1,992✔
152	return FourierIntegration(GaussLaguerre,
153	ext::make_shared<GaussLaguerreIntegration>(intOrder));	3,984✔
154	}
155
156	FourierIntegration FourierIntegration::gaussLegendre(Size intOrder) {	8✔
157	return FourierIntegration(GaussLegendre,
158	ext::make_shared<GaussLegendreIntegration>(intOrder));	16✔
159	}
160
161	FourierIntegration FourierIntegration::gaussChebyshev(Size intOrder) {	8✔
162	return FourierIntegration(GaussChebyshev,
163	ext::make_shared<GaussChebyshevIntegration>(intOrder));	16✔
164	}
165
166	FourierIntegration FourierIntegration::gaussChebyshev2nd(Size intOrder) {	8✔
167	return FourierIntegration(GaussChebyshev2nd,
168	ext::make_shared<GaussChebyshev2ndIntegration>(intOrder));	16✔
169	}
170
171	FourierIntegration FourierIntegration::discreteSimpson(Size evaluations) {	3✔
172	return FourierIntegration(
173	DiscreteSimpson, ext::make_shared<DiscreteSimpsonIntegrator>(evaluations));	6✔
174	}
175
176	FourierIntegration FourierIntegration::discreteTrapezoid(Size evaluations) {	13✔
177	return FourierIntegration(
178	DiscreteTrapezoid, ext::make_shared<DiscreteTrapezoidIntegrator>(evaluations));	26✔
179	}
180
NEW 181	FourierIntegration FourierIntegration::expSinh(Real relTolerance) {	×
182	return FourierIntegration(
NEW 183	ExpSinh, ext::make_shared<ExpSinhIntegral>(relTolerance));	×
184	}
185
186	Size FourierIntegration::numberOfEvaluations() const {	307,910✔
187	if (integrator_ != nullptr) {	307,910✔
188	return integrator_->numberOfEvaluations();	58,936✔
189	} else if (gaussianQuadrature_ != nullptr) {	248,974✔
190	return gaussianQuadrature_->order();	248,974✔
191	} else {
NEW 192	QL_FAIL("neither Integrator nor GaussianQuadrature given");	×
193	}
194	}
195
196	bool FourierIntegration::isAdaptiveIntegration() const {	30✔
197	return intAlgo_ == GaussLobatto	30✔
198	\|\| intAlgo_ == GaussKronrod
199	\|\| intAlgo_ == Simpson
200	\|\| intAlgo_ == Trapezoid
201	\|\| intAlgo_ == ExpSinh;	30✔
202	}
203
204	Real FourierIntegration::calculate(	307,912✔
205	Real c_inf,
206	const std::function<Real(Real)>& f,
207	const std::function<Real()>& maxBound,
208	const Real scaling) const {
209
210	Real retVal;
211
212	switch(intAlgo_) {	307,912✔
213	case GaussLaguerre:	246,439✔
214	retVal = (*gaussianQuadrature_)(f);	246,439✔
215	break;	246,439✔
216	case GaussLegendre:	2,535✔
217	case GaussChebyshev:
218	case GaussChebyshev2nd:
219	retVal = (*gaussianQuadrature_)(integrand1(c_inf, f));	7,605✔
220	break;	2,535✔
NEW 221	case ExpSinh:	×
NEW 222	retVal = scaling(integrator_)(	×
NEW 223	[scaling, f](Real x) -> Real { return f(scaling*x);},	×
224	0.0, std::numeric_limits<Real>::max());
NEW 225	break;	×
226	case Simpson:
227	case Trapezoid:
228	case GaussLobatto:
229	case GaussKronrod:
230	if (maxBound && maxBound() != Null<Real>())	620✔
231	retVal = (*integrator_)(f, 0.0, maxBound());	610✔
232	else
233	retVal = (*integrator_)(integrand2(c_inf, f), 0.0, 1.0);	30✔
234	break;
235	case DiscreteTrapezoid:
236	case DiscreteSimpson:
237	if (maxBound && maxBound() != Null<Real>())	58,318✔
238	retVal = (*integrator_)(f, 0.0, maxBound());	58,314✔
239	else
240	retVal = (*integrator_)(integrand3(c_inf, f), 0.0, 1.0);	8✔
241	break;
NEW 242	default:	×
NEW 243	QL_FAIL("unknwon integration algorithm");	×
244	}
245
246	return retVal;	307,912✔
247	}
248
249	Real FourierIntegration::calculate(	1✔
250	Real c_inf,
251	const std::function<Real(Real)>& f,
252	Real maxBound) const {
253
254	return FourierIntegration::calculate(	1✔
255	c_inf, f, [=](){ return maxBound; });	1✔
256	}
257
258	Real FourierIntegration::andersenPiterbargIntegrationLimit(	117,845✔
259	Real c_inf, Real epsilon, Real v0, Real t) {
260
261	const Real uMaxGuess = -std::log(epsilon)/c_inf;	117,845✔
262	const Real uMaxStep = 0.1*uMaxGuess;	117,845✔
263
264	const Real uMax = Brent().solve(u_Max(c_inf, epsilon),	117,845✔
265	QL_EPSILON*uMaxGuess, uMaxGuess, uMaxStep);	117,845✔
266
267	try {
268	const Real uHatMaxGuess = std::sqrt(-std::log(epsilon)/(0.5v0t));	117,845✔
269	const Real uHatMax = Brent().solve(uHat_Max(0.5v0t, epsilon),	117,845✔
270	QL_EPSILONuHatMaxGuess, uHatMaxGuess, 0.001uHatMaxGuess);	117,845✔
271
272	return std::max(uMax, uHatMax);	117,845✔
273	}
NEW 274	catch (const Error&) {	×
275	return uMax;
NEW 276	}	×
277	}
278	}

lballabio / QuantLib / 30110271200

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous