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