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) 2010 Kakhkhor Abdijalilov 7 8 This file is part of QuantLib, a free-software/open-source library 9 for financial quantitative analysts and developers - http://quantlib.org/ 10 11 QuantLib is free software: you can redistribute it and/or modify it 12 under the terms of the QuantLib license. You should have received a 13 copy of the license along with this program; if not, please email 14 <quantlib-dev@lists.sf.net>. The license is also available online at 15 <http://quantlib.org/license.shtml>. 16 17 This program is distributed in the hope that it will be useful, but WITHOUT 18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 19 FOR A PARTICULAR PURPOSE. See the license for more details. 20 */ 21 22 /*! \file normaldistribution.hpp 23 \brief normal, cumulative and inverse cumulative distributions 24 */ 25 26 #ifndef quantlib_normal_distribution_hpp 27 #define quantlib_normal_distribution_hpp 28 29 #include <ql/math/errorfunction.hpp> 30 #include <ql/errors.hpp> 31 32 namespace QuantLib { 33 34 //! Normal distribution function 35 /*! Given x, it returns its probability in a Gaussian normal distribution. 36 It provides the first derivative too. 37 38 \test the correctness of the returned value is tested by 39 checking it against numerical calculations. Cross-checks 40 are also performed against the 41 CumulativeNormalDistribution and InverseCumulativeNormal 42 classes. 43 */ 44 class NormalDistribution { 45 public: 46 typedef Real argument_type; 47 typedef Real result_type; 48 49 NormalDistribution(Real average = 0.0, 50 Real sigma = 1.0); 51 // function 52 Real operator()(Real x) const; 53 Real derivative(Real x) const; 54 private: 55 Real average_, sigma_, normalizationFactor_, denominator_, 56 derNormalizationFactor_; 57 }; 58 59 typedef NormalDistribution GaussianDistribution; 60 61 62 //! Cumulative normal distribution function 63 /*! Given x it provides an approximation to the 64 integral of the gaussian normal distribution: 65 formula here ... 66 67 For this implementation see M. Abramowitz and I. Stegun, 68 Handbook of Mathematical Functions, 69 Dover Publications, New York (1972) 70 */ 71 class CumulativeNormalDistribution { 72 public: 73 typedef Real argument_type; 74 typedef Real result_type; 75 76 CumulativeNormalDistribution(Real average = 0.0, 77 Real sigma = 1.0); 78 // function 79 Real operator()(Real x) const; 80 Real derivative(Real x) const; 81 private: 82 Real average_, sigma_; 83 NormalDistribution gaussian_; 84 ErrorFunction errorFunction_; 85 }; 86 87 88 //! Inverse cumulative normal distribution function 89 /*! Given x between zero and one as 90 the integral value of a gaussian normal distribution 91 this class provides the value y such that 92 formula here ... 93 94 It use Acklam's approximation: 95 by Peter J. Acklam, University of Oslo, Statistics Division. 96 URL: http://home.online.no/~pjacklam/notes/invnorm/index.html 97 98 This class can also be used to generate a gaussian normal 99 distribution from a uniform distribution. 100 This is especially useful when a gaussian normal distribution 101 is generated from a low discrepancy uniform distribution: 102 in this case the traditional Box-Muller approach and its 103 variants would not preserve the sequence's low-discrepancy. 104 105 */ 106 class InverseCumulativeNormal { 107 public: 108 typedef Real argument_type; 109 typedef Real result_type; 110 111 InverseCumulativeNormal(Real average = 0.0, 112 Real sigma = 1.0); 113 // function operator ()(Real x) const114 Real operator()(Real x) const { 115 return average_ + sigma_*standard_value(x); 116 } 117 // value for average=0, sigma=1 118 /* Compared to operator(), this method avoids 2 floating point 119 operations (we use average=0 and sigma=1 most of the 120 time). The speed difference is noticeable. 121 */ standard_value(Real x)122 static Real standard_value(Real x) { 123 Real z; 124 if (x < x_low_ || x_high_ < x) { 125 z = tail_value(x); 126 } else { 127 z = x - 0.5; 128 Real r = z*z; 129 z = (((((a1_*r+a2_)*r+a3_)*r+a4_)*r+a5_)*r+a6_)*z / 130 (((((b1_*r+b2_)*r+b3_)*r+b4_)*r+b5_)*r+1.0); 131 } 132 133 // The relative error of the approximation has absolute value less 134 // than 1.15e-9. One iteration of Halley's rational method (third 135 // order) gives full machine precision. 136 // #define REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD 137 #ifdef REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD 138 // error (f_(z) - x) divided by the cumulative's derivative 139 const Real r = (f_(z) - x) * M_SQRT2 * M_SQRTPI * exp(0.5 * z*z); 140 // Halley's method 141 z -= r/(1+0.5*z*r); 142 #endif 143 144 return z; 145 } 146 private: 147 /* Handling tails moved into a separate method, which should 148 make the inlining of operator() and standard_value method 149 easier. tail_value is called rarely and doesn't need to be 150 inlined. 151 */ 152 static Real tail_value(Real x); 153 #if defined(QL_PATCH_SOLARIS) 154 CumulativeNormalDistribution f_; 155 #else 156 static const CumulativeNormalDistribution f_; 157 #endif 158 Real average_, sigma_; 159 static const Real a1_; 160 static const Real a2_; 161 static const Real a3_; 162 static const Real a4_; 163 static const Real a5_; 164 static const Real a6_; 165 static const Real b1_; 166 static const Real b2_; 167 static const Real b3_; 168 static const Real b4_; 169 static const Real b5_; 170 static const Real c1_; 171 static const Real c2_; 172 static const Real c3_; 173 static const Real c4_; 174 static const Real c5_; 175 static const Real c6_; 176 static const Real d1_; 177 static const Real d2_; 178 static const Real d3_; 179 static const Real d4_; 180 static const Real x_low_; 181 static const Real x_high_; 182 }; 183 184 // backward compatibility 185 typedef InverseCumulativeNormal InvCumulativeNormalDistribution; 186 187 //! Moro Inverse cumulative normal distribution class 188 /*! Given x between zero and one as 189 the integral value of a gaussian normal distribution 190 this class provides the value y such that 191 formula here ... 192 193 It uses Beasly and Springer approximation, with an improved 194 approximation for the tails. See Boris Moro, 195 "The Full Monte", 1995, Risk Magazine. 196 197 This class can also be used to generate a gaussian normal 198 distribution from a uniform distribution. 199 This is especially useful when a gaussian normal distribution 200 is generated from a low discrepancy uniform distribution: 201 in this case the traditional Box-Muller approach and its 202 variants would not preserve the sequence's low-discrepancy. 203 204 Peter J. Acklam's approximation is better and is available 205 as QuantLib::InverseCumulativeNormal 206 */ 207 class MoroInverseCumulativeNormal { 208 public: 209 typedef Real argument_type; 210 typedef Real result_type; 211 212 MoroInverseCumulativeNormal(Real average = 0.0, 213 Real sigma = 1.0); 214 // function 215 Real operator()(Real x) const; 216 private: 217 Real average_, sigma_; 218 static const Real a0_; 219 static const Real a1_; 220 static const Real a2_; 221 static const Real a3_; 222 static const Real b0_; 223 static const Real b1_; 224 static const Real b2_; 225 static const Real b3_; 226 static const Real c0_; 227 static const Real c1_; 228 static const Real c2_; 229 static const Real c3_; 230 static const Real c4_; 231 static const Real c5_; 232 static const Real c6_; 233 static const Real c7_; 234 static const Real c8_; 235 }; 236 237 //! Maddock's Inverse cumulative normal distribution class 238 /*! Given x between zero and one as 239 the integral value of a gaussian normal distribution 240 this class provides the value y such that 241 formula here ... 242 243 From the boost documentation: 244 These functions use a rational approximation devised by 245 John Maddock to calculate an initial approximation to the 246 result that is accurate to ~10^-19, then only if that has 247 insufficient accuracy compared to the epsilon for type double, 248 do we clean up the result using Halley iteration. 249 */ 250 class MaddockInverseCumulativeNormal { 251 public: 252 typedef Real argument_type; 253 typedef Real result_type; 254 255 MaddockInverseCumulativeNormal(Real average = 0.0, 256 Real sigma = 1.0); 257 Real operator()(Real x) const; 258 259 private: 260 const Real average_, sigma_; 261 }; 262 263 //! Maddock's cumulative normal distribution class 264 class MaddockCumulativeNormal { 265 public: 266 typedef Real argument_type; 267 typedef Real result_type; 268 269 MaddockCumulativeNormal(Real average = 0.0, 270 Real sigma = 1.0); 271 Real operator()(Real x) const; 272 273 private: 274 const Real average_, sigma_; 275 }; 276 277 278 // inline definitions 279 NormalDistribution(Real average,Real sigma)280 inline NormalDistribution::NormalDistribution(Real average, 281 Real sigma) 282 : average_(average), sigma_(sigma) { 283 284 QL_REQUIRE(sigma_>0.0, 285 "sigma must be greater than 0.0 (" 286 << sigma_ << " not allowed)"); 287 288 normalizationFactor_ = M_SQRT_2*M_1_SQRTPI/sigma_; 289 derNormalizationFactor_ = sigma_*sigma_; 290 denominator_ = 2.0*derNormalizationFactor_; 291 } 292 operator ()(Real x) const293 inline Real NormalDistribution::operator()(Real x) const { 294 Real deltax = x-average_; 295 Real exponent = -(deltax*deltax)/denominator_; 296 // debian alpha had some strange problem in the very-low range 297 return exponent <= -690.0 ? 0.0 : // exp(x) < 1.0e-300 anyway 298 normalizationFactor_*std::exp(exponent); 299 } 300 derivative(Real x) const301 inline Real NormalDistribution::derivative(Real x) const { 302 return ((*this)(x) * (average_ - x)) / derNormalizationFactor_; 303 } 304 CumulativeNormalDistribution(Real average,Real sigma)305 inline CumulativeNormalDistribution::CumulativeNormalDistribution( 306 Real average, Real sigma) 307 : average_(average), sigma_(sigma) { 308 309 QL_REQUIRE(sigma_>0.0, 310 "sigma must be greater than 0.0 (" 311 << sigma_ << " not allowed)"); 312 } 313 derivative(Real x) const314 inline Real CumulativeNormalDistribution::derivative(Real x) const { 315 Real xn = (x - average_) / sigma_; 316 return gaussian_(xn) / sigma_; 317 } 318 InverseCumulativeNormal(Real average,Real sigma)319 inline InverseCumulativeNormal::InverseCumulativeNormal( 320 Real average, Real sigma) 321 : average_(average), sigma_(sigma) { 322 323 QL_REQUIRE(sigma_>0.0, 324 "sigma must be greater than 0.0 (" 325 << sigma_ << " not allowed)"); 326 } 327 MoroInverseCumulativeNormal(Real average,Real sigma)328 inline MoroInverseCumulativeNormal::MoroInverseCumulativeNormal( 329 Real average, Real sigma) 330 : average_(average), sigma_(sigma) { 331 332 QL_REQUIRE(sigma_>0.0, 333 "sigma must be greater than 0.0 (" 334 << sigma_ << " not allowed)"); 335 } 336 337 } 338 339 340 #endif 341