23#include <ql/math/distributions/bivariatenormaldistribution.hpp>
24#include <ql/math/integrals/gaussianquadratures.hpp>
48 : rho_(rho), rho2_(rho*rho) {
51 "rho must be >= -1.0 (" << rho <<
" not allowed)");
53 "rho must be <= 1.0 (" << rho <<
" not allowed)");
60 Real CumNormDistA = cumNormalDist(a);
61 Real CumNormDistB = cumNormalDist(b);
62 Real MaxCumNormDistAB = std::max(CumNormDistA, CumNormDistB);
63 Real MinCumNormDistAB = std::min(CumNormDistA, CumNormDistB);
65 if (1.0-MaxCumNormDistAB<1e-15)
66 return MinCumNormDistAB;
68 if (MinCumNormDistAB<1e-15)
69 return MinCumNormDistAB;
71 Real a1 = a / std::sqrt(2.0 * (1.0 -
rho2_));
72 Real b1 = b / std::sqrt(2.0 * (1.0 -
rho2_));
76 if (a<=0.0 && b<=0 &&
rho_<=0) {
78 for (
Size i=0; i<5; i++) {
79 for (
Size j=0;j<5; j++) {
81 std::exp(a1*(2.0*
y_[i]-a1)+b1*(2.0*
y_[j]-b1)
85 result = std::sqrt(1.0 -
rho2_)/M_PI*sum;
86 }
else if (a<=0 && b>=0 &&
rho_>=0) {
88 result= CumNormDistA - bivCumNormalDist(a, -b);
89 }
else if (a>=0.0 && b<=0.0 && rho_>=0.0) {
91 result= CumNormDistB - bivCumNormalDist(-a, b);
92 }
else if (a>=0.0 && b>=0.0 &&
rho_<=0.0) {
93 result= CumNormDistA + CumNormDistB -1.0 + (*this)(-a, -b);
94 }
else if (a*b*
rho_>0.0) {
95 Real rho1 = (
rho_*a-b)*(a>0.0 ? 1.0: -1.0)/
96 std::sqrt(a*a-2.0*
rho_*a*b+b*b);
99 Real rho2 = (
rho_*b-a)*(b>0.0 ? 1.0: -1.0)/
100 std::sqrt(a*a-2.0*
rho_*a*b+b*b);
103 Real delta = (1.0-(a>0.0 ? 1.0: -1.0)*(b>0.0 ? 1.0: -1.0))/4.0;
105 result= bivCumNormalDist(a, 0.0) + CBND2(b, 0.0) - delta;
107 QL_FAIL(
"case not handled");
121 : hk_(h * k), asr_(asr), hs_((h * h + k * k) / 2) {}
124 Real sn = std::sin(asr_ * (-x + 1) * 0.5);
125 return std::exp((sn * hk_ - hs_) / (1.0 - sn * sn));
134 : a_(a), c_(c), d_(d), bs_(bs), hk_(hk) {}
136 Real xs = a_ * (-x + 1);
137 xs = std::fabs(xs*xs);
138 Real rs = std::sqrt(1 - xs);
139 Real asr = -(bs_ / xs + hk_) / 2;
141 return (a_ * std::exp(asr) *
142 (std::exp(-hk_ * (1 - rs) / (2 * (1 + rs))) / rs -
143 (1 + c_ * xs * (1 + d_ * xs))));
149 Real a_, c_, d_, bs_, hk_;
156 : correlation_(rho) {
158 QL_REQUIRE(rho>=-1.0,
159 "rho must be >= -1.0 (" << rho <<
" not allowed)");
161 "rho must be <= 1.0 (" << rho <<
" not allowed)");
185 gaussLegendreQuad.
order(6);
187 gaussLegendreQuad.
order(12);
201 BVN = gaussLegendreQuad(f);
202 BVN *= asr * (0.25 / M_PI);
216 Real a = std::sqrt(Ass);
217 Real bs = (h-k)*(h-k);
218 Real c = (4 - hk) / 8;
219 Real d = (12 - hk) / 16;
220 Real asr = -(bs / Ass + hk) / 2;
223 BVN = a * std::exp(asr) *
224 (1 - c * (bs - Ass) * (1 - d * bs / 5) / 3 +
225 c * d * Ass * Ass / 5);
229 Real B = std::sqrt(bs);
230 BVN -= std::exp(-hk / 2) * 2.506628274631 *
232 (1 - c * bs * (1 - d * bs / 5) / 3);
236 BVN += gaussLegendreQuad(f);
237 BVN /= (-2.0 * M_PI);
Cumulative bivariate normal distribution function.
BivariateCumulativeNormalDistributionDr78(Real rho)
Real operator()(Real a, Real b) const
CumulativeNormalDistribution cumnorm_
Real operator()(Real a, Real b) const
BivariateCumulativeNormalDistributionWe04DP(Real rho)
Cumulative normal distribution function.
tabulated Gauss-Legendre quadratures
std::size_t Size
size of a container