1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003, 2004 Ferdinando Ametrano
5 Copyright (C) 2003 StatPro Italia srl
6 Copyright (C) 2005 Gary Kennedy
7 Copyright (C) 2013 Fabien Le Floc'h
8 Copyright (C) 2016 Klaus Spanderen
9
10
11 This file is part of QuantLib, a free-software/open-source library
12 for financial quantitative analysts and developers - http://quantlib.org/
13
14 QuantLib is free software: you can redistribute it and/or modify it
15 under the terms of the QuantLib license. You should have received a
16 copy of the license along with this program; if not, please email
17 <quantlib-dev@lists.sf.net>. The license is also available online at
18 <http://quantlib.org/license.shtml>.
19
20 This program is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
22 FOR A PARTICULAR PURPOSE. See the license for more details.
23*/
24
25#include "distributions.hpp"
26#include "utilities.hpp"
27#include <ql/math/distributions/normaldistribution.hpp>
28#include <ql/math/distributions/bivariatenormaldistribution.hpp>
29#include <ql/math/distributions/bivariatestudenttdistribution.hpp>
30#include <ql/math/distributions/chisquaredistribution.hpp>
31#include <ql/math/distributions/poissondistribution.hpp>
32#include <ql/math/randomnumbers/stochasticcollocationinvcdf.hpp>
33#include <ql/math/comparison.hpp>
34#include <boost/math/distributions/non_central_chi_squared.hpp>
35
36using namespace QuantLib;
37using namespace boost::unit_test_framework;
38
39namespace distributions_test {
40
41 Real average = 1.0, sigma = 2.0;
42
43 Real gaussian(Real x) {
44 Real normFact = sigma*std::sqrt(x: 2*M_PI);
45 Real dx = x-average;
46 return std::exp(x: -dx*dx/(2.0*sigma*sigma))/normFact;
47 }
48
49 Real gaussianDerivative(Real x) {
50 Real normFact = sigma*sigma*sigma*std::sqrt(x: 2*M_PI);
51 Real dx = x-average;
52 return -dx*std::exp(x: -dx*dx/(2.0*sigma*sigma))/normFact;
53 }
54
55 struct BivariateTestData {
56 Real a;
57 Real b;
58 Real rho;
59 Real result;
60 };
61
62 template <class Bivariate>
63 void checkBivariate(const char* tag) {
64
65 BivariateTestData values[] = {
66 /* The data below are from
67 "Option pricing formulas", E.G. Haug, McGraw-Hill 1998
68 pag 193
69 */
70 { .a: 0.0, .b: 0.0, .rho: 0.0, .result: 0.250000 },
71 { .a: 0.0, .b: 0.0, .rho: -0.5, .result: 0.166667 },
72 { .a: 0.0, .b: 0.0, .rho: 0.5, .result: 1.0/3 },
73 { .a: 0.0, .b: -0.5, .rho: 0.0, .result: 0.154269 },
74 { .a: 0.0, .b: -0.5, .rho: -0.5, .result: 0.081660 },
75 { .a: 0.0, .b: -0.5, .rho: 0.5, .result: 0.226878 },
76 { .a: 0.0, .b: 0.5, .rho: 0.0, .result: 0.345731 },
77 { .a: 0.0, .b: 0.5, .rho: -0.5, .result: 0.273122 },
78 { .a: 0.0, .b: 0.5, .rho: 0.5, .result: 0.418340 },
79
80 { .a: -0.5, .b: 0.0, .rho: 0.0, .result: 0.154269 },
81 { .a: -0.5, .b: 0.0, .rho: -0.5, .result: 0.081660 },
82 { .a: -0.5, .b: 0.0, .rho: 0.5, .result: 0.226878 },
83 { .a: -0.5, .b: -0.5, .rho: 0.0, .result: 0.095195 },
84 { .a: -0.5, .b: -0.5, .rho: -0.5, .result: 0.036298 },
85 { .a: -0.5, .b: -0.5, .rho: 0.5, .result: 0.163319 },
86 { .a: -0.5, .b: 0.5, .rho: 0.0, .result: 0.213342 },
87 { .a: -0.5, .b: 0.5, .rho: -0.5, .result: 0.145218 },
88 { .a: -0.5, .b: 0.5, .rho: 0.5, .result: 0.272239 },
89
90 { .a: 0.5, .b: 0.0, .rho: 0.0, .result: 0.345731 },
91 { .a: 0.5, .b: 0.0, .rho: -0.5, .result: 0.273122 },
92 { .a: 0.5, .b: 0.0, .rho: 0.5, .result: 0.418340 },
93 { .a: 0.5, .b: -0.5, .rho: 0.0, .result: 0.213342 },
94 { .a: 0.5, .b: -0.5, .rho: -0.5, .result: 0.145218 },
95 { .a: 0.5, .b: -0.5, .rho: 0.5, .result: 0.272239 },
96 { .a: 0.5, .b: 0.5, .rho: 0.0, .result: 0.478120 },
97 { .a: 0.5, .b: 0.5, .rho: -0.5, .result: 0.419223 },
98 { .a: 0.5, .b: 0.5, .rho: 0.5, .result: 0.546244 },
99
100 // known analytical values
101 { .a: 0.0, .b: 0.0, .rho: std::sqrt(x: 1/2.0), .result: 3.0/8},
102
103 // { 0.0, big, any, 0.500000 },
104 { .a: 0.0, .b: 30, .rho: -1.0, .result: 0.500000 },
105 { .a: 0.0, .b: 30, .rho: 0.0, .result: 0.500000 },
106 { .a: 0.0, .b: 30, .rho: 1.0, .result: 0.500000 },
107
108 // { big, big, any, 1.000000 },
109 { .a: 30, .b: 30, .rho: -1.0, .result: 1.000000 },
110 { .a: 30, .b: 30, .rho: 0.0, .result: 1.000000 },
111 { .a: 30, .b: 30, .rho: 1.0, .result: 1.000000 },
112
113 // {-big, any, any, 0.000000 }
114 { .a: -30, .b: -1.0, .rho: -1.0, .result: 0.000000 },
115 { .a: -30, .b: 0.0, .rho: -1.0, .result: 0.000000 },
116 { .a: -30, .b: 1.0, .rho: -1.0, .result: 0.000000 },
117 { .a: -30, .b: -1.0, .rho: 0.0, .result: 0.000000 },
118 { .a: -30, .b: 0.0, .rho: 0.0, .result: 0.000000 },
119 { .a: -30, .b: 1.0, .rho: 0.0, .result: 0.000000 },
120 { .a: -30, .b: -1.0, .rho: 1.0, .result: 0.000000 },
121 { .a: -30, .b: 0.0, .rho: 1.0, .result: 0.000000 },
122 { .a: -30, .b: 1.0, .rho: 1.0, .result: 0.000000 }
123 };
124
125 for (Size i=0; i<LENGTH(values); i++) {
126 Bivariate bcd(values[i].rho);
127 Real value = bcd(values[i].a, values[i].b);
128
129 Real tolerance = 1.0e-6;
130 if (std::fabs(x: value-values[i].result) >= tolerance) {
131 BOOST_ERROR(tag << " bivariate cumulative distribution\n"
132 << " case: " << i+1 << "\n"
133 << std::fixed
134 << " a: " << values[i].a << "\n"
135 << " b: " << values[i].b << "\n"
136 << " rho: " << values[i].rho <<"\n"
137 << std::scientific
138 << " tabulated value: "
139 << values[i].result << "\n"
140 << " result: " << value);
141 }
142 }
143 }
144
145 template <class Bivariate>
146 void checkBivariateAtZero(const char* tag, Real tolerance) {
147
148 /*
149 BVN(0.0,0.0,rho) = 1/4 + arcsin(rho)/(2*M_PI)
150 "Handbook of the Normal Distribution",
151 J.K. Patel & C.B.Read, 2nd Ed, 1996
152 */
153 const Real rho[] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5,
154 0.6, 0.7, 0.8, 0.9, 0.99999 };
155 const Real x(0.0);
156 const Real y(0.0);
157
158 for (Real i : rho) {
159 for (Integer sgn=-1; sgn < 2; sgn+=2) {
160 Bivariate bvn(sgn * i);
161 Real expected = 0.25 + std::asin(x: sgn * i) / (2 * M_PI);
162 Real realised = bvn(x,y);
163
164 if (std::fabs(x: realised-expected)>=tolerance) {
165 BOOST_ERROR(tag << " bivariate cumulative distribution\n"
166 << std::scientific << " rho: " << sgn * i << "\n"
167 << " expected: " << expected << "\n"
168 << " realised: " << realised << "\n"
169 << " tolerance: " << tolerance);
170 }
171 }
172 }
173 }
174
175 template <class Bivariate>
176 void checkBivariateTail(const char* tag, Real tolerance) {
177
178 /* make sure numerical greeks are sensible, numerical error in
179 * the tails can make garbage greeks for partial time barrier
180 * option */
181 Real x = -6.9;
182 Real y = 6.9;
183 Real corr = -0.999;
184 Bivariate bvn(corr);
185 for (int i = 0; i<10;i++) {
186 Real cdf0 = bvn(x,y);
187 y = y + tolerance;
188 Real cdf1 = bvn(x,y);
189 if (cdf0 > cdf1) {
190 BOOST_ERROR(tag << " cdf must be decreasing in the tails\n"
191 << std::scientific
192 << " cdf0: " << cdf0 << "\n"
193 << " cdf1: " << cdf1 << "\n"
194 << " x: " << x << "\n"
195 << " y: " << y << "\n"
196 << " rho: " << corr);
197 }
198 }
199 }
200
201 struct BivariateStudentTestData {
202 Natural n;
203 Real rho;
204 Real x;
205 Real y;
206 Real result;
207 };
208
209}
210
211void DistributionTest::testNormal() {
212
213 BOOST_TEST_MESSAGE("Testing normal distributions...");
214
215 using namespace distributions_test;
216
217 InverseCumulativeNormal invCumStandardNormal;
218 Real check = invCumStandardNormal(0.5);
219 if (check != 0.0e0) {
220 BOOST_ERROR("C++ inverse cumulative of the standard normal at 0.5 is "
221 << std::scientific << check
222 << "\n instead of zero: something is wrong!");
223 }
224
225 NormalDistribution normal(average,sigma);
226 CumulativeNormalDistribution cum(average,sigma);
227 InverseCumulativeNormal invCum(average,sigma);
228
229 Size numberOfStandardDeviation = 6;
230 Real xMin = average - numberOfStandardDeviation*sigma,
231 xMax = average + numberOfStandardDeviation*sigma;
232 Size N = 100001;
233 Real h = (xMax-xMin)/(N-1);
234
235 std::vector<Real> x(N), y(N), yd(N), temp(N), diff(N);
236
237 Size i;
238 for (i=0; i<N; i++)
239 x[i] = xMin+h*i;
240 std::transform(first: x.begin(), last: x.end(), result: y.begin(), unary_op: gaussian);
241 std::transform(first: x.begin(), last: x.end(), result: yd.begin(), unary_op: gaussianDerivative);
242
243 // check that normal = Gaussian
244 std::transform(first: x.begin(), last: x.end(), result: temp.begin(), unary_op: normal);
245 std::transform(first1: y.begin(), last1: y.end(), first2: temp.begin(), result: diff.begin(), binary_op: std::minus<>());
246 Real e = norm(begin: diff.begin(), end: diff.end(), h);
247 if (e > 1.0e-16) {
248 BOOST_ERROR("norm of C++ NormalDistribution minus analytic Gaussian: "
249 << std::scientific << e << "\n"
250 << "tolerance exceeded");
251 }
252
253 // check that invCum . cum = identity
254 std::transform(first: x.begin(), last: x.end(), result: temp.begin(), unary_op: cum);
255 std::transform(first: temp.begin(), last: temp.end(), result: temp.begin(), unary_op: invCum);
256 std::transform(first1: x.begin(), last1: x.end(), first2: temp.begin(), result: diff.begin(), binary_op: std::minus<>());
257 e = norm(begin: diff.begin(), end: diff.end(), h);
258 if (e > 1.0e-7) {
259 BOOST_ERROR("norm of invCum . cum minus identity: "
260 << std::scientific << e << "\n"
261 << "tolerance exceeded");
262 }
263
264 MaddockInverseCumulativeNormal mInvCum(average, sigma);
265 std::transform(first: x.begin(), last: x.end(), result: diff.begin(),
266 unary_op: [&](Real x) -> Real {
267 return x - mInvCum(cum(x));
268 });
269
270 e = norm(begin: diff.begin(), end: diff.end(), h);
271 if (e > 1.0e-7) {
272 BOOST_ERROR("norm of MaddokInvCum . cum minus identity: "
273 << std::scientific << e << "\n"
274 << "tolerance exceeded");
275 }
276
277 // check that cum.derivative = Gaussian
278 for (i=0; i<x.size(); i++)
279 temp[i] = cum.derivative(x: x[i]);
280 std::transform(first1: y.begin(), last1: y.end(), first2: temp.begin(), result: diff.begin(), binary_op: std::minus<>());
281 e = norm(begin: diff.begin(), end: diff.end(), h);
282 if (e > 1.0e-16) {
283 BOOST_ERROR(
284 "norm of C++ Cumulative.derivative minus analytic Gaussian: "
285 << std::scientific << e << "\n"
286 << "tolerance exceeded");
287 }
288
289 // check that normal.derivative = gaussianDerivative
290 for (i=0; i<x.size(); i++)
291 temp[i] = normal.derivative(x: x[i]);
292 std::transform(first1: yd.begin(), last1: yd.end(), first2: temp.begin(), result: diff.begin(), binary_op: std::minus<>());
293 e = norm(begin: diff.begin(), end: diff.end(), h);
294 if (e > 1.0e-16) {
295 BOOST_ERROR("norm of C++ Normal.derivative minus analytic derivative: "
296 << std::scientific << e << "\n"
297 << "tolerance exceeded");
298 }
299}
300
301void DistributionTest::testBivariate() {
302
303 BOOST_TEST_MESSAGE("Testing bivariate cumulative normal distribution...");
304
305 using namespace distributions_test;
306
307 checkBivariateAtZero<BivariateCumulativeNormalDistributionDr78>(
308 tag: "Drezner 1978", tolerance: 1.0e-6);
309 checkBivariate<BivariateCumulativeNormalDistributionDr78>(tag: "Drezner 1978");
310
311 // due to relative low accuracy of Dr78, it does not pass with a
312 // smaller perturbation
313 checkBivariateTail<BivariateCumulativeNormalDistributionDr78>(
314 tag: "Drezner 1978", tolerance: 1.0e-5);
315
316 checkBivariateAtZero<BivariateCumulativeNormalDistributionWe04DP>(
317 tag: "West 2004", tolerance: 1.0e-15);
318 checkBivariate<BivariateCumulativeNormalDistributionWe04DP>(tag: "West 2004");
319
320 checkBivariateTail<BivariateCumulativeNormalDistributionWe04DP>(
321 tag: "West 2004", tolerance: 1.0e-6);
322 checkBivariateTail<BivariateCumulativeNormalDistributionWe04DP>(
323 tag: "West 2004", tolerance: 1.0e-8);
324}
325
326
327void DistributionTest::testPoisson() {
328
329 BOOST_TEST_MESSAGE("Testing Poisson distribution...");
330
331 for (Real mean=0.0; mean<=10.0; mean+=0.5) {
332 BigNatural i = 0;
333 PoissonDistribution pdf(mean);
334 Real calculated = pdf(i);
335 Real logHelper = -mean;
336 Real expected = std::exp(x: logHelper);
337 Real error = std::fabs(x: calculated-expected);
338 if (error > 1.0e-16)
339 BOOST_ERROR("Poisson pdf(" << mean << ")(" << i << ")\n"
340 << std::setprecision(16)
341 << " calculated: " << calculated << "\n"
342 << " expected: " << expected << "\n"
343 << " error: " << error);
344
345 for (i=1; i<25; i++) {
346 calculated = pdf(i);
347 if (mean == 0.0) {
348 expected = 0.0;
349 } else {
350 logHelper = logHelper+std::log(x: mean)-std::log(x: Real(i));
351 expected = std::exp(x: logHelper);
352 }
353 error = std::fabs(x: calculated-expected);
354 if (error>1.0e-13)
355 BOOST_ERROR("Poisson pdf(" << mean << ")(" << i << ")\n"
356 << std::setprecision(13)
357 << " calculated: " << calculated << "\n"
358 << " expected: " << expected << "\n"
359 << " error: " << error);
360 }
361 }
362}
363
364void DistributionTest::testCumulativePoisson() {
365
366 BOOST_TEST_MESSAGE("Testing cumulative Poisson distribution...");
367
368 for (Real mean=0.0; mean<=10.0; mean+=0.5) {
369 BigNatural i = 0;
370 CumulativePoissonDistribution cdf(mean);
371 Real cumCalculated = cdf(i);
372 Real logHelper = -mean;
373 Real cumExpected = std::exp(x: logHelper);
374 Real error = std::fabs(x: cumCalculated-cumExpected);
375 if (error>1.0e-13)
376 BOOST_ERROR("Poisson cdf(" << mean << ")(" << i << ")\n"
377 << std::setprecision(13)
378 << " calculated: " << cumCalculated << "\n"
379 << " expected: " << cumExpected << "\n"
380 << " error: " << error);
381 for (i=1; i<25; i++) {
382 cumCalculated = cdf(i);
383 if (mean == 0.0) {
384 cumExpected = 1.0;
385 } else {
386 logHelper = logHelper+std::log(x: mean)-std::log(x: Real(i));
387 cumExpected += std::exp(x: logHelper);
388 }
389 error = std::fabs(x: cumCalculated-cumExpected);
390 if (error>1.0e-12)
391 BOOST_ERROR("Poisson cdf(" << mean << ")(" << i << ")\n"
392 << std::setprecision(12)
393 << " calculated: " << cumCalculated << "\n"
394 << " expected: " << cumExpected << "\n"
395 << " error: " << error);
396 }
397 }
398}
399
400void DistributionTest::testInverseCumulativePoisson() {
401
402 BOOST_TEST_MESSAGE("Testing inverse cumulative Poisson distribution...");
403
404 InverseCumulativePoisson icp(1.0);
405
406 Real data[] = { 0.2,
407 0.5,
408 0.9,
409 0.98,
410 0.99,
411 0.999,
412 0.9999,
413 0.99995,
414 0.99999,
415 0.999999,
416 0.9999999,
417 0.99999999
418 };
419
420 for (Size i=0; i<LENGTH(data); i++) {
421 if (!close(x: icp(data[i]), y: static_cast<Real>(i))) {
422 BOOST_ERROR(std::setprecision(8)
423 << "failed to reproduce known value for x = "
424 << data[i] << "\n"
425 << " calculated: " << icp(data[i]) << "\n"
426 << " expected: " << Real(i));
427 }
428 }
429}
430
431
432void DistributionTest::testBivariateCumulativeStudent() {
433 BOOST_TEST_MESSAGE(
434 "Testing bivariate cumulative Student t distribution...");
435
436 using namespace distributions_test;
437
438 Real xs[14] = { 0.00, 0.50, 1.00, 1.50, 2.00, 2.50, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00 };
439 Natural ns[20] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 60, 90, 120, 150, 300, 600 };
440 // Part of table 1 from the reference paper
441 Real expected1[280] = {
442 0.33333, 0.50000, 0.63497, 0.72338, 0.78063, 0.81943, 0.84704, 0.88332, 0.90590, 0.92124, 0.93231, 0.94066, 0.94719, 0.95243,
443 0.33333, 0.52017, 0.68114, 0.78925, 0.85607, 0.89754, 0.92417, 0.95433, 0.96978, 0.97862, 0.98411, 0.98774, 0.99026, 0.99208,
444 0.33333, 0.52818, 0.70018, 0.81702, 0.88720, 0.92812, 0.95238, 0.97667, 0.98712, 0.99222, 0.99497, 0.99657, 0.99756, 0.99821,
445 0.33333, 0.53245, 0.71052, 0.83231, 0.90402, 0.94394, 0.96612, 0.98616, 0.99353, 0.99664, 0.99810, 0.99885, 0.99927, 0.99951,
446 0.33333, 0.53510, 0.71701, 0.84196, 0.91449, 0.95344, 0.97397, 0.99095, 0.99637, 0.99836, 0.99918, 0.99956, 0.99975, 0.99985,
447 0.33333, 0.53689, 0.72146, 0.84862, 0.92163, 0.95972, 0.97893, 0.99365, 0.99779, 0.99913, 0.99962, 0.99982, 0.99990, 0.99995,
448 0.33333, 0.53819, 0.72470, 0.85348, 0.92679, 0.96415, 0.98230, 0.99531, 0.99857, 0.99950, 0.99981, 0.99992, 0.99996, 0.99998,
449 0.33333, 0.53917, 0.72716, 0.85719, 0.93070, 0.96743, 0.98470, 0.99639, 0.99903, 0.99970, 0.99990, 0.99996, 0.99998, 0.99999,
450 0.33333, 0.53994, 0.72909, 0.86011, 0.93375, 0.96995, 0.98650, 0.99713, 0.99931, 0.99981, 0.99994, 0.99998, 0.99999, 1.00000,
451 0.33333, 0.54056, 0.73065, 0.86247, 0.93621, 0.97194, 0.98788, 0.99766, 0.99950, 0.99988, 0.99996, 0.99999, 1.00000, 1.00000,
452 0.33333, 0.54243, 0.73540, 0.86968, 0.94362, 0.97774, 0.99168, 0.99890, 0.99985, 0.99998, 1.00000, 1.00000, 1.00000, 1.00000,
453 0.33333, 0.54338, 0.73781, 0.87336, 0.94735, 0.98053, 0.99337, 0.99932, 0.99993, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000,
454 0.33333, 0.54395, 0.73927, 0.87560, 0.94959, 0.98216, 0.99430, 0.99952, 0.99996, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
455 0.33333, 0.54433, 0.74025, 0.87709, 0.95108, 0.98322, 0.99489, 0.99963, 0.99998, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
456 0.33333, 0.54528, 0.74271, 0.88087, 0.95482, 0.98580, 0.99623, 0.99983, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
457 0.33333, 0.54560, 0.74354, 0.88215, 0.95607, 0.98663, 0.99664, 0.99987, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
458 0.33333, 0.54576, 0.74396, 0.88279, 0.95669, 0.98704, 0.99683, 0.99989, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
459 0.33333, 0.54586, 0.74420, 0.88317, 0.95706, 0.98729, 0.99695, 0.99990, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
460 0.33333, 0.54605, 0.74470, 0.88394, 0.95781, 0.98777, 0.99717, 0.99992, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
461 0.33333, 0.54615, 0.74495, 0.88432, 0.95818, 0.98801, 0.99728, 0.99993, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000
462 };
463 // Part of table 2 from the reference paper
464 Real expected2[280] = {
465 0.16667, 0.36554, 0.54022, 0.65333, 0.72582, 0.77465, 0.80928, 0.85466, 0.88284, 0.90196, 0.91575, 0.92616, 0.93429, 0.94081,
466 0.16667, 0.38889, 0.59968, 0.73892, 0.82320, 0.87479, 0.90763, 0.94458, 0.96339, 0.97412, 0.98078, 0.98518, 0.98823, 0.99044,
467 0.16667, 0.39817, 0.62478, 0.77566, 0.86365, 0.91391, 0.94330, 0.97241, 0.98483, 0.99086, 0.99410, 0.99598, 0.99714, 0.99790,
468 0.16667, 0.40313, 0.63863, 0.79605, 0.88547, 0.93396, 0.96043, 0.98400, 0.99256, 0.99614, 0.99782, 0.99868, 0.99916, 0.99944,
469 0.16667, 0.40620, 0.64740, 0.80900, 0.89902, 0.94588, 0.97007, 0.98972, 0.99591, 0.99816, 0.99909, 0.99951, 0.99972, 0.99983,
470 0.16667, 0.40829, 0.65345, 0.81794, 0.90820, 0.95368, 0.97607, 0.99290, 0.99755, 0.99904, 0.99958, 0.99980, 0.99989, 0.99994,
471 0.16667, 0.40980, 0.65788, 0.82449, 0.91482, 0.95914, 0.98010, 0.99482, 0.99844, 0.99946, 0.99979, 0.99991, 0.99996, 0.99998,
472 0.16667, 0.41095, 0.66126, 0.82948, 0.91981, 0.96314, 0.98295, 0.99605, 0.99895, 0.99968, 0.99989, 0.99996, 0.99998, 0.99999,
473 0.16667, 0.41185, 0.66393, 0.83342, 0.92369, 0.96619, 0.98506, 0.99689, 0.99926, 0.99980, 0.99994, 0.99998, 0.99999, 1.00000,
474 0.16667, 0.41257, 0.66608, 0.83661, 0.92681, 0.96859, 0.98667, 0.99748, 0.99946, 0.99987, 0.99996, 0.99999, 1.00000, 1.00000,
475 0.16667, 0.41476, 0.67268, 0.84633, 0.93614, 0.97550, 0.99103, 0.99884, 0.99984, 0.99998, 1.00000, 1.00000, 1.00000, 1.00000,
476 0.16667, 0.41586, 0.67605, 0.85129, 0.94078, 0.97877, 0.99292, 0.99930, 0.99993, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000,
477 0.16667, 0.41653, 0.67810, 0.85430, 0.94356, 0.98066, 0.99396, 0.99950, 0.99996, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
478 0.16667, 0.41698, 0.67947, 0.85632, 0.94540, 0.98189, 0.99461, 0.99962, 0.99998, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
479 0.16667, 0.41810, 0.68294, 0.86141, 0.94998, 0.98483, 0.99607, 0.99982, 0.99999, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
480 0.16667, 0.41847, 0.68411, 0.86312, 0.95149, 0.98577, 0.99651, 0.99987, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
481 0.16667, 0.41866, 0.68470, 0.86398, 0.95225, 0.98623, 0.99672, 0.99989, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
482 0.16667, 0.41877, 0.68505, 0.86449, 0.95270, 0.98650, 0.99684, 0.99990, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
483 0.16667, 0.41900, 0.68576, 0.86552, 0.95360, 0.98705, 0.99707, 0.99992, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,
484 0.16667, 0.41911, 0.68612, 0.86604, 0.95405, 0.98731, 0.99719, 0.99993, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000
485 };
486
487 Real tolerance = 1.0e-5;
488 for (Size i=0; i < LENGTH(ns); ++i) {
489 BivariateCumulativeStudentDistribution f1(ns[i], 0.5);
490 BivariateCumulativeStudentDistribution f2(ns[i], -0.5);
491 for (Size j=0; j < LENGTH(xs); ++j) {
492 Real calculated1 = f1(xs[j], xs[j]);
493 Real reference1 = expected1[i*LENGTH(xs)+j];
494 Real calculated2 = f2(xs[j], xs[j]);
495 Real reference2 = expected2[i*LENGTH(xs)+j];
496 if (std::fabs(x: calculated1 - reference1) > tolerance)
497 BOOST_ERROR("Failed to reproduce CDF value at " << xs[j] <<
498 "\n calculated: " << calculated1 <<
499 "\n expected: " << reference1);
500 if (std::fabs(x: calculated2 - reference2) > tolerance)
501 BOOST_ERROR("Failed to reproduce CDF value at " << xs[j] <<
502 "\n calculated: " << calculated2 <<
503 "\n expected: " << reference1);
504 }
505 }
506
507 // a few more random cases
508 BivariateStudentTestData cases[] = {
509 {.n: 2, .rho: -1.0, .x: 5.0, .y: 8.0, .result: 0.973491},
510 {.n: 2, .rho: 1.0, .x: -2.0, .y: 8.0, .result: 0.091752},
511 {.n: 2, .rho: 1.0, .x: 5.25, .y: -9.5, .result: 0.005450},
512 {.n: 3, .rho: -0.5, .x: -5.0, .y: -5.0, .result: 0.000220},
513 {.n: 4, .rho: -1.0, .x: -8.0, .y: 7.5, .result: 0.0},
514 {.n: 4, .rho: 0.5, .x: -5.5, .y: 10.0, .result: 0.002655},
515 {.n: 4, .rho: 1.0, .x: -5.0, .y: 6.0, .result: 0.003745},
516 {.n: 4, .rho: 1.0, .x: 6.0, .y: 5.5, .result: 0.997336},
517 {.n: 5, .rho: -0.5, .x: -7.0, .y: -6.25, .result: 0.000004},
518 {.n: 5, .rho: -0.5, .x: 3.75, .y: -7.25, .result: 0.000166},
519 {.n: 5, .rho: -0.5, .x: 7.75, .y: -1.25, .result: 0.133073},
520 {.n: 6, .rho: 0.0, .x: 7.5, .y: 3.25, .result: 0.991149},
521 {.n: 7, .rho: -0.5, .x: -1.0, .y: -8.5, .result: 0.000001},
522 {.n: 7, .rho: -1.0, .x: -4.25, .y: -4.0, .result: 0.0},
523 {.n: 7, .rho: 0.0, .x: 0.5, .y: -2.25, .result: 0.018819},
524 {.n: 8, .rho: -1.0, .x: 8.25, .y: 1.75, .result: 0.940866},
525 {.n: 8, .rho: 0.0, .x: 2.25, .y: 4.75, .result: 0.972105},
526 {.n: 9, .rho: -0.5, .x: -4.0, .y: 8.25, .result: 0.001550},
527 {.n: 9, .rho: -1.0, .x: -1.25, .y: -8.75, .result: 0.0},
528 {.n: 9, .rho: -1.0, .x: 5.75, .y: -6.0, .result: 0.0},
529 {.n: 9, .rho: 0.5, .x: -6.5, .y: -9.5, .result: 0.000001},
530 {.n: 9, .rho: 1.0, .x: -2.0, .y: 9.25, .result: 0.038276},
531 {.n: 10, .rho: -1.0, .x: -0.5, .y: 6.0, .result: 0.313881},
532 {.n: 10, .rho: 0.5, .x: 0.0, .y: 9.25, .result: 0.5},
533 {.n: 10, .rho: 0.5, .x: 6.75, .y: -2.25, .result: 0.024090},
534 {.n: 10, .rho: 1.0, .x: -1.75, .y: -1.0, .result: 0.055341},
535 {.n: 15, .rho: 0.0, .x: -1.25, .y: -4.75, .result: 0.000029},
536 {.n: 15, .rho: 0.0, .x: -2.0, .y: -1.5, .result: 0.003411},
537 {.n: 15, .rho: 0.5, .x: 3.0, .y: -3.25, .result: 0.002691},
538 {.n: 20, .rho: -0.5, .x: 2.0, .y: -1.25, .result: 0.098333},
539 {.n: 20, .rho: -1.0, .x: 3.0, .y: 8.0, .result: 0.996462},
540 {.n: 20, .rho: 0.0, .x: -7.5, .y: 1.5, .result: 0.0},
541 {.n: 20, .rho: 0.5, .x: 1.25, .y: 9.75, .result: 0.887136},
542 {.n: 25, .rho: -1.0, .x: -4.25, .y: 5.0, .result: 0.000111},
543 {.n: 25, .rho: 0.5, .x: 9.5, .y: -1.5, .result: 0.073069},
544 {.n: 25, .rho: 1.0, .x: -6.5, .y: -3.25, .result: 0.0},
545 {.n: 30, .rho: -1.0, .x: -7.75, .y: 10.0, .result: 0.0},
546 {.n: 30, .rho: 1.0, .x: 0.5, .y: 9.5, .result: 0.689638},
547 {.n: 60, .rho: -1.0, .x: -3.5, .y: -8.25, .result: 0.0},
548 {.n: 60, .rho: -1.0, .x: 4.25, .y: 0.75, .result: 0.771869},
549 {.n: 60, .rho: -1.0, .x: 5.75, .y: 3.75, .result: 0.9998},
550 {.n: 60, .rho: 0.5, .x: -4.5, .y: 8.25, .result: 0.000016},
551 {.n: 60, .rho: 1.0, .x: 6.5, .y: -4.0, .result: 0.000088},
552 {.n: 90, .rho: -0.5, .x: -3.75, .y: -2.75, .result: 0.0},
553 {.n: 90, .rho: 0.5, .x: 8.75, .y: -7.0, .result: 0.0},
554 {.n: 120, .rho: 0.0, .x: -3.5, .y: -9.25, .result: 0.0},
555 {.n: 120, .rho: 0.0, .x: -8.25, .y: 5.0, .result: 0.0},
556 {.n: 120, .rho: 1.0, .x: -0.75, .y: 3.75, .result: 0.227361},
557 {.n: 120, .rho: 1.0, .x: -3.5, .y: -8.0, .result: 0.0},
558 {.n: 150, .rho: 0.0, .x: 10.0, .y: -1.75, .result: 0.041082},
559 {.n: 300, .rho: -0.5, .x: -6.0, .y: 3.75, .result: 0.0},
560 {.n: 300, .rho: -0.5, .x: 3.5, .y: -4.5, .result: 0.000004},
561 {.n: 300, .rho: 0.0, .x: 6.5, .y: -5.0, .result: 0.0},
562 {.n: 600, .rho: -0.5, .x: 9.25, .y: 1.5, .result: 0.93293},
563 {.n: 600, .rho: -1.0, .x: -9.25, .y: 1.5, .result: 0.0},
564 {.n: 600, .rho: 0.5, .x: -5.0, .y: 8.0, .result: 0.0},
565 {.n: 600, .rho: 1.0, .x: -2.75, .y: -9.0, .result: 0.0},
566 {.n: 1000, .rho: -0.5, .x: -2.5, .y: 0.25, .result: 0.000589},
567 {.n: 1000, .rho: -0.5, .x: 3.0, .y: 1.0, .result: 0.839842},
568 {.n: 2000, .rho: -1.0, .x: 9.0, .y: -4.75, .result: 0.000001},
569 {.n: 2000, .rho: 0.5, .x: 9.75, .y: 7.25, .result: 1.0},
570 {.n: 2000, .rho: 1.0, .x: 0.75, .y: -9.0, .result: 0.0},
571 {.n: 5000, .rho: -0.5, .x: 9.75, .y: 5.5, .result: 1.0},
572 {.n: 5000, .rho: -1.0, .x: 6.0, .y: 1.0, .result: 0.841321},
573 {.n: 5000, .rho: 1.0, .x: 4.0, .y: -7.75, .result: 0.0},
574 {.n: 10000, .rho: 0.5, .x: 1.5, .y: 6.0, .result: 0.933177}
575 };
576
577 tolerance = 1.0e-6;
578 for (auto& i : cases) {
579 BivariateCumulativeStudentDistribution f(i.n, i.rho);
580 Real calculated = f(i.x, i.y);
581 Real expected = i.result;
582 if (std::fabs(x: calculated - expected) > tolerance)
583 BOOST_ERROR("Failed to reproduce CDF value:"
584 << "\n n: " << i.n << "\n rho: " << i.rho << "\n x: " << i.x
585 << "\n y: " << i.y << "\n calculated: " << calculated
586 << "\n expected: " << expected);
587 }
588}
589
590void DistributionTest::testBivariateCumulativeStudentVsBivariate() {
591 BOOST_TEST_MESSAGE(
592 "Testing bivariate cumulative Student t distribution for large N...");
593
594 Natural n = 10000; // for this value, the distribution should be
595 // close to a bivariate normal distribution.
596
597 for (Real rho = -1.0; rho < 1.01; rho += 0.25) {
598 BivariateCumulativeStudentDistribution T(n, rho);
599 BivariateCumulativeNormalDistribution N(rho);
600
601 Real avgDiff = 0.0;
602 Size m = 0;
603 Real tolerance = 4.0e-5;
604 for (Real x = -10; x < 10.1; x += 0.5) {
605 for (Real y = -10; y < 10.1; y += 0.5) {
606 Real calculated = T(x, y);
607 Real expected = N(x, y);
608 Real diff = std::fabs(x: calculated - expected);
609 if (diff > tolerance)
610 BOOST_ERROR("Failed to reproduce limit value:" <<
611 "\n rho: " << rho <<
612 "\n x: " << x <<
613 "\n y: " << y <<
614 "\n calculated: " << calculated <<
615 "\n expected: " << expected);
616
617 avgDiff += diff;
618 ++m;
619 }
620 }
621 avgDiff /= m;
622 if (avgDiff > 3.0e-6)
623 BOOST_ERROR("Failed to reproduce average limit value:" <<
624 "\n rho: " << rho <<
625 "\n average error: " << avgDiff);
626 }
627}
628
629
630namespace distributions_test {
631 class InverseNonCentralChiSquared {
632 public:
633 InverseNonCentralChiSquared(Real df, Real ncp)
634 : dist_(df, ncp) {}
635
636 Real operator()(Real x) const {
637 return boost::math::quantile(dist: dist_, p: x);
638 }
639 private:
640 const boost::math::non_central_chi_squared_distribution<Real> dist_;
641 };
642}
643
644void DistributionTest::testInvCDFviaStochasticCollocation() {
645 BOOST_TEST_MESSAGE(
646 "Testing inverse CDF based on stochastic collocation...");
647
648 using namespace distributions_test;
649
650 const Real k = 3.0;
651 const Real lambda = 1.0;
652
653 const InverseCumulativeNormal invNormalCDF;
654 const CumulativeNormalDistribution normalCDF;
655 const InverseNonCentralChiSquared invCDF(k, lambda);
656
657 const StochasticCollocationInvCDF scInvCDF10(invCDF, 10);
658
659 // low precision
660 for (Real x=-3.0; x < 3.0; x+=0.1) {
661 const Real u = normalCDF(x);
662
663 const Real calculated1 = scInvCDF10(u);
664 const Real calculated2 = scInvCDF10.value(x);
665 const Real expected = invCDF(u);
666
667 if (std::fabs(x: calculated1 - calculated2) > 1e-6) {
668 BOOST_FAIL("Failed to reproduce equal stochastic collocation "
669 "inverse CDF" <<
670 "\n x: " << x <<
671 "\n calculated via normal distribution : "
672 << calculated2 <<
673 "\n calculated via uniform distribution: "
674 << calculated1 <<
675 "\n diff: " << calculated1 - calculated2);
676 }
677
678 const Real tol = 1e-2;
679 if (std::fabs(x: calculated2 - expected) > tol) {
680 BOOST_FAIL("Failed to reproduce invCDF with "
681 "stochastic collocation method" <<
682 "\n x: " << x <<
683 "\n invCDF :" << expected <<
684 "\n scInvCDF: " << calculated2 <<
685 "\n diff : " << std::fabs(expected-calculated2) <<
686 "\n tol : " << tol);
687 }
688 }
689
690 // high precision
691 const StochasticCollocationInvCDF scInvCDF30(invCDF, 30, 0.9999999);
692 for (Real x=-4.0; x < 4.0; x+=0.1) {
693 const Real u = normalCDF(x);
694
695 const Real expected = invCDF(u);
696 const Real calculated = scInvCDF30(u);
697
698 const Real tol = 1e-6;
699 if (std::fabs(x: calculated - expected) > tol) {
700 BOOST_FAIL("Failed to reproduce invCDF with "
701 "stochastic collocation method" <<
702 "\n x: " << x <<
703 "\n invCDF :" << expected <<
704 "\n scInvCDF: " << calculated <<
705 "\n diff : " << std::fabs(expected-calculated) <<
706 "\n tol : " << tol);
707 }
708 }
709}
710
711void DistributionTest::testSankaranApproximation() {
712 BOOST_TEST_MESSAGE("Testing Sankaran approximation for the "
713 "non-central cumulative chi-square distribution...");
714
715 const Real dfs[] = {2,2,2,4,4};
716 const Real ncps[] = {1,2,3,1,2,3};
717
718 const Real tol = 0.01;
719 for (Real df : dfs) {
720 for (Real ncp : ncps) {
721 const NonCentralCumulativeChiSquareDistribution d(df, ncp);
722 const NonCentralCumulativeChiSquareSankaranApprox sankaran(df, ncp);
723
724 for (Real x=0.25; x < 10; x+=0.1) {
725 const Real expected = d(x);
726 const Real calculated = sankaran(x);
727 const Real diff = std::fabs(x: expected - calculated);
728
729 if (diff > tol) {
730 BOOST_ERROR("Failed to match accuracy of Sankaran approximation"""
731 "\n df : " << df <<
732 "\n ncp : " << ncp <<
733 "\n x : " << x <<
734 "\n expected : " << expected <<
735 "\n calculated: " << calculated <<
736 "\n diff : " << diff <<
737 "\n tol : " << tol);
738 }
739 }
740 }
741 }
742}
743
744test_suite* DistributionTest::suite(SpeedLevel) {
745 auto* suite = BOOST_TEST_SUITE("Distribution tests");
746
747 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testNormal));
748 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testBivariate));
749 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testPoisson));
750 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testCumulativePoisson));
751 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testInverseCumulativePoisson));
752 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testBivariateCumulativeStudent));
753 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testInvCDFviaStochasticCollocation));
754 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testSankaranApproximation));
755 suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testBivariateCumulativeStudentVsBivariate));
756
757 return suite;
758}
759
760

source code of quantlib/test-suite/distributions.cpp

Morty Proxy This is a proxified and sanitized view of the page, visit original site.