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 | Copyright (C) 2015 Johannes Göttker-Schnetmann |
7 | |
8 | This file is part of QuantLib, a free-software/open-source library |
9 | for financial quantitative analysts and developers - http://quantlib.org/ |
10 | |
11 | QuantLib is free software: you can redistribute it and/or modify it |
12 | under the terms of the QuantLib license. You should have received a |
13 | copy of the license along with this program; if not, please email |
14 | <quantlib-dev@lists.sf.net>. The license is also available online at |
15 | <http://quantlib.org/license.shtml>. |
16 | |
17 | This program is distributed in the hope that it will be useful, but WITHOUT |
18 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
19 | FOR A PARTICULAR PURPOSE. See the license for more details. |
20 | */ |
21 | |
22 | #include "functions.hpp" |
23 | #include "utilities.hpp" |
24 | #include <ql/math/comparison.hpp> |
25 | #include <ql/math/factorial.hpp> |
26 | #include <ql/math/distributions/gammadistribution.hpp> |
27 | #include <ql/math/modifiedbessel.hpp> |
28 | |
29 | using namespace QuantLib; |
30 | using namespace boost::unit_test_framework; |
31 | |
32 | using std::exp; |
33 | |
34 | void FunctionsTest::testFactorial() { |
35 | |
36 | BOOST_TEST_MESSAGE("Testing factorial numbers..." ); |
37 | |
38 | Real expected = 1.0; |
39 | Real calculated = Factorial::get(n: 0); |
40 | if (calculated!=expected) |
41 | BOOST_FAIL("Factorial(0) = " << calculated); |
42 | |
43 | for (Natural i=1; i<171; ++i) { |
44 | expected *= i; |
45 | calculated = Factorial::get(n: i); |
46 | if (std::fabs(x: calculated-expected)/expected > 1.0e-9) |
47 | BOOST_FAIL("Factorial(" << i << ")" << |
48 | std::setprecision(16) << std::scientific << |
49 | "\n calculated: " << calculated << |
50 | "\n expected: " << expected << |
51 | "\n rel. error: " << |
52 | std::fabs(calculated-expected)/expected); |
53 | } |
54 | } |
55 | |
56 | void FunctionsTest::testGammaFunction() { |
57 | |
58 | BOOST_TEST_MESSAGE("Testing Gamma function..." ); |
59 | |
60 | Real expected = 0.0; |
61 | Real calculated = GammaFunction().logValue(x: 1); |
62 | if (std::fabs(x: calculated) > 1.0e-15) |
63 | BOOST_ERROR("GammaFunction(1)\n" |
64 | << std::setprecision(16) << std::scientific |
65 | << " calculated: " << calculated << "\n" |
66 | << " expected: " << expected); |
67 | |
68 | for (Size i=2; i<9000; i++) { |
69 | expected += std::log(x: Real(i)); |
70 | calculated = GammaFunction().logValue(x: static_cast<Real>(i+1)); |
71 | if (std::fabs(x: calculated-expected)/expected > 1.0e-9) |
72 | BOOST_ERROR("GammaFunction(" << i << ")\n" |
73 | << std::setprecision(16) << std::scientific |
74 | << " calculated: " << calculated << "\n" |
75 | << " expected: " << expected << "\n" |
76 | << " rel. error: " |
77 | << std::fabs(calculated-expected)/expected); |
78 | } |
79 | } |
80 | |
81 | void FunctionsTest::testGammaValues() { |
82 | |
83 | BOOST_TEST_MESSAGE("Testing Gamma values..." ); |
84 | |
85 | // reference results are calculated with R |
86 | Real tasks[][3] = { |
87 | { 0.0001, 9999.422883231624, 1e3}, |
88 | { 1.2, 0.9181687423997607, 1e3}, |
89 | { 7.3, 1271.4236336639089586, 1e3}, |
90 | {-1.1, 9.7148063829028946, 1e3}, |
91 | {-4.001,-41.6040228304425312, 1e3}, |
92 | {-4.999, -8.347576090315059, 1e3}, |
93 | {-19.000001, 8.220610833201313e-12, 1e8}, |
94 | {-19.5, 5.811045977502255e-18, 1e3}, |
95 | {-21.000001, 1.957288098276488e-14, 1e8}, |
96 | {-21.5, 1.318444918321553e-20, 1e6} |
97 | }; |
98 | |
99 | for (auto& task : tasks) { |
100 | const Real x = task[0]; |
101 | const Real expected = task[1]; |
102 | const Real calculated = GammaFunction().value(x); |
103 | const Real tol = task[2] * QL_EPSILON * std::fabs(x: expected); |
104 | |
105 | if (std::fabs(x: calculated - expected) > tol) { |
106 | BOOST_ERROR("GammaFunction(" << x << ")\n" |
107 | << std::setprecision(16) << std::scientific |
108 | << " calculated: " << calculated << "\n" |
109 | << " expected: " << expected << "\n" |
110 | << " rel. error: " |
111 | << std::fabs(calculated-expected)/expected); |
112 | } |
113 | } |
114 | } |
115 | |
116 | void FunctionsTest::testModifiedBesselFunctions() { |
117 | BOOST_TEST_MESSAGE("Testing modified Bessel function of first and second kind..." ); |
118 | |
119 | /* reference values are computed with R and the additional package Bessel |
120 | * http://cran.r-project.org/web/packages/Bessel |
121 | */ |
122 | |
123 | Real r[][4] = { |
124 | {-1.3, 2.0, 1.2079888436539505, 0.1608243636110430}, |
125 | { 1.3, 2.0, 1.2908192151358788, 0.1608243636110430}, |
126 | { 0.001, 2.0, 2.2794705965773794, 0.1138938963603362}, |
127 | { 1.2, 0.5, 0.1768918783499572, 2.1086579232338192}, |
128 | { 2.3, 0.1, 0.00037954958988425198, 572.096866928290183}, |
129 | {-2.3, 1.1, 1.07222017902746969, 1.88152553684107371}, |
130 | {-10.0001, 1.1, 13857.7715614282552, 69288858.9474423379} |
131 | }; |
132 | |
133 | for (auto& i : r) { |
134 | const Real nu = i[0]; |
135 | const Real x = i[1]; |
136 | const Real expected_i = i[2]; |
137 | const Real expected_k = i[3]; |
138 | const Real tol_i = 5e4 * QL_EPSILON*std::fabs(x: expected_i); |
139 | const Real tol_k = 5e4 * QL_EPSILON*std::fabs(x: expected_k); |
140 | |
141 | const Real calculated_i = modifiedBesselFunction_i(nu, x); |
142 | const Real calculated_k = modifiedBesselFunction_k(nu, x); |
143 | |
144 | if (std::fabs(x: expected_i - calculated_i) > tol_i) { |
145 | BOOST_ERROR("failed to reproduce modified Bessel " |
146 | << "function of first kind" |
147 | << "\n order : " << nu |
148 | << "\n argument : " << x |
149 | << "\n calculated: " << calculated_i |
150 | << "\n expected : " << expected_i); |
151 | } |
152 | if (std::fabs(x: expected_k - calculated_k) > tol_k) { |
153 | BOOST_ERROR("failed to reproduce modified Bessel " |
154 | << "function of second kind" |
155 | << "\n order : " << nu |
156 | << "\n argument : " << x |
157 | << "\n calculated: " << calculated_k |
158 | << "\n expected : " << expected_k); |
159 | } |
160 | } |
161 | |
162 | Real c[][7] = { |
163 | {-1.3, 2.0, 0.0, 1.2079888436539505, 0.0, |
164 | 0.1608243636110430, 0.0}, |
165 | { 1.2, 1.5, 0.3, 0.7891550871263575, 0.2721408731632123, |
166 | 0.275126507673411, -0.1316314405663727}, |
167 | { 1.2, -1.5,0.0,-0.6650597524355781, -0.4831941938091643, |
168 | -0.251112360556051, -2.400130904230102}, |
169 | {-11.2, 1.5, 0.3,12780719.20252659, 16401053.26770633, |
170 | -34155172.65672453, -43830147.36759921}, |
171 | { 1.2, -1.5,2.0,-0.3869803778520574, 0.9756701796853728, |
172 | -3.111629716783005, 0.6307859871879062}, |
173 | { 1.2, 0.0, 9.9999,-0.03507838078252647, 0.1079601550451466, |
174 | -0.05979939995451453, 0.3929814473878203}, |
175 | { 1.2, 0.0, 10.1, -0.02782046891519293, 0.08562259917678558, |
176 | -0.02035685034691133, 0.3949834389686676}, |
177 | { 1.2, 0.0, 12.1, 0.07092110620741207, -0.2182727210128104, |
178 | 0.3368505862966958, -0.1299038064313366}, |
179 | { 1.2, 0.0, 14.1,-0.03014378676768797, 0.09277303628303372, |
180 | -0.237531022649052, -0.2351923034581644}, |
181 | { 1.2, 0.0, 16.1,-0.03823210284792657, 0.1176663135266562, |
182 | -0.1091239402448228, 0.2930535651966139}, |
183 | { 1.2, 0.0, 18.1,0.05626742394733754, -0.173173324361983, |
184 | 0.2941636588154642, -0.02023355577954348}, |
185 | { 1.2, 0.0, 180.1,-0.001230682086826484, 0.003787649998122361, |
186 | 0.02284509628723454, 0.09055419580980778}, |
187 | { 1.2, 0.0, 21.0,-0.04746415965014021, 0.1460796627610969, |
188 | -0.2693825171336859, -0.04830804448126782}, |
189 | { 1.2, 10.0, 0.0, 2609.784936867044, 0, 1.904394919838336e-05, 0}, |
190 | { 1.2, 14.0, 0.0, 122690.4873454286, 0, 2.902060692576643e-07, 0}, |
191 | { 1.2, 20.0, 10.0, -37452017.91168936, -13917587.22151363, |
192 | -3.821534367487143e-10, 4.083211255351664e-10}, |
193 | { 1.2, 9.0, 9.0, -621.7335051293694, 618.1455736670332, |
194 | -4.480795479964915e-05, -3.489034389148745e-08} |
195 | }; |
196 | |
197 | for (auto& i : c) { |
198 | const Real nu = i[0]; |
199 | const std::complex<Real> z = std::complex<Real>(i[1], i[2]); |
200 | const std::complex<Real> expected_i = std::complex<Real>(i[3], i[4]); |
201 | const std::complex<Real> expected_k = std::complex<Real>(i[5], i[6]); |
202 | |
203 | const Real tol_i = 5e4*QL_EPSILON*std::abs(z: expected_i); |
204 | const Real tol_k = 1e6*QL_EPSILON*std::abs(z: expected_k); |
205 | |
206 | const std::complex<Real> calculated_i=modifiedBesselFunction_i(nu, z); |
207 | const std::complex<Real> calculated_k=modifiedBesselFunction_k(nu, z); |
208 | |
209 | if (std::abs(z: expected_i - calculated_i) > tol_i) { |
210 | BOOST_ERROR("failed to reproduce modified Bessel " |
211 | << "function of first kind" |
212 | << "\n order : " << nu |
213 | << "\n argument : " << z |
214 | << "\n calculated: " << calculated_i |
215 | << "\n expected : " << expected_i); |
216 | } |
217 | if ( std::abs(z: expected_k) > 1e-4 // do not check small values |
218 | && std::abs(z: expected_k - calculated_k) > tol_k) { |
219 | BOOST_ERROR("failed to reproduce modified Bessel " |
220 | << "function of second kind" |
221 | << "\n order : " << nu |
222 | << "\n argument : " << z |
223 | << "\n diff : " << calculated_k-expected_k |
224 | << "\n calculated: " << calculated_k |
225 | << "\n expected : " << expected_k); |
226 | } |
227 | } |
228 | } |
229 | |
230 | void FunctionsTest::testWeightedModifiedBesselFunctions() { |
231 | BOOST_TEST_MESSAGE("Testing weighted modified Bessel functions..." ); |
232 | for (Real nu = -5.0; nu <= 5.0; nu += 0.5) { |
233 | for (Real x = 0.1; x <= 15.0; x += 0.5) { |
234 | Real calculated_i = modifiedBesselFunction_i_exponentiallyWeighted(nu, x); |
235 | Real expected_i = modifiedBesselFunction_i(nu, x) * exp(x: -x); |
236 | Real calculated_k = modifiedBesselFunction_k_exponentiallyWeighted(nu, x); |
237 | Real expected_k = |
238 | M_PI_2 * (modifiedBesselFunction_i(nu: -nu, x) - modifiedBesselFunction_i(nu, x)) * |
239 | exp(x: -x) / std::sin(M_PI * nu); |
240 | Real tol_i = 1e3 * QL_EPSILON * std::fabs(x: expected_i) * std::max(a: exp(x: x), b: 1.0); |
241 | Real tol_k = std::max(QL_EPSILON, b: 1e3 * QL_EPSILON * std::fabs(x: expected_k) * |
242 | std::max(a: exp(x: x), b: 1.0)); |
243 | if (std::abs(x: expected_i - calculated_i) > tol_i) { |
244 | BOOST_ERROR("failed to verify exponentially weighted" |
245 | << "modified Bessel function of first kind" |
246 | << "\n order : " << nu << "\n argument : " << x |
247 | << "\n calculated : " << calculated_i << "\n expected : " |
248 | << expected_i << "\n difference : " << (expected_i - calculated_i)); |
249 | } |
250 | if (std::abs(x: expected_k - calculated_k) > tol_k) { |
251 | BOOST_ERROR("failed to verify exponentially weighted" |
252 | << "modified Bessel function of second kind" |
253 | << "\n order : " << nu << "\n argument : " << x |
254 | << "\n calculated : " << calculated_k << "\n expected : " |
255 | << expected_k << "\n difference : " << (expected_k - calculated_k)); |
256 | } |
257 | } |
258 | } |
259 | for (Real nu = -5.0; nu <= 5.0; nu += 0.5) { |
260 | for (Real x = -5.0; x <= 5.0; x += 0.5) { |
261 | for (Real y = -5.0; y <= 5.0; y += 0.5) { |
262 | std::complex<Real> z(x, y); |
263 | std::complex<Real> calculated_i = |
264 | modifiedBesselFunction_i_exponentiallyWeighted(nu, z); |
265 | std::complex<Real> expected_i = modifiedBesselFunction_i(nu, z) * exp(z: -z); |
266 | std::complex<Real> calculated_k = |
267 | modifiedBesselFunction_k_exponentiallyWeighted(nu, z); |
268 | std::complex<Real> expected_k = M_PI_2 * |
269 | (modifiedBesselFunction_i(nu: -nu, z) * exp(z: -z) - |
270 | modifiedBesselFunction_i(nu, z) * exp(z: -z)) / |
271 | std::sin(M_PI * nu); |
272 | Real tol_i = 1e3 * QL_EPSILON * std::abs(z: calculated_i); |
273 | Real tol_k = 1e3 * QL_EPSILON * std::abs(z: calculated_k); |
274 | if (std::abs(z: calculated_i - expected_i) > tol_i) { |
275 | BOOST_ERROR("failed to verify exponentially weighted" |
276 | << "modified Bessel function of first kind" |
277 | << "\n order : " << nu << "\n argument : " << x |
278 | << "\n calculated : " << calculated_i << "\n expected : " |
279 | << expected_i << "\n difference : " << (expected_i - calculated_i)); |
280 | } |
281 | if (std::abs(z: expected_k - calculated_k) > tol_k) { |
282 | BOOST_ERROR("failed to verify exponentially weighted" |
283 | << "modified Bessel function of second kind" |
284 | << "\n order : " << nu << "\n argument : " << x |
285 | << "\n calculated : " << calculated_k << "\n expected : " |
286 | << expected_k << "\n difference : " << (expected_k - calculated_k)); |
287 | } |
288 | } |
289 | } |
290 | } |
291 | } |
292 | |
293 | test_suite* FunctionsTest::suite() { |
294 | auto* suite = BOOST_TEST_SUITE("Factorial tests" ); |
295 | suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testFactorial)); |
296 | suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testGammaFunction)); |
297 | suite->add(QUANTLIB_TEST_CASE(&FunctionsTest::testGammaValues)); |
298 | suite->add(QUANTLIB_TEST_CASE( |
299 | &FunctionsTest::testModifiedBesselFunctions)); |
300 | suite->add(QUANTLIB_TEST_CASE( |
301 | &FunctionsTest::testWeightedModifiedBesselFunctions)); |
302 | return suite; |
303 | } |
304 | |