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