20#include <ql/math/distributions/bivariatestudenttdistribution.hpp>
26 Real epsilon = 1.0e-8;
29 return val == 0.0 ? 0.0
30 : (val < 0.0 ? -1.0 : 1.0);
37 Real res = std::atan2(x, y);
38 return res >= 0.0 ? res : res + 2 * M_PI;
43 Real unCor = 1 - rho*rho;
44 Real sub = std::pow(h - rho * k, 2);
45 Real denom = sub + unCor * (m + k*k);
48 return sub / (sub + unCor * (m + k*k));
53 Real unCor = 1.0 - rho*rho;
55 Real div = 4 * std::sqrt(n * M_PI);
56 Real xHK = f_x(n, h, k, rho);
57 Real xKH = f_x(n, k, h, rho);
58 Real divH = 1 + h*h / n;
59 Real divK = 1 + k*k / n;
60 Real sgnHK = sign(h - rho * k);
61 Real sgnKH = sign(k - rho * h);
65 Real res = arctan(std::sqrt(unCor), -rho) / M_TWOPI;
68 Real dgM = 2 * (1 - xHK);
69 Real gjM = sgnHK * 2 / M_PI;
71 Real f_j = std::sqrt(M_PI / divK);
72 Real g_j = 1 + gjM * arctan(std::sqrt(xHK), std::sqrt(1 - xHK));
77 Real dgj = gjM * std::sqrt(xHK * (1 - xHK));
81 for (
Natural j = 3; j <= n / 2; ++j) {
82 f_j *= (j - 1.5) / (
Real) (j - 1) / divK;
83 dgj *= (
Real) (j - 2) / (2 * j - 3) * dgM;
92 gjM = sgnKH * 2 / M_PI;
94 f_j = std::sqrt(M_PI / divH);
95 g_j = 1 + gjM * arctan(std::sqrt(xKH), std::sqrt(1 - xKH));
100 Real dgj = gjM * std::sqrt(xKH * (1 - xKH));
104 for (
Natural j = 3; j <= n / 2; ++j) {
105 f_j *= (j - 1.5) / (
Real) (j - 1) / divH;
106 dgj *= (
Real) (j - 2) / (2 * j - 3) * dgM;
111 res += h / div * sum;
117 Real hkcn = hk + rho * n;
118 Real sqrtExpr = std::sqrt(h*h - 2 * rho * hk + k*k + n * unCor);
119 Real res = arctan(std::sqrt(
Real(n)) * (-(h + k) * hkcn - (hk - n) * sqrtExpr),
120 (hk - n) * hkcn - n * (h + k) * sqrtExpr ) / M_TWOPI;
124 Real mult = (1 - xHK) / 2;
126 Real f_j = 2 / std::sqrt(M_PI) / divK;
127 Real dgj = sgnHK * std::sqrt(xHK);
129 Real sum = f_j * g_j;
131 for (
Natural j = 2; j <= (n - 1) / 2; ++j) {
132 f_j *= (
Real) (j - 1) / (j - 0.5) / divK;
133 dgj *= (
Real) (2 * j - 3) / (j - 1) * mult;
137 res += k / div * sum;
140 mult = (1 - xKH) / 2;
142 f_j = 2 / std::sqrt(M_PI) / divH;
143 dgj = sgnKH * std::sqrt(xKH);
147 for (
Natural j = 2; j <= (n - 1) / 2; ++j) {
148 f_j *= (
Real) (j - 1) / (j - 0.5) / divH;
149 dgj *= (
Real) (2 * j - 3) / (j - 1) * mult;
153 res += h / div * sum;
165 : n_(n), rho_(rho) {}
169 return P_n(x, y,
n_,
rho_);
Real operator()(Real x, Real y) const
BivariateCumulativeStudentDistribution(Natural n, Real rho)
unsigned QL_INTEGER Natural
positive integer