1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
5 Copyright (C) 2002, 2003 Ferdinando Ametrano
6 Copyright (C) 2008 StatPro Italia srl
7 Copyright (C) 2010 Kakhkhor Abdijalilov
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
23#include <ql/math/distributions/normaldistribution.hpp>
24#include <ql/math/comparison.hpp>
25
26#include <boost/math/distributions/normal.hpp>
27
28namespace QuantLib {
29
30 Real CumulativeNormalDistribution::operator()(Real z) const {
31 //QL_REQUIRE(!(z >= average_ && 2.0*average_-z > average_),
32 // "not a real number. ");
33 z = (z - average_) / sigma_;
34
35 Real result = 0.5 * ( 1.0 + errorFunction_( z*M_SQRT_2 ) );
36 if (result<=1e-8) { //todo: investigate the threshold level
37 // Asymptotic expansion for very negative z following (26.2.12)
38 // on page 408 in M. Abramowitz and A. Stegun,
39 // Pocketbook of Mathematical Functions, ISBN 3-87144818-4.
40 Real sum=1.0, zsqr=z*z, i=1.0, g=1.0, x, y,
41 a=QL_MAX_REAL, lasta;
42 do {
43 lasta=a;
44 x = (4.0*i-3.0)/zsqr;
45 y = x*((4.0*i-1)/zsqr);
46 a = g*(x-y);
47 sum -= a;
48 g *= y;
49 ++i;
50 a = std::fabs(x: a);
51 } while (lasta>a && a>=std::fabs(x: sum*QL_EPSILON));
52 result = -gaussian_(z)/z*sum;
53 }
54 return result;
55 }
56
57 #if !defined(QL_PATCH_SOLARIS)
58 const CumulativeNormalDistribution InverseCumulativeNormal::f_;
59 #endif
60
61 // Coefficients for the rational approximation.
62 const Real InverseCumulativeNormal::a1_ = -3.969683028665376e+01;
63 const Real InverseCumulativeNormal::a2_ = 2.209460984245205e+02;
64 const Real InverseCumulativeNormal::a3_ = -2.759285104469687e+02;
65 const Real InverseCumulativeNormal::a4_ = 1.383577518672690e+02;
66 const Real InverseCumulativeNormal::a5_ = -3.066479806614716e+01;
67 const Real InverseCumulativeNormal::a6_ = 2.506628277459239e+00;
68
69 const Real InverseCumulativeNormal::b1_ = -5.447609879822406e+01;
70 const Real InverseCumulativeNormal::b2_ = 1.615858368580409e+02;
71 const Real InverseCumulativeNormal::b3_ = -1.556989798598866e+02;
72 const Real InverseCumulativeNormal::b4_ = 6.680131188771972e+01;
73 const Real InverseCumulativeNormal::b5_ = -1.328068155288572e+01;
74
75 const Real InverseCumulativeNormal::c1_ = -7.784894002430293e-03;
76 const Real InverseCumulativeNormal::c2_ = -3.223964580411365e-01;
77 const Real InverseCumulativeNormal::c3_ = -2.400758277161838e+00;
78 const Real InverseCumulativeNormal::c4_ = -2.549732539343734e+00;
79 const Real InverseCumulativeNormal::c5_ = 4.374664141464968e+00;
80 const Real InverseCumulativeNormal::c6_ = 2.938163982698783e+00;
81
82 const Real InverseCumulativeNormal::d1_ = 7.784695709041462e-03;
83 const Real InverseCumulativeNormal::d2_ = 3.224671290700398e-01;
84 const Real InverseCumulativeNormal::d3_ = 2.445134137142996e+00;
85 const Real InverseCumulativeNormal::d4_ = 3.754408661907416e+00;
86
87 // Limits of the approximation regions
88 const Real InverseCumulativeNormal::x_low_ = 0.02425;
89 const Real InverseCumulativeNormal::x_high_= 1.0 - x_low_;
90
91 Real InverseCumulativeNormal::tail_value(Real x) {
92 if (x <= 0.0 || x >= 1.0) {
93 // try to recover if due to numerical error
94 if (close_enough(x, y: 1.0)) {
95 return QL_MAX_REAL; // largest value available
96 } else if (std::fabs(x: x) < QL_EPSILON) {
97 return QL_MIN_REAL; // largest negative value available
98 } else {
99 QL_FAIL("InverseCumulativeNormal(" << x
100 << ") undefined: must be 0 < x < 1");
101 }
102 }
103
104 Real z;
105 if (x < x_low_) {
106 // Rational approximation for the lower region 0<x<u_low
107 z = std::sqrt(x: -2.0*std::log(x: x));
108 z = (((((c1_*z+c2_)*z+c3_)*z+c4_)*z+c5_)*z+c6_) /
109 ((((d1_*z+d2_)*z+d3_)*z+d4_)*z+1.0);
110 } else {
111 // Rational approximation for the upper region u_high<x<1
112 z = std::sqrt(x: -2.0*std::log(x: 1.0-x));
113 z = -(((((c1_*z+c2_)*z+c3_)*z+c4_)*z+c5_)*z+c6_) /
114 ((((d1_*z+d2_)*z+d3_)*z+d4_)*z+1.0);
115 }
116
117 return z;
118 }
119
120 const Real MoroInverseCumulativeNormal::a0_ = 2.50662823884;
121 const Real MoroInverseCumulativeNormal::a1_ =-18.61500062529;
122 const Real MoroInverseCumulativeNormal::a2_ = 41.39119773534;
123 const Real MoroInverseCumulativeNormal::a3_ =-25.44106049637;
124
125 const Real MoroInverseCumulativeNormal::b0_ = -8.47351093090;
126 const Real MoroInverseCumulativeNormal::b1_ = 23.08336743743;
127 const Real MoroInverseCumulativeNormal::b2_ =-21.06224101826;
128 const Real MoroInverseCumulativeNormal::b3_ = 3.13082909833;
129
130 const Real MoroInverseCumulativeNormal::c0_ = 0.3374754822726147;
131 const Real MoroInverseCumulativeNormal::c1_ = 0.9761690190917186;
132 const Real MoroInverseCumulativeNormal::c2_ = 0.1607979714918209;
133 const Real MoroInverseCumulativeNormal::c3_ = 0.0276438810333863;
134 const Real MoroInverseCumulativeNormal::c4_ = 0.0038405729373609;
135 const Real MoroInverseCumulativeNormal::c5_ = 0.0003951896511919;
136 const Real MoroInverseCumulativeNormal::c6_ = 0.0000321767881768;
137 const Real MoroInverseCumulativeNormal::c7_ = 0.0000002888167364;
138 const Real MoroInverseCumulativeNormal::c8_ = 0.0000003960315187;
139
140 Real MoroInverseCumulativeNormal::operator()(Real x) const {
141 QL_REQUIRE(x > 0.0 && x < 1.0,
142 "MoroInverseCumulativeNormal(" << x
143 << ") undefined: must be 0<x<1");
144
145 Real result;
146 Real temp=x-0.5;
147
148 if (std::fabs(x: temp) < 0.42) {
149 // Beasley and Springer, 1977
150 result=temp*temp;
151 result=temp*
152 (((a3_*result+a2_)*result+a1_)*result+a0_) /
153 ((((b3_*result+b2_)*result+b1_)*result+b0_)*result+1.0);
154 } else {
155 // improved approximation for the tail (Moro 1995)
156 if (x<0.5)
157 result = x;
158 else
159 result=1.0-x;
160 result = std::log(x: -std::log(x: result));
161 result = c0_+result*(c1_+result*(c2_+result*(c3_+result*
162 (c4_+result*(c5_+result*(c6_+result*
163 (c7_+result*c8_)))))));
164 if (x<0.5)
165 result=-result;
166 }
167
168 return average_ + result*sigma_;
169 }
170
171 MaddockInverseCumulativeNormal::MaddockInverseCumulativeNormal(
172 Real average, Real sigma)
173 : average_(average), sigma_(sigma) {}
174
175 Real MaddockInverseCumulativeNormal::operator()(Real x) const {
176 return boost::math::quantile(
177 dist: boost::math::normal_distribution<Real>(average_, sigma_), p: x);
178 }
179
180 MaddockCumulativeNormal::MaddockCumulativeNormal(
181 Real average, Real sigma)
182 : average_(average), sigma_(sigma) {}
183
184 Real MaddockCumulativeNormal::operator()(Real x) const {
185 return boost::math::cdf(
186 dist: boost::math::normal_distribution<Real>(average_, sigma_), x);
187 }
188}
189

source code of quantlib/ql/math/distributions/normaldistribution.cpp