Skip to content

Commit

Permalink
Merge pull request #802.
Browse files Browse the repository at this point in the history
added Laguerre sine and cosine quadrature
  • Loading branch information
lballabio committed Apr 22, 2020
2 parents 0bfd8b5 + d00ab74 commit af45ef3
Show file tree
Hide file tree
Showing 13 changed files with 520 additions and 205 deletions.
4 changes: 3 additions & 1 deletion QuantLib.vcxproj
Original file line number Diff line number Diff line change
Expand Up @@ -1010,9 +1010,11 @@
<ClInclude Include="ql\math\integrals\filonintegral.hpp" />
<ClInclude Include="ql\math\integrals\gaussianorthogonalpolynomial.hpp" />
<ClInclude Include="ql\math\integrals\gaussianquadratures.hpp" />
<ClInclude Include="ql\math\integrals\gausslaguerrecosinepolynomial.hpp" />
<ClInclude Include="ql\math\integrals\gausslobattointegral.hpp" />
<ClInclude Include="ql\math\integrals\integral.hpp" />
<ClInclude Include="ql\math\integrals\kronrodintegral.hpp" />
<ClInclude Include="ql\math\integrals\momentbasedgaussianpolynomial.hpp" />
<ClInclude Include="ql\math\integrals\segmentintegral.hpp" />
<ClInclude Include="ql\math\integrals\simpsonintegral.hpp" />
<ClInclude Include="ql\math\integrals\trapezoidintegral.hpp" />
Expand Down Expand Up @@ -2741,4 +2743,4 @@
<Import Project=".\Build.props" Condition="Exists('.\Build.props')" />
<ImportGroup Label="ExtensionTargets">
</ImportGroup>
</Project>
</Project>
8 changes: 7 additions & 1 deletion QuantLib.vcxproj.filters
Original file line number Diff line number Diff line change
Expand Up @@ -4294,6 +4294,12 @@
<ClInclude Include="ql\methods\finitedifferences\schemes\cranknicolsonscheme.hpp">
<Filter>methods\finitedifferences\schemes</Filter>
</ClInclude>
<ClInclude Include="ql\math\integrals\gausslaguerrecosinepolynomial.hpp">
<Filter>math\integrals</Filter>
</ClInclude>
<ClInclude Include="ql\math\integrals\momentbasedgaussianpolynomial.hpp">
<Filter>math\integrals</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="ql\methods\montecarlo\brownianbridge.cpp">
Expand Down Expand Up @@ -6931,4 +6937,4 @@
<Filter>methods\finitedifferences\schemes</Filter>
</ClCompile>
</ItemGroup>
</Project>
</Project>
2 changes: 2 additions & 0 deletions ql/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -1432,9 +1432,11 @@ set(QuantLib_HDR
math/integrals/filonintegral.hpp
math/integrals/gaussianorthogonalpolynomial.hpp
math/integrals/gaussianquadratures.hpp
math/integrals/gausslaguerrecosinepolynomial.hpp
math/integrals/gausslobattointegral.hpp
math/integrals/integral.hpp
math/integrals/kronrodintegral.hpp
math/integrals/momentbasedgaussianpolynomial.hpp
math/integrals/segmentintegral.hpp
math/integrals/simpsonintegral.hpp
math/integrals/trapezoidintegral.hpp
Expand Down
184 changes: 51 additions & 133 deletions ql/experimental/math/gaussiannoncentralchisquaredpolynomial.cpp

Large diffs are not rendered by default.

55 changes: 5 additions & 50 deletions ql/experimental/math/gaussiannoncentralchisquaredpolynomial.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -24,68 +24,23 @@
#ifndef quantlib_gaussian_non_central_chi_squared_polynomial_hpp
#define quantlib_gaussian_non_central_chi_squared_polynomial_hpp

#include <ql/math/integrals/gaussianorthogonalpolynomial.hpp>

/* Define this to improve the precision of the non central chi squared
gaussian quadrature by using the boost multiprecision library.
Needs boost version > 1.52
*/
#ifndef MULTIPRECISION_NON_CENTRAL_CHI_SQUARED_QUADRATURE
//# define MULTIPRECISION_NON_CENTRAL_CHI_SQUARED_QUADRATURE
#endif

#ifdef MULTIPRECISION_NON_CENTRAL_CHI_SQUARED_QUADRATURE
#if BOOST_VERSION < 105300
#error This boost version is too old for boost multi precision support
#endif

#include <boost/multiprecision/cpp_dec_float.hpp>
#endif

#include <vector>
#include <ql/functional.hpp>
#include <ql/math/integrals/momentbasedgaussianpolynomial.hpp>

namespace QuantLib {
//! Gauss polynomial for the non central chi squared distribution
/*! References:
Gauss quadratures and orthogonal polynomials
G.H. Gloub and J.H. Welsch: Calculation of Gauss quadrature rule.
Math. Comput. 23 (1986), 221-230,
http://web.stanford.edu/class/cme335/spr11/S0025-5718-69-99647-1.pdf
M. Morandi Cecchi and M. Redivo Zaglia, Computing the coefficients
of a recurrence formula for numerical integration by moments and
modified moments.
http://ac.els-cdn.com/0377042793901522/1-s2.0-0377042793901522-main.pdf?_tid=643d5dca-a05d-11e6-9a56-00000aab0f27&acdnat=1478023545_cf7c87cba4cc9e37a136e68a2564d411
*/

class GaussNonCentralChiSquaredPolynomial
: public GaussianOrthogonalPolynomial {
: public MomentBasedGaussianPolynomial<Real> {
public:
GaussNonCentralChiSquaredPolynomial(Real nu, Real lambda);

Real mu_0() const;
Real alpha(Size i) const;
Real beta(Size i) const;
Real w(Real x) const;

#ifdef MULTIPRECISION_NON_CENTRAL_CHI_SQUARED_QUADRATURE
typedef boost::multiprecision::number<
boost::multiprecision::cpp_dec_float<100> > mp_float;
#else
typedef Real mp_float;
#endif
Real moment(Size i) const;

private:
mp_float alpha_(Size i) const;
mp_float beta_(Size i) const;

mp_float z(Integer k, Integer i) const;

const Real nu_, lambda_;

mutable std::vector<mp_float> b_, c_;
mutable std::vector<std::vector<mp_float> > z_;
static std::vector<ext::function<Real(Real, Real)> > moments;
};
}

Expand Down
2 changes: 2 additions & 0 deletions ql/math/integrals/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,13 @@ this_include_HEADERS = \
all.hpp \
discreteintegrals.hpp \
filonintegral.hpp \
gausslaguerrecosinepolynomial.hpp \
gausslobattointegral.hpp \
gaussianorthogonalpolynomial.hpp \
gaussianquadratures.hpp \
integral.hpp \
kronrodintegral.hpp \
momentbasedgaussianpolynomial.hpp \
segmentintegral.hpp \
simpsonintegral.hpp \
trapezoidintegral.hpp \
Expand Down
2 changes: 2 additions & 0 deletions ql/math/integrals/all.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,13 @@

#include <ql/math/integrals/discreteintegrals.hpp>
#include <ql/math/integrals/filonintegral.hpp>
#include <ql/math/integrals/gausslaguerrecosinepolynomial.hpp>
#include <ql/math/integrals/gausslobattointegral.hpp>
#include <ql/math/integrals/gaussianorthogonalpolynomial.hpp>
#include <ql/math/integrals/gaussianquadratures.hpp>
#include <ql/math/integrals/integral.hpp>
#include <ql/math/integrals/kronrodintegral.hpp>
#include <ql/math/integrals/momentbasedgaussianpolynomial.hpp>
#include <ql/math/integrals/segmentintegral.hpp>
#include <ql/math/integrals/simpsonintegral.hpp>
#include <ql/math/integrals/trapezoidintegral.hpp>
Expand Down
156 changes: 156 additions & 0 deletions ql/math/integrals/gausslaguerrecosinepolynomial.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,156 @@
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
Copyright (C) 2020 Klaus Spanderen
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 gausslaguerrecosinepolynomial.hpp
\brief Laguerre-Cosine Gaussian quadrature
*/

#ifndef quantlib_gauss_laguerre_cosine_polynomial_hpp
#define quantlib_gauss_laguerre_cosine_polynomial_hpp

#include <ql/math/functional.hpp>
#include <ql/math/integrals/momentbasedgaussianpolynomial.hpp>

namespace QuantLib {

template <class mp_real>
class GaussLaguerreTrigonometricBase
: public MomentBasedGaussianPolynomial<mp_real> {
public:
explicit GaussLaguerreTrigonometricBase(Real u)
: u_(u) { }

protected:
virtual mp_real m0() const = 0;
virtual mp_real m1() const = 0;

mp_real moment_(Size n) const {
if (m_.size() <= n)
m_.resize(n+1, std::numeric_limits<mp_real>::quiet_NaN());

if (boost::math::isnan(m_[n])) {
if (n == 0)
m_[0] = m0();
else if (n == 1)
m_[1] = m1();
else
m_[n] = (2*n*moment_(n-1)
- n*(n-1)*moment_(n-2))/(1+u_*u_);
}

return m_[n];
}
mp_real fact(Size n) const {
if (f_.size() <= n)
f_.resize(n+1, std::numeric_limits<mp_real>::quiet_NaN());

if (boost::math::isnan(f_[n])) {
if (n == 0)
f_[0] = 1.0;
else
f_[n] = n*fact(n-1);
}
return f_[n];

}
const Real u_;

private:
mutable std::vector<mp_real> m_, f_;
};

//! Gauss-Laguerre Cosine integration

/*! This class performs a 1-dimensional Gauss-Laguerre-Cosine integration.
\f[
\int_{0}^{\inf} f(x) \mathrm{d}x
\f]
The weighting function is
\f[
w(x;u)=e^{-x}*\cos{u*x}
\f]
*/
template <class mp_real>
class GaussLaguerreCosinePolynomial
: public GaussLaguerreTrigonometricBase<mp_real> {
public:
explicit GaussLaguerreCosinePolynomial(Real u)
: GaussLaguerreTrigonometricBase<mp_real>(u),
m0_(1.0+1.0/(1.0+u*u)) { }

mp_real moment(Size n) const {
return (this->moment_(n) + this->fact(n))/m0_;
}
Real w(Real x) const {
return std::exp(-x)*(1 + std::cos(this->u_*x))/m0_;
}

protected:
mp_real m0() const {
return 1/(1 + this->u_*this->u_);
}
mp_real m1() const {
return (1 - this->u_*this->u_)
/square<mp_real>()(1 + this->u_*this->u_);
}

private:
const Real m0_;
};

//! Gauss-Laguerre Sine integration

/*! This class performs a 1-dimensional Gauss-Laguerre-Cosine integration.
\f[
\int_{0}^{\inf} f(x) \mathrm{d}x
\f]
The weighting function is
\f[
w(x;u)=e^{-x}*\sin{u*x}
\f]
*/
template <class mp_real>
class GaussLaguerreSinePolynomial
: public GaussLaguerreTrigonometricBase<mp_real> {
public:
explicit GaussLaguerreSinePolynomial(Real u)
: GaussLaguerreTrigonometricBase<mp_real>(u),
m0_(1.0+u/(1.0+u*u)) { }

mp_real moment(Size n) const {
return (this->moment_(n) + this->fact(n))/m0_;
}
Real w(Real x) const {
return std::exp(-x)*(1 + std::sin(this->u_*x))/m0_;
}

protected:
mp_real m0() const {
return this->u_/(1 + this->u_*this->u_);
}
mp_real m1() const {
return 2*this->u_/square<mp_real>()(1 + this->u_*this->u_);
}

private:
const Real m0_;
};
}

#endif
Loading

0 comments on commit af45ef3

Please sign in to comment.