Use Boost for inverse cumulative Student t distribution

[mpwaser]

Proposal

This improves InverseCumulativeStudent by replacing the current unbracketed Newton iteration with Boost.Math's Student t quantile implementation.

The existing implementation already contains a TODO asking to replace the Newton iteration with a more efficient inverse Student t algorithm:

QuantLib/ql/math/distributions/studenttdistribution.hpp

Lines 75 to 78 in 4567efa

	/*! \todo Find/implement an efficient algorithm for evaluating the
	cumulative Student t-distribution, replacing the Newton
	iteration
	*/

This PR proposes to replace the algorithm with a Boost.Math based implementation:

• Boost.Math implementation, this PR: https://github.com/mpwaser/quant-lib/tree/feature/student-t-inverse-cdf-boost
• Boost.Math diff against QuantLib master: master...mpwaser:feature/student-t-inverse-cdf-boost

The proposed Boost.Math based implementation:

• fixes failures and large tail errors in the current Newton implementation
• is substantially faster in the tested grids
• keeps the existing Student t input range policy (see implementation note at the end)

In case you prefer not to delegate this calculation to Boost.Math, I also prepared a native QuantLib-only Brent-based alternative. It is also a clear improvement over the current Newton implementation, though the Boost.Math implementation performs better in the benchmarks below.

• Brent alternative: https://github.com/mpwaser/quant-lib/tree/feature/student-t-inverse-cdf
• Brent diff against QuantLib master: master...mpwaser:feature/student-t-inverse-cdf

The following tables compare all three implementations in terms of accuracy, robustness, and efficiency.

Benchmark results

Accuracy and robustness

Absolute and relative errors are measured against Boost.Math quantile. Since this PR delegates to Boost.Math, its absolute and relative errors are zero by construction in this comparison. The CDF residual is included as an additional check using QuantLib's cumulative Student t distribution.

Method	OK	Failed	max abs error	max rel error	max CDF residual
Newton (current)	200	25	3.180619e+14	9.984223e-01	9.773293e-13
Brent	225	0	1.396984e-09	1.878553e-12	2.997602e-14
Boost.Math	225	0	0.000000e+00	0.000000e+00	3.140960e-13

Timing (ns/call)

Scenario	Newton (old)	Brent	Boost.Math
Central	886.4	655.4	73.7
Common	1255.3	874.6	81.7
Moderate tails	2621.6	1309.2	99.5
Mixed deep tails	3406.1	1199.8	97.4

Validation

• ql_test_suite --run_test=QuantLibTests/DistributionTests/testInverseCumulativeStudent
• ql_test_suite --run_test=QuantLibTests/DistributionTests
• git diff --check

Implementation note: Endpoint handling

The previous Student t inverse accepted inputs in [0, 1]. At the exact endpoints, the mathematical inverse is not finite. The old Newton iteration therefore returned finite, accuracy-dependent tail values as an artifact of its stopping criterion.

This implementation handles exact and numerically near-endpoint inputs before calling Boost.Math, following the convention already used by QuantLib's InverseCumulativeNormal::tail_value:

• inputs close to 0 return QL_MIN_REAL
• inputs close to 1 return QL_MAX_REAL
• inputs genuinely outside [0, 1] still throw

This avoids exposing Boost.Math's endpoint overflow behavior while aligning the Student t inverse with an existing QuantLib inverse-distribution convention.

One detail worth noting is that the inherited normal-inverse convention is slightly asymmetric: the upper endpoint is detected with close_enough(y, 1.0), while the lower endpoint is detected with std::fabs(y) < QL_EPSILON. I kept that behavior here for consistency with InverseCumulativeNormal.

A cleaner alternative would be a symmetric endpoint policy, for instance using the same style of tolerance check at both endpoints. I did not include that broader policy change in this PR in order to keep the change focused on replacing the Student t Newton iteration.

Appendix: Benchmark details

The benchmark was run outside the PR code, against a Release build with -O3 -DNDEBUG.

The implementations compared were:

• current QuantLib Newton implementation, copied into the benchmark
• native QuantLib-only Brent alternative
• this PR's InverseCumulativeStudent implementation using Boost.Math
• direct Boost.Math quantile with the distribution object cached per degree of freedom

For the QuantLib implementations, the benchmark used accuracy = 1e-12 and maxIterations = 100.

The accuracy grid used degrees of freedom: 1, 2, 3, 4, 5, 10, 30, 100, 1000.

The accuracy grid used probabilities: 1e-15, 1e-12, 1e-10, 1e-9, 1e-8, 1e-6, 1e-4, 0.001, 0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99, 0.999, 1-1e-4, 1-1e-6, 1-1e-8, 1-1e-9, 1-1e-10, 1-1e-12, 1-1e-15.

In the accuracy table, Failed means that the implementation either threw or returned a non-finite value for the tested input. In this benchmark, the 25 old-Newton failures were all exceptions for n = 1000, caused by non-finite intermediate values in the Newton iteration. This high number of degrees of freedom was included as a stress case to check robustness outside the common parameter range. The large tail errors shown separately are finite results, not counted as failures.

Timing used std::chrono::steady_clock; each row reports nanoseconds per call. A volatile checksum was accumulated to prevent calls from being optimized away. Direct Boost.Math with the distribution object cached per degree of freedom gave very similar timings to this PR implementation, so distribution construction was not a material cost in this benchmark.

[CLAassistant]

All committers have signed the CLA.

[boring-cyborg]

Thanks for opening this pull request! It might take a while before we look at it, so don't worry if there seems to be no feedback. We'll get to it.

[coveralls]

coverage: 74.954% (+0.07%) from 74.885% — mpwaser:feature/student-t-inverse-cdf-boost into lballabio:master