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_));
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)
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);
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;
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));
137 xs = std::fabs(xs*xs);
138 Real rs = std::sqrt(1 - xs);
141 return (a_ * std::exp(asr) *
142 (std::exp(-hk_ * (1 - rs) / (2 * (1 + rs))) / rs -
143 (1 +
c_ * xs * (1 +
d_ * xs))));
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);
bivariate cumulative normal distribution
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
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
ext::function< Real(Real)> b
Integral of a 1-dimensional function using the Gauss quadratures.
std::size_t Size
size of a container