QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
lossdistribution.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Roland Lichters
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
22
23using namespace std;
24
25namespace QuantLib {
26
27 //--------------------------------------------------------------------------
29 //--------------------------------------------------------------------------
30 BinomialDistribution binomial (p[0], p.size());
31 return binomial(n);
32 }
33
34 //--------------------------------------------------------------------------
36 //--------------------------------------------------------------------------
37 CumulativeBinomialDistribution binomial(p[0], p.size());
38 return 1.0 - binomial(n-1);
39 /*
40 Real defp = 0;
41 for (Size i = n; i <= p.size(); i++)
42 defp += binomialProbabilityOfNEvents (i, p);
43
44 return defp;
45 */
46 }
47
48 //--------------------------------------------------------------------------
49 vector<Real> LossDist::probabilityOfNEvents(vector<Real>& p) {
50 //--------------------------------------------------------------------------
51 Size n = p.size();
52 vector<Real> probability(n+1, 0.0);
53 vector<Real> prev;
54 probability[0] = 1.0;
55 for (Size j = 0; j < n; j++) {
56 prev = probability;
57 probability[0] = prev[0] * (1.0 - p[j]);
58 for (Size i = 1; i <= j; i++)
59 probability[i] = prev[i-1] * p[j] + prev[i] * (1.0 - p[j]);
60 probability[j+1] = prev[j] * p[j];
61 }
62
63 return probability;
64 }
65
66 //--------------------------------------------------------------------------
67 Real LossDist::probabilityOfNEvents(int k, vector<Real>& p) {
68 //--------------------------------------------------------------------------
69 return probabilityOfNEvents(p)[k];
70
71// vector<Real> w (p.size(), 0);
72// vector<Real> u (k+1, 0);
73// vector<Real> v (k+1, 0);
74
75// Real pZero = 1.0;
76// for (Size i = 0; i < w.size(); i++) {
77// pZero *= (1.0 - p[i]);
78// w[i] = p[i] / (1.0 - p[i]);
79// }
80
81// if (k == 0) return pZero;
82
83// int kk = k;
84// Real prodw = 1.0;
85
86// Cumulated probability of up to n events:
87// Cut off when the cumulated probability reaches 1,
88// i.e. set all following probabilities of exactly n events to zero.
89// Real sum = 1.0;
90
91// u[0] = 1.0;
92// for (int i = 1; i <= kk; i++) {
93// v[i] = 0;
94// for (Size j = 0; j < w.size(); j++)
95// v[i] += pow (w[j], i);
96// u[i] = 0;
97// for (int j = 1; j <= i; j++)
98// u[i] += pow (-1.0, j+1) * v[j] * u[i-j];
99// u[i] /= i;
100
101// cut off
102// if (sum * pZero >= 1.0 || u[i] < 0 || u[i] * pZero >= 1.0)
103// u[i] = 0;
104
105// sum += u[i];
106// }
107
108// return pZero * prodw * u[kk];
109 }
110
111 //--------------------------------------------------------------------------
113 //--------------------------------------------------------------------------
114 vector<Real> probability = probabilityOfNEvents(p);
115 Real sum = 1.0;
116 for (int j = 0; j < k; j++)
117 sum -= probability[j];
118 return sum;
119 /*
120 Real sum = 0;
121 for (Size i = k; i <= p.size(); i++)
122 sum += probabilityOfNEvents (i, p);
123 return sum;
124 */
125 }
126
127 //--------------------------------------------------------------------------
129 //--------------------------------------------------------------------------
131 }
132
133 //--------------------------------------------------------------------------
135 //--------------------------------------------------------------------------
137 }
138
139 //--------------------------------------------------------------------------
141 //--------------------------------------------------------------------------
143 }
144
145 //--------------------------------------------------------------------------
147 Real probability) const {
148 //--------------------------------------------------------------------------
149 n_ = n;
150 probability_.clear();
151 probability_.resize(n_+1, 0.0);
152 Distribution dist (nBuckets_, 0.0, maximum_);
154 for (Size i = 0; i <= n; i++) {
155 if (volume_ * i <= maximum_) {
156 probability_[i] = binomial(i);
157 Size bucket = dist.locate(volume * i);
158 dist.addDensity (bucket, probability_[i] / dist.dx(bucket));
159 dist.addAverage (bucket, volume * i);
160 }
161 }
162
163 excessProbability_.clear();
164 excessProbability_.resize(n_+1, 0.0);
166 for (int k = n_-1; k >= 0; k--)
168
169 dist.normalize();
170
171 return dist;
172 }
173
174 //--------------------------------------------------------------------------
175 Distribution LossDistBinomial::operator()(const vector<Real>& nominals,
176 const vector<Real>& probabilities) const {
177 //--------------------------------------------------------------------------
178 return operator()(nominals.size(), nominals[0], probabilities[0]);
179 }
180
181 //--------------------------------------------------------------------------
183 const vector<Real>& p) const {
184 //--------------------------------------------------------------------------
185 volume_ = volume;
186 n_ = p.size();
187 probability_.clear();
188 probability_.resize(n_+1, 0.0);
189 vector<Real> prev;
190 probability_[0] = 1.0;
191 for (Size k = 0; k < n_; k++) {
192 prev = probability_;
193 probability_[0] = prev[0] * (1.0 - p[k]);
194 for (Size i = 1; i <= k; i++)
195 probability_[i] = prev[i-1] * p[k] + prev[i] * (1.0 - p[k]);
196 probability_[k+1] = prev[k] * p[k];
197 }
198
199 excessProbability_.clear();
200 excessProbability_.resize(n_+1, 0.0);
202 for (int k = n_ - 1; k >= 0; k--)
204
205 Distribution dist (nBuckets_, 0.0, maximum_);
206 for (Size i = 0; i <= n_; i++) {
207 if (volume * i <= maximum_) {
208 Size bucket = dist.locate(volume * i);
209 dist.addDensity (bucket, probability_[i] / dist.dx(bucket));
210 dist.addAverage (bucket, volume*i);
211 }
212 }
213
214 dist.normalize();
215
216 return dist;
217 }
218
219 //--------------------------------------------------------------------------
220 Distribution LossDistHomogeneous::operator()(const vector<Real>& nominals,
221 const vector<Real>& probabilities) const {
222 //--------------------------------------------------------------------------
223 return operator()(nominals[0], probabilities);
224 }
225
226 //--------------------------------------------------------------------------
227 Distribution LossDistBucketing::operator()(const vector<Real>& nominals,
228 const vector<Real>& probabilities) const {
229 //--------------------------------------------------------------------------
230 QL_REQUIRE (nominals.size() == probabilities.size(), "sizes differ: "
231 << nominals.size() << " vs " << probabilities.size());
232
233 vector<Real> p (nBuckets_, 0.0);
234 vector<Real> a (nBuckets_, 0.0);
235 vector<Real> ap (nBuckets_, 0.0);
236
237 p[0] = 1.0;
238 a[0] = 0.0;
239 Real dx = maximum_ / nBuckets_;
240 for (Size k = 1; k < nBuckets_; k++)
241 a[k] = dx * k + dx/2;
242
243 for (Size i = 0; i < nominals.size(); i++) {
244 Real L = nominals[i];
245 Real P = probabilities[i];
246 for (int k = a.size()-1; k >= 0; k--) {
247 if (p[k] > 0) {
248 int u = locateTargetBucket (a[k] + L, k);
249 QL_REQUIRE (u >= 0, "u=" << u << " at i=" << i << " k=" << k);
250 QL_REQUIRE (u >= k, "u=" << u << "<k=" << k << " at i=" << i);
251
252 Real dp = p[k] * P;
253 if (u == k)
254 a[k] += P * L;
255 else {
256 // no update of a[u] and p[u] if u is beyond grid end
257 if (u < int(nBuckets_)) {
258 // a[u] remains unchanged, if dp = 0
259 if (dp > 0.0) {
260 // on Windows, p[u]/dp could cause a NaN for
261 // some very small values of p[k].
262 // Writing the above as (p[u]/p[k])/P prevents
263 // the NaN. What can I say?
264 Real f = 1.0 / (1.0 + (p[u]/p[k]) / P);
265 a[u] = (1.0 - f) * a[u] + f * (a[k] + L);
266 }
267 /* formulation of Hull-White:
268 if (p[u] + dp > 0)
269 a[u] = (p[u] * a[u] + dp * (a[k] + L))
270 / (p[u] + dp);
271 */
272 p[u] += dp;
273 }
274 p[k] -= dp;
275 }
276 }
277 QL_REQUIRE(a[k] + epsilon_ >= dx * k && a[k] < dx * (k+1),
278 "a out of range at k=" << k << ", contract " << i);
279 }
280 }
281
282 Distribution dist (nBuckets_, 0.0, maximum_);
283 for (Size i = 0; i < nBuckets_; i++) {
284 dist.addDensity (i, p[i] / dx);
285 dist.addAverage (i, a[i]);
286 }
287
288 return dist;
289 }
290
291 //--------------------------------------------------------------------------
293 //--------------------------------------------------------------------------
294 QL_REQUIRE (loss >= 0, "loss " << loss << " must be >= 0");
295 Real dx = maximum_ / nBuckets_;
296 for (Size i = i0; i < nBuckets_; i++)
297 if (dx * i > loss + epsilon_) return i - 1;
298 return nBuckets_;
299 }
300
301 //--------------------------------------------------------------------------
302 Distribution LossDistMonteCarlo::operator()(const vector<Real>& nominals,
303 const vector<Real>& probabilities) const {
304 //--------------------------------------------------------------------------
305 Distribution dist (nBuckets_, 0.0, maximum_);
306 // KnuthUniformRng rng(seed_);
307 // LecuyerUniformRng rng(seed_);
309 for (Size i = 0; i < simulations_; i++) {
310 Real e = 0;
311 for (Size j = 0; j < nominals.size(); j++) {
312 Real r = rng.next().value;
313 if (r <= probabilities[j])
314 e += nominals[j];
315 }
316 dist.add (e + epsilon_);
317 }
318
319 dist.normalize();
320
321 return dist;
322 }
323
324}
Binomial probability distribution function.
Real operator()(std::vector< Real > p) const
Cumulative binomial distribution function.
void add(Real value)
void addDensity(int bucket, Real value)
void addAverage(int bucket, Real value)
std::vector< Real > excessProbability_
std::vector< Real > probability_
std::vector< Real > probability() const
Distribution operator()(Size n, Real volume, Real probability) const
Distribution operator()(const std::vector< Real > &volumes, const std::vector< Real > &probabilities) const override
int locateTargetBucket(Real loss, Size i0=0) const
std::vector< Real > excessProbability_
Distribution operator()(Real volume, const std::vector< Real > &probabilities) const
static Real binomialProbabilityOfAtLeastNEvents(int n, std::vector< Real > &p)
static std::vector< Real > probabilityOfNEvents(std::vector< Real > &p)
static Real probabilityOfAtLeastNEvents(int n, std::vector< Real > &p)
static Real binomialProbabilityOfNEvents(int n, std::vector< Real > &p)
Distribution operator()(const std::vector< Real > &volumes, const std::vector< Real > &probabilities) const override
Uniform random number generator.
Real operator()(std::vector< Real > p) const
Real operator()(std::vector< Real > p) const
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Loss distributions and probability of n defaults.
Mersenne Twister uniform random number generator.
Definition: any.hpp:35
STL namespace.
ext::shared_ptr< YieldTermStructure > r