QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
bivariatenormaldistribution.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2002, 2003 Ferdinando Ametrano
5 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
6 Copyright (C) 2003 StatPro Italia srl
7 Copyright (C) 2005 Gary Kennedy
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
25
26namespace QuantLib {
27
28 // Drezner 1978
29
31 0.24840615,
32 0.39233107,
33 0.21141819,
34 0.03324666,
35 0.00082485334
36 };
37
39 0.10024215,
40 0.48281397,
41 1.06094980,
42 1.77972940,
43 2.66976040000
44 };
45
48 : rho_(rho), rho2_(rho*rho) {
49
50 QL_REQUIRE(rho>=-1.0,
51 "rho must be >= -1.0 (" << rho << " not allowed)");
52 QL_REQUIRE(rho<=1.0,
53 "rho must be <= 1.0 (" << rho << " not allowed)");
54 }
55
57 Real b) const {
58
59 CumulativeNormalDistribution cumNormalDist;
60 Real CumNormDistA = cumNormalDist(a);
61 Real CumNormDistB = cumNormalDist(b);
62 Real MaxCumNormDistAB = std::max(CumNormDistA, CumNormDistB);
63 Real MinCumNormDistAB = std::min(CumNormDistA, CumNormDistB);
64
65 if (1.0-MaxCumNormDistAB<1e-15)
66 return MinCumNormDistAB;
67
68 if (MinCumNormDistAB<1e-15)
69 return MinCumNormDistAB;
70
71 Real a1 = a / std::sqrt(2.0 * (1.0 - rho2_));
72 Real b1 = b / std::sqrt(2.0 * (1.0 - rho2_));
73
74 Real result=-1.0;
75
76 if (a<=0.0 && b<=0 && rho_<=0) {
77 Real sum=0.0;
78 for (Size i=0; i<5; i++) {
79 for (Size j=0;j<5; j++) {
80 sum += x_[i]*x_[j]*
81 std::exp(a1*(2.0*y_[i]-a1)+b1*(2.0*y_[j]-b1)
82 +2.0*rho_*(y_[i]-a1)*(y_[j]-b1));
83 }
84 }
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);
97 BivariateCumulativeNormalDistributionDr78 bivCumNormalDist(rho1);
98
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);
102
103 Real delta = (1.0-(a>0.0 ? 1.0: -1.0)*(b>0.0 ? 1.0: -1.0))/4.0;
104
105 result= bivCumNormalDist(a, 0.0) + CBND2(b, 0.0) - delta;
106 } else {
107 QL_FAIL("case not handled");
108 }
109
110 return result;
111 }
112
113
114 // West 2004
115
116 namespace {
117
118 class eqn3 { /* Relates to eqn3 Genz 2004 */
119 public:
120 eqn3(Real h, Real k, Real asr)
121 : hk_(h * k), asr_(asr), hs_((h * h + k * k) / 2) {}
122
123 Real operator()(Real x) const {
124 Real sn = std::sin(asr_ * (-x + 1) * 0.5);
125 return std::exp((sn * hk_ - hs_) / (1.0 - sn * sn));
126 }
127 private:
129 };
130
131 class eqn6 { /* Relates to eqn6 Genz 2004 */
132 public:
133 eqn6(Real a, Real c, Real d, Real bs, Real hk)
134 : a_(a), c_(c), d_(d), bs_(bs), hk_(hk) {}
135 Real operator()(Real x) const {
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;
140 if (asr > -100.0) {
141 return (a_ * std::exp(asr) *
142 (std::exp(-hk_ * (1 - rs) / (2 * (1 + rs))) / rs -
143 (1 + c_ * xs * (1 + d_ * xs))));
144 } else {
145 return 0.0;
146 }
147 }
148 private:
150 };
151
152 }
153
156 : correlation_(rho) {
157
158 QL_REQUIRE(rho>=-1.0,
159 "rho must be >= -1.0 (" << rho << " not allowed)");
160 QL_REQUIRE(rho<=1.0,
161 "rho must be <= 1.0 (" << rho << " not allowed)");
162 }
163
164
166 Real x, Real y) const {
167
168 /* The implementation is described at section 2.4 "Hybrid
169 Numerical Integration Algorithms" of "Numerical Computation
170 of Rectangular Bivariate an Trivariate Normal and t
171 Probabilities", Genz (2004), Statistics and Computing 14,
172 151-160. (available at
173 www.sci.wsu.edu/math/faculty/henz/homepage)
174
175 The Gauss-Legendre quadrature have been extracted to
176 TabulatedGaussLegendre (x,w zero-based)
177
178 Tthe functions ot be integrated numerically have been moved
179 to classes eqn3 and eqn6
180
181 Change some magic numbers to M_PI */
182
183 TabulatedGaussLegendre gaussLegendreQuad(20);
184 if (std::fabs(correlation_) < 0.3) {
185 gaussLegendreQuad.order(6);
186 } else if (std::fabs(correlation_) < 0.75) {
187 gaussLegendreQuad.order(12);
188 }
189
190 Real h = -x;
191 Real k = -y;
192 Real hk = h * k;
193 Real BVN = 0.0;
194
195 if (std::fabs(correlation_) < 0.925)
196 {
197 if (std::fabs(correlation_) > 0)
198 {
199 Real asr = std::asin(correlation_);
200 eqn3 f(h,k,asr);
201 BVN = gaussLegendreQuad(f);
202 BVN *= asr * (0.25 / M_PI);
203 }
204 BVN += cumnorm_(-h) * cumnorm_(-k);
205 }
206 else
207 {
208 if (correlation_ < 0)
209 {
210 k *= -1;
211 hk *= -1;
212 }
213 if (std::fabs(correlation_) < 1)
214 {
215 Real Ass = (1 - correlation_) * (1 + correlation_);
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;
221 if (asr > -100)
222 {
223 BVN = a * std::exp(asr) *
224 (1 - c * (bs - Ass) * (1 - d * bs / 5) / 3 +
225 c * d * Ass * Ass / 5);
226 }
227 if (-hk < 100)
228 {
229 Real B = std::sqrt(bs);
230 BVN -= std::exp(-hk / 2) * 2.506628274631 *
231 cumnorm_(-B / a) * B *
232 (1 - c * bs * (1 - d * bs / 5) / 3);
233 }
234 a /= 2;
235 eqn6 f(a,c,d,bs,hk);
236 BVN += gaussLegendreQuad(f);
237 BVN /= (-2.0 * M_PI);
238 }
239
240 if (correlation_ > 0) {
241 BVN += cumnorm_(-std::max(h, k));
242 } else {
243 BVN *= -1;
244 if (k > h) {
245 // evaluate cumnorm where it is most precise, that
246 // is in the lower tail because of double accuracy
247 // around 0.0 vs around 1.0
248 if (h >= 0) {
249 BVN += cumnorm_(-h) - cumnorm_(-k);
250 } else {
251 BVN += cumnorm_(k) - cumnorm_(h);
252 }
253 }
254 }
255 }
256 return BVN;
257 }
258
259}
const Real correlation_
bivariate cumulative normal distribution
Cumulative bivariate normal distribution function.
Cumulative normal distribution function.
tabulated Gauss-Legendre quadratures
Real a_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
ext::function< Real(Real)> b
Integral of a 1-dimensional function using the Gauss quadratures.
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Real rho
#define M_PI
Definition: any.hpp:35