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
20#include <ql/experimental/credit/lossdistribution.hpp>
21#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
22
23using namespace std;
24
25namespace QuantLib {
26
27 //--------------------------------------------------------------------------
28 Real LossDist::binomialProbabilityOfNEvents(int n, vector<Real>& p) {
29 //--------------------------------------------------------------------------
30 BinomialDistribution binomial (p[0], p.size());
31 return binomial(n);
32 }
33
34 //--------------------------------------------------------------------------
35 Real LossDist::binomialProbabilityOfAtLeastNEvents(int n, vector<Real>& p) {
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 //--------------------------------------------------------------------------
112 Real LossDist::probabilityOfAtLeastNEvents (int k, vector<Real>& p) {
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 //--------------------------------------------------------------------------
128 Real ProbabilityOfNEvents::operator()(vector<Real> p) const {
129 //--------------------------------------------------------------------------
130 return LossDist::probabilityOfNEvents (k: n_, p);
131 }
132
133 //--------------------------------------------------------------------------
134 Real ProbabilityOfAtLeastNEvents::operator()(vector<Real> p) const {
135 //--------------------------------------------------------------------------
136 return LossDist::probabilityOfAtLeastNEvents (k: n_, p);
137 }
138
139 //--------------------------------------------------------------------------
140 Real BinomialProbabilityOfAtLeastNEvents::operator()(vector<Real> p) const {
141 //--------------------------------------------------------------------------
142 return LossDist::binomialProbabilityOfAtLeastNEvents(n: n_, p);
143 }
144
145 //--------------------------------------------------------------------------
146 Distribution LossDistBinomial::operator()(Size n, Real volume,
147 Real probability) const {
148 //--------------------------------------------------------------------------
149 n_ = n;
150 probability_.clear();
151 probability_.resize(new_size: n_+1, x: 0.0);
152 Distribution dist (nBuckets_, 0.0, maximum_);
153 BinomialDistribution binomial (probability, n);
154 for (Size i = 0; i <= n; i++) {
155 if (volume_ * i <= maximum_) {
156 probability_[i] = binomial(i);
157 Size bucket = dist.locate(x: volume * i);
158 dist.addDensity (bucket, value: probability_[i] / dist.dx(k: bucket));
159 dist.addAverage (bucket, value: volume * i);
160 }
161 }
162
163 excessProbability_.clear();
164 excessProbability_.resize(new_size: n_+1, x: 0.0);
165 excessProbability_[n_] = probability_[n_];
166 for (int k = n_-1; k >= 0; k--)
167 excessProbability_[k] = excessProbability_[k+1] + probability_[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()(n: nominals.size(), volume: nominals[0], probability: probabilities[0]);
179 }
180
181 //--------------------------------------------------------------------------
182 Distribution LossDistHomogeneous::operator()(Real volume,
183 const vector<Real>& p) const {
184 //--------------------------------------------------------------------------
185 volume_ = volume;
186 n_ = p.size();
187 probability_.clear();
188 probability_.resize(new_size: n_+1, x: 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(new_size: n_+1, x: 0.0);
201 excessProbability_[n_] = probability_[n_];
202 for (int k = n_ - 1; k >= 0; k--)
203 excessProbability_[k] = excessProbability_[k+1] + probability_[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(x: volume * i);
209 dist.addDensity (bucket, value: probability_[i] / dist.dx(k: bucket));
210 dist.addAverage (bucket, value: 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()(volume: nominals[0], p: 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 (loss: a[k] + L, i0: 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 (bucket: i, value: p[i] / dx);
285 dist.addAverage (bucket: i, value: a[i]);
286 }
287
288 return dist;
289 }
290
291 //--------------------------------------------------------------------------
292 int LossDistBucketing::locateTargetBucket (Real loss, Size i0) const {
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_);
308 MersenneTwisterUniformRng 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 (value: e + epsilon_);
317 }
318
319 dist.normalize();
320
321 return dist;
322 }
323
324}
325

source code of quantlib/ql/experimental/credit/lossdistribution.cpp