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
32namespace 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 /*! \deprecated Use `auto` or `decltype` instead.
47 Deprecated in version 1.29.
48 */
49 QL_DEPRECATED
50 typedef Real argument_type;
51
52 /*! \deprecated Use `auto` or `decltype` instead.
53 Deprecated in version 1.29.
54 */
55 QL_DEPRECATED
56 typedef Real result_type;
57
58 NormalDistribution(Real average = 0.0,
59 Real sigma = 1.0);
60 // function
61 Real operator()(Real x) const;
62 Real derivative(Real x) const;
63 private:
64 Real average_, sigma_, normalizationFactor_, denominator_,
65 derNormalizationFactor_;
66 };
67
68 typedef NormalDistribution GaussianDistribution;
69
70
71 //! Cumulative normal distribution function
72 /*! Given x it provides an approximation to the
73 integral of the gaussian normal distribution:
74 formula here ...
75
76 For this implementation see M. Abramowitz and I. Stegun,
77 Handbook of Mathematical Functions,
78 Dover Publications, New York (1972)
79 */
80 class CumulativeNormalDistribution {
81 public:
82 /*! \deprecated Use `auto` or `decltype` instead.
83 Deprecated in version 1.29.
84 */
85 QL_DEPRECATED
86 typedef Real argument_type;
87
88 /*! \deprecated Use `auto` or `decltype` instead.
89 Deprecated in version 1.29.
90 */
91 QL_DEPRECATED
92 typedef Real result_type;
93
94 CumulativeNormalDistribution(Real average = 0.0,
95 Real sigma = 1.0);
96 // function
97 Real operator()(Real x) const;
98 Real derivative(Real x) const;
99 private:
100 Real average_, sigma_;
101 NormalDistribution gaussian_;
102 ErrorFunction errorFunction_;
103 };
104
105
106 //! Inverse cumulative normal distribution function
107 /*! Given x between zero and one as
108 the integral value of a gaussian normal distribution
109 this class provides the value y such that
110 formula here ...
111
112 It use Acklam's approximation:
113 by Peter J. Acklam, University of Oslo, Statistics Division.
114 URL: http://home.online.no/~pjacklam/notes/invnorm/index.html
115
116 This class can also be used to generate a gaussian normal
117 distribution from a uniform distribution.
118 This is especially useful when a gaussian normal distribution
119 is generated from a low discrepancy uniform distribution:
120 in this case the traditional Box-Muller approach and its
121 variants would not preserve the sequence's low-discrepancy.
122
123 */
124 class InverseCumulativeNormal {
125 public:
126 /*! \deprecated Use `auto` or `decltype` instead.
127 Deprecated in version 1.29.
128 */
129 QL_DEPRECATED
130 typedef Real argument_type;
131
132 /*! \deprecated Use `auto` or `decltype` instead.
133 Deprecated in version 1.29.
134 */
135 QL_DEPRECATED
136 typedef Real result_type;
137
138 InverseCumulativeNormal(Real average = 0.0,
139 Real sigma = 1.0);
140 // function
141 Real operator()(Real x) const {
142 return average_ + sigma_*standard_value(x);
143 }
144 // value for average=0, sigma=1
145 /* Compared to operator(), this method avoids 2 floating point
146 operations (we use average=0 and sigma=1 most of the
147 time). The speed difference is noticeable.
148 */
149 static Real standard_value(Real x) {
150 Real z;
151 if (x < x_low_ || x_high_ < x) {
152 z = tail_value(x);
153 } else {
154 z = x - 0.5;
155 Real r = z*z;
156 z = (((((a1_*r+a2_)*r+a3_)*r+a4_)*r+a5_)*r+a6_)*z /
157 (((((b1_*r+b2_)*r+b3_)*r+b4_)*r+b5_)*r+1.0);
158 }
159
160 // The relative error of the approximation has absolute value less
161 // than 1.15e-9. One iteration of Halley's rational method (third
162 // order) gives full machine precision.
163 // #define REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD
164 #ifdef REFINE_TO_FULL_MACHINE_PRECISION_USING_HALLEYS_METHOD
165 // error (f_(z) - x) divided by the cumulative's derivative
166 const Real r = (f_(z) - x) * M_SQRT2 * M_SQRTPI * exp(0.5 * z*z);
167 // Halley's method
168 z -= r/(1+0.5*z*r);
169 #endif
170
171 return z;
172 }
173 private:
174 /* Handling tails moved into a separate method, which should
175 make the inlining of operator() and standard_value method
176 easier. tail_value is called rarely and doesn't need to be
177 inlined.
178 */
179 static Real tail_value(Real x);
180 #if defined(QL_PATCH_SOLARIS)
181 CumulativeNormalDistribution f_;
182 #else
183 static const CumulativeNormalDistribution f_;
184 #endif
185 Real average_, sigma_;
186 static const Real a1_;
187 static const Real a2_;
188 static const Real a3_;
189 static const Real a4_;
190 static const Real a5_;
191 static const Real a6_;
192 static const Real b1_;
193 static const Real b2_;
194 static const Real b3_;
195 static const Real b4_;
196 static const Real b5_;
197 static const Real c1_;
198 static const Real c2_;
199 static const Real c3_;
200 static const Real c4_;
201 static const Real c5_;
202 static const Real c6_;
203 static const Real d1_;
204 static const Real d2_;
205 static const Real d3_;
206 static const Real d4_;
207 static const Real x_low_;
208 static const Real x_high_;
209 };
210
211 // backward compatibility
212 typedef InverseCumulativeNormal InvCumulativeNormalDistribution;
213
214 //! Moro Inverse cumulative normal distribution class
215 /*! Given x between zero and one as
216 the integral value of a gaussian normal distribution
217 this class provides the value y such that
218 formula here ...
219
220 It uses Beasly and Springer approximation, with an improved
221 approximation for the tails. See Boris Moro,
222 "The Full Monte", 1995, Risk Magazine.
223
224 This class can also be used to generate a gaussian normal
225 distribution from a uniform distribution.
226 This is especially useful when a gaussian normal distribution
227 is generated from a low discrepancy uniform distribution:
228 in this case the traditional Box-Muller approach and its
229 variants would not preserve the sequence's low-discrepancy.
230
231 Peter J. Acklam's approximation is better and is available
232 as QuantLib::InverseCumulativeNormal
233 */
234 class MoroInverseCumulativeNormal {
235 public:
236 /*! \deprecated Use `auto` or `decltype` instead.
237 Deprecated in version 1.29.
238 */
239 QL_DEPRECATED
240 typedef Real argument_type;
241
242 /*! \deprecated Use `auto` or `decltype` instead.
243 Deprecated in version 1.29.
244 */
245 QL_DEPRECATED
246 typedef Real result_type;
247
248 MoroInverseCumulativeNormal(Real average = 0.0,
249 Real sigma = 1.0);
250 // function
251 Real operator()(Real x) const;
252 private:
253 Real average_, sigma_;
254 static const Real a0_;
255 static const Real a1_;
256 static const Real a2_;
257 static const Real a3_;
258 static const Real b0_;
259 static const Real b1_;
260 static const Real b2_;
261 static const Real b3_;
262 static const Real c0_;
263 static const Real c1_;
264 static const Real c2_;
265 static const Real c3_;
266 static const Real c4_;
267 static const Real c5_;
268 static const Real c6_;
269 static const Real c7_;
270 static const Real c8_;
271 };
272
273 //! Maddock's Inverse cumulative normal distribution class
274 /*! Given x between zero and one as
275 the integral value of a gaussian normal distribution
276 this class provides the value y such that
277 formula here ...
278
279 From the boost documentation:
280 These functions use a rational approximation devised by
281 John Maddock to calculate an initial approximation to the
282 result that is accurate to ~10^-19, then only if that has
283 insufficient accuracy compared to the epsilon for type double,
284 do we clean up the result using Halley iteration.
285 */
286 class MaddockInverseCumulativeNormal {
287 public:
288 /*! \deprecated Use `auto` or `decltype` instead.
289 Deprecated in version 1.29.
290 */
291 QL_DEPRECATED
292 typedef Real argument_type;
293
294 /*! \deprecated Use `auto` or `decltype` instead.
295 Deprecated in version 1.29.
296 */
297 QL_DEPRECATED
298 typedef Real result_type;
299
300 MaddockInverseCumulativeNormal(Real average = 0.0,
301 Real sigma = 1.0);
302 Real operator()(Real x) const;
303
304 private:
305 const Real average_, sigma_;
306 };
307
308 //! Maddock's cumulative normal distribution class
309 class MaddockCumulativeNormal {
310 public:
311 /*! \deprecated Use `auto` or `decltype` instead.
312 Deprecated in version 1.29.
313 */
314 QL_DEPRECATED
315 typedef Real argument_type;
316
317 /*! \deprecated Use `auto` or `decltype` instead.
318 Deprecated in version 1.29.
319 */
320 QL_DEPRECATED
321 typedef Real result_type;
322
323 MaddockCumulativeNormal(Real average = 0.0,
324 Real sigma = 1.0);
325 Real operator()(Real x) const;
326
327 private:
328 const Real average_, sigma_;
329 };
330
331
332 // inline definitions
333
334 inline NormalDistribution::NormalDistribution(Real average,
335 Real sigma)
336 : average_(average), sigma_(sigma) {
337
338 QL_REQUIRE(sigma_>0.0,
339 "sigma must be greater than 0.0 ("
340 << sigma_ << " not allowed)");
341
342 normalizationFactor_ = M_SQRT_2*M_1_SQRTPI/sigma_;
343 derNormalizationFactor_ = sigma_*sigma_;
344 denominator_ = 2.0*derNormalizationFactor_;
345 }
346
347 inline Real NormalDistribution::operator()(Real x) const {
348 Real deltax = x-average_;
349 Real exponent = -(deltax*deltax)/denominator_;
350 // debian alpha had some strange problem in the very-low range
351 return exponent <= -690.0 ? 0.0 : // exp(x) < 1.0e-300 anyway
352 Real(normalizationFactor_*std::exp(x: exponent));
353 }
354
355 inline Real NormalDistribution::derivative(Real x) const {
356 return ((*this)(x) * (average_ - x)) / derNormalizationFactor_;
357 }
358
359 inline CumulativeNormalDistribution::CumulativeNormalDistribution(
360 Real average, Real sigma)
361 : average_(average), sigma_(sigma) {
362
363 QL_REQUIRE(sigma_>0.0,
364 "sigma must be greater than 0.0 ("
365 << sigma_ << " not allowed)");
366 }
367
368 inline Real CumulativeNormalDistribution::derivative(Real x) const {
369 Real xn = (x - average_) / sigma_;
370 return gaussian_(xn) / sigma_;
371 }
372
373 inline InverseCumulativeNormal::InverseCumulativeNormal(
374 Real average, Real sigma)
375 : average_(average), sigma_(sigma) {
376
377 QL_REQUIRE(sigma_>0.0,
378 "sigma must be greater than 0.0 ("
379 << sigma_ << " not allowed)");
380 }
381
382 inline MoroInverseCumulativeNormal::MoroInverseCumulativeNormal(
383 Real average, Real sigma)
384 : average_(average), sigma_(sigma) {
385
386 QL_REQUIRE(sigma_>0.0,
387 "sigma must be greater than 0.0 ("
388 << sigma_ << " not allowed)");
389 }
390
391}
392
393
394#endif
395

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