Skip to content

Commit 43de3fc

Browse files
lballabiolballabio
committed
Moved modified Bessel tests with tests for other functions
1 parent d55bf3f commit 43de3fc

15 files changed

+261
-260
lines changed

‎test-suite/Makefile.am

+1-1
Original file line numberDiff line numberDiff line change
@@ -40,11 +40,11 @@ QL_TESTS = \
4040
everestoption.hpp everestoption.cpp \
4141
exchangerate.hpp exchangerate.cpp \
4242
extendedtrees.hpp extendedtrees.cpp \
43-
factorial.hpp factorial.cpp \
4443
fastfouriertransform.hpp fastfouriertransform.cpp \
4544
fdheston.hpp fdheston.cpp \
4645
fdmlinearop.hpp fdmlinearop.cpp \
4746
forwardoption.hpp forwardoption.cpp \
47+
functions.hpp functions.cpp \
4848
garch.hpp garch.cpp \
4949
gaussianquadratures.hpp gaussianquadratures.cpp \
5050
gjrgarchmodel.hpp gjrgarchmodel.cpp \

‎test-suite/factorial.cpp

-117
This file was deleted.

‎test-suite/functions.cpp

+237
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,237 @@
1+
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2+
3+
/*
4+
Copyright (C) 2003 Ferdinando Ametrano
5+
Copyright (C) 2014 Klaus Spanderen
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+
#include "functions.hpp"
22+
#include "utilities.hpp"
23+
#include <ql/math/factorial.hpp>
24+
#include <ql/math/distributions/gammadistribution.hpp>
25+
#include <ql/math/modifiedbessel.hpp>
26+
27+
using namespace QuantLib;
28+
using namespace boost::unit_test_framework;
29+
30+
void FunctionsTest::testFactorial() {
31+
32+
BOOST_TEST_MESSAGE("Testing factorial numbers...");
33+
34+
Real expected = 1.0;
35+
Real calculated = Factorial::get(0);
36+
if (calculated!=expected)
37+
BOOST_FAIL("Factorial(0) = " << calculated);
38+
39+
for (Natural i=1; i<171; ++i) {
40+
expected *= i;
41+
calculated = Factorial::get(i);
42+
if (std::fabs(calculated-expected)/expected > 1.0e-9)
43+
BOOST_FAIL("Factorial(" << i << ")" <<
44+
std::setprecision(16) << QL_SCIENTIFIC <<
45+
"\n calculated: " << calculated <<
46+
"\n expected: " << expected <<
47+
"\n rel. error: " <<
48+
std::fabs(calculated-expected)/expected);
49+
}
50+
}
51+
52+
void FunctionsTest::testGammaFunction() {
53+
54+
BOOST_TEST_MESSAGE("Testing Gamma function...");
55+
56+
Real expected = 0.0;
57+
Real calculated = GammaFunction().logValue(1);
58+
if (std::fabs(calculated) > 1.0e-15)
59+
BOOST_ERROR("GammaFunction(1)\n"
60+
<< std::setprecision(16) << QL_SCIENTIFIC
61+
<< " calculated: " << calculated << "\n"
62+
<< " expected: " << expected);
63+
64+
for (Size i=2; i<9000; i++) {
65+
expected += std::log(Real(i));
66+
calculated = GammaFunction().logValue(static_cast<Real>(i+1));
67+
if (std::fabs(calculated-expected)/expected > 1.0e-9)
68+
BOOST_ERROR("GammaFunction(" << i << ")\n"
69+
<< std::setprecision(16) << QL_SCIENTIFIC
70+
<< " calculated: " << calculated << "\n"
71+
<< " expected: " << expected << "\n"
72+
<< " rel. error: "
73+
<< std::fabs(calculated-expected)/expected);
74+
}
75+
}
76+
77+
void FunctionsTest::testGammaValues() {
78+
79+
BOOST_TEST_MESSAGE("Testing Gamma values...");
80+
81+
// reference results are calculated with R
82+
Real tasks[][3] = {
83+
{ 0.0001, 9999.422883231624, 1e3},
84+
{ 1.2, 0.9181687423997607, 1e3},
85+
{ 7.3, 1271.4236336639089586, 1e3},
86+
{-1.1, 9.7148063829028946, 1e3},
87+
{-4.001,-41.6040228304425312, 1e3},
88+
{-4.999, -8.347576090315059, 1e3},
89+
{-19.000001, 8.220610833201313e-12, 1e8},
90+
{-19.5, 5.811045977502255e-18, 1e3},
91+
{-21.000001, 1.957288098276488e-14, 1e8},
92+
{-21.5, 1.318444918321553e-20, 1e6}
93+
};
94+
95+
for (Size i=0; i < LENGTH(tasks); ++i) {
96+
const Real x = tasks[i][0];
97+
const Real expected = tasks[i][1];
98+
const Real calculated = GammaFunction().value(x);
99+
const Real tol = tasks[i][2] * QL_EPSILON*std::fabs(expected);
100+
101+
if (std::fabs(calculated - expected) > tol) {
102+
BOOST_ERROR("GammaFunction(" << x << ")\n"
103+
<< std::setprecision(16) << QL_SCIENTIFIC
104+
<< " calculated: " << calculated << "\n"
105+
<< " expected: " << expected << "\n"
106+
<< " rel. error: "
107+
<< std::fabs(calculated-expected)/expected);
108+
}
109+
}
110+
}
111+
112+
void FunctionsTest::testModifiedBesselFunctions() {
113+
BOOST_TEST_MESSAGE("Testing Modified Bessel function of first kind ...");
114+
115+
/* reference values are computed with R and the additional package Bessel
116+
* http://cran.r-project.org/web/packages/Bessel
117+
*/
118+
119+
Real r[][4] = {
120+
{-1.3, 2.0, 1.2079888436539505, 0.1608243636110430},
121+
{ 1.3, 2.0, 1.2908192151358788, 0.1608243636110430},
122+
{ 0.001, 2.0, 2.2794705965773794, 0.1138938963603362},
123+
{ 1.2, 0.5, 0.1768918783499572, 2.1086579232338192},
124+
{ 2.3, 0.1, 0.00037954958988425198, 572.096866928290183},
125+
{-2.3, 1.1, 1.07222017902746969, 1.88152553684107371},
126+
{-10.0001, 1.1, 13857.7715614282552, 69288858.9474423379}
127+
};
128+
129+
for (Size i=0; i < LENGTH(r); ++i) {
130+
const Real nu = r[i][0];
131+
const Real x = r[i][1];
132+
const Real expected_i = r[i][2];
133+
const Real expected_k = r[i][3];
134+
const Real tol_i = 5e4 * QL_EPSILON*std::fabs(expected_i);
135+
const Real tol_k = 5e4 * QL_EPSILON*std::fabs(expected_k);
136+
137+
const Real calculated_i = modifiedBesselFunction_i(nu, x);
138+
const Real calculated_k = modifiedBesselFunction_k(nu, x);
139+
140+
if (std::fabs(expected_i - calculated_i) > tol_i) {
141+
BOOST_ERROR("failed to reproduce modified Bessel "
142+
<< "function of first kind"
143+
<< "\n order : " << nu
144+
<< "\n argument : " << x
145+
<< "\n calculated: " << calculated_i
146+
<< "\n expected : " << expected_i);
147+
}
148+
if (std::fabs(expected_k - calculated_k) > tol_k) {
149+
BOOST_ERROR("failed to reproduce modified Bessel "
150+
<< "function of second kind"
151+
<< "\n order : " << nu
152+
<< "\n argument : " << x
153+
<< "\n calculated: " << calculated_k
154+
<< "\n expected : " << expected_k);
155+
}
156+
}
157+
158+
Real c[][7] = {
159+
{-1.3, 2.0, 0.0, 1.2079888436539505, 0.0,
160+
0.1608243636110430, 0.0},
161+
{ 1.2, 1.5, 0.3, 0.7891550871263575, 0.2721408731632123,
162+
0.275126507673411, -0.1316314405663727},
163+
{ 1.2, -1.5,0.0,-0.6650597524355781, -0.4831941938091643,
164+
-0.251112360556051, -2.400130904230102},
165+
{-11.2, 1.5, 0.3,12780719.20252659, 16401053.26770633,
166+
-34155172.65672453, -43830147.36759921},
167+
{ 1.2, -1.5,2.0,-0.3869803778520574, 0.9756701796853728,
168+
-3.111629716783005, 0.6307859871879062},
169+
{ 1.2, 0.0, 9.9999,-0.03507838078252647, 0.1079601550451466,
170+
-0.05979939995451453, 0.3929814473878203},
171+
{ 1.2, 0.0, 10.1, -0.02782046891519293, 0.08562259917678558,
172+
-0.02035685034691133, 0.3949834389686676},
173+
{ 1.2, 0.0, 12.1, 0.07092110620741207, -0.2182727210128104,
174+
0.3368505862966958, -0.1299038064313366},
175+
{ 1.2, 0.0, 14.1,-0.03014378676768797, 0.09277303628303372,
176+
-0.237531022649052, -0.2351923034581644},
177+
{ 1.2, 0.0, 16.1,-0.03823210284792657, 0.1176663135266562,
178+
-0.1091239402448228, 0.2930535651966139},
179+
{ 1.2, 0.0, 18.1,0.05626742394733754, -0.173173324361983,
180+
0.2941636588154642, -0.02023355577954348},
181+
{ 1.2, 0.0, 180.1,-0.001230682086826484, 0.003787649998122361,
182+
0.02284509628723454, 0.09055419580980778},
183+
{ 1.2, 0.0, 21.0,-0.04746415965014021, 0.1460796627610969,
184+
-0.2693825171336859, -0.04830804448126782},
185+
{ 1.2, 10.0, 0.0, 2609.784936867044, 0, 1.904394919838336e-05, 0},
186+
{ 1.2, 14.0, 0.0, 122690.4873454286, 0, 2.902060692576643e-07, 0},
187+
{ 1.2, 20.0, 10.0, -37452017.91168936, -13917587.22151363,
188+
-3.821534367487143e-10, 4.083211255351664e-10},
189+
{ 1.2, 9.0, 9.0, -621.7335051293694, 618.1455736670332,
190+
-4.480795479964915e-05, -3.489034389148745e-08}
191+
};
192+
193+
for (Size i=0; i < LENGTH(c); ++i) {
194+
const Real nu = c[i][0];
195+
const std::complex<Real> z = std::complex<Real>(c[i][1], c[i][2]);
196+
const std::complex<Real> expected_i
197+
= std::complex<Real>(c[i][3],c[i][4]);
198+
const std::complex<Real> expected_k
199+
= std::complex<Real>(c[i][5],c[i][6]);
200+
201+
const Real tol_i = 5e4*QL_EPSILON*std::abs(expected_i);
202+
const Real tol_k = 1e6*QL_EPSILON*std::abs(expected_k);
203+
204+
const std::complex<Real> calculated_i=modifiedBesselFunction_i(nu, z);
205+
const std::complex<Real> calculated_k=modifiedBesselFunction_k(nu, z);
206+
207+
if (std::abs(expected_i - calculated_i) > tol_i) {
208+
BOOST_ERROR("failed to reproduce modified Bessel "
209+
<< "function of first kind"
210+
<< "\n order : " << nu
211+
<< "\n argument : " << z
212+
<< "\n calculated: " << calculated_i
213+
<< "\n expected : " << expected_i);
214+
}
215+
if ( std::abs(expected_k) > 1e-4 // do not check small values
216+
&& std::abs(expected_k - calculated_k) > tol_k) {
217+
BOOST_ERROR("failed to reproduce modified Bessel "
218+
<< "function of second kind"
219+
<< "\n order : " << nu
220+
<< "\n argument : " << z
221+
<< "\n diff : " << calculated_k-expected_k
222+
<< "\n calculated: " << calculated_k
223+
<< "\n expected : " << expected_k);
224+
}
225+
}
226+
}
227+
228+
test_suite* FunctionsTest::suite() {
229+
test_suite* suite = BOOST_TEST_SUITE("Factorial tests");
230+
suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testFactorial));
231+
suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testGammaFunction));
232+
suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testGammaValues));
233+
suite->add(QUANTLIB_TEST_CASE(
234+
&FunctionsTest::testModifiedBesselFunctions));
235+
return suite;
236+
}
237+

0 commit comments

Comments
 (0)