/src/quantlib/ql/math/integrals/fourierintegration.cpp

Source
/* -*- 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/tanhsinhintegral.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) {
        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));
    }

    FourierIntegration FourierIntegration::tanhSinh(Real relTolerance) {
        return FourierIntegration(
            TanhSinh, ext::make_shared<TanhSinhIntegral>(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
            || intAlgo_ == TanhSinh;
    }

    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 TanhSinh:
            retVal = (*integrator_)(f, 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;
        }
    }
}

Coverage Report

Created: 2026-10-31 06:47

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