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/*! \file distribution.cpp
21 \brief Discretized probability density and cumulative probability
22*/
23
24#include <ql/experimental/credit/distribution.hpp>
25#include <ql/math/comparison.hpp>
26#include <ql/errors.hpp>
27#include <algorithm>
28#include <functional>
29
30namespace QuantLib {
31
32 //-------------------------------------------------------------------------
33 Distribution::Distribution (int nBuckets, Real xmin, Real xmax)
34 //-------------------------------------------------------------------------
35 : size_(nBuckets),
36 xmin_(xmin), xmax_(xmax), count_(nBuckets),
37 x_(nBuckets,0), dx_(nBuckets,0),
38 density_(nBuckets,0),
39 cumulativeDensity_(nBuckets,0),
40 excessProbability_(nBuckets,0),
41 cumulativeExcessProbability_(nBuckets,0),
42 average_(nBuckets,0),
43 overFlow_(0), underFlow_(0),
44 isNormalized_(false) {
45 for (int i = 0; i < nBuckets; i++) {
46 dx_[i] = (xmax - xmin) / nBuckets;
47 x_[i] = (i == 0 ? xmin : x_[i-1] + dx_[i-1]);
48 }
49 // ensure we match exactly the domain, otherwise we might fail the
50 // locate test because of precission mismatches
51 dx_.back() = xmax - x_.back();
52 }
53
54 //-------------------------------------------------------------------------
55 int Distribution::locate (Real x) {
56 //-------------------------------------------------------------------------
57 QL_REQUIRE ((x >= x_.front() || close(x, x_.front())) &&
58 (x <= x_.back() + dx_.back()
59 || close(x, x_.back() + dx_.back())),
60 "coordinate " << x
61 << " out of range [" << x_.front() << "; "
62 << x_.back() + dx_.back() << "]");
63 for (Size i = 0; i < x_.size(); i++) {
64 if (x_[i] > x)
65 return i - 1;
66 }
67 return x_.size() - 1;
68 }
69
70 //-------------------------------------------------------------------------
71 Real Distribution::dx (Real x) {
72 //-------------------------------------------------------------------------
73 int i = locate (x);
74 return dx_[i];
75 }
76
77 //-------------------------------------------------------------------------
78 void Distribution::add (Real value) {
79 //-------------------------------------------------------------------------
80 isNormalized_ = false;
81 if (value < x_.front()) underFlow_++;
82 else {
83 for (Size i = 0; i < count_.size(); i++) {
84 if (x_[i] + dx_[i] > value) {
85 count_[i]++;
86 average_[i] += value;
87 return;
88 }
89 }
90 overFlow_++;
91 }
92 }
93
94 //-------------------------------------------------------------------------
95 void Distribution::addDensity (int bucket, Real value) {
96 //-------------------------------------------------------------------------
97 QL_REQUIRE (bucket >= 0 && bucket < size_, "bucket out of range");
98 isNormalized_ = false;
99 density_[bucket] += value;
100 }
101
102 //-------------------------------------------------------------------------
103 void Distribution::addAverage (int bucket, Real value) {
104 //-------------------------------------------------------------------------
105 QL_REQUIRE (bucket >= 0 && bucket < size_, "bucket out of range");
106 isNormalized_ = false;
107 average_[bucket] += value;
108 }
109
110 //-------------------------------------------------------------------------
111 void Distribution::normalize () {
112 //-------------------------------------------------------------------------
113 if (isNormalized_)
114 return;
115
116 int count = underFlow_ + overFlow_;
117 for (int i = 0; i < size_; i++)
118 count += count_[i];
119
120 excessProbability_[0] = 1.0;
121 cumulativeExcessProbability_[0] = 0.0;
122 for (int i = 0; i < size_; i++) {
123 if (count > 0) {
124 density_[i] = 1.0 / dx_[i] * count_[i] / count;
125 if (count_[i] > 0)
126 average_[i] /= count_[i];
127 }
128 if (density_[i] == 0.0)
129 average_[i] = x_[i] + dx_[i]/2;
130
131 cumulativeDensity_[i] = density_[i] * dx_[i];
132 if (i > 0) {
133 cumulativeDensity_[i] += cumulativeDensity_[i-1];
134 excessProbability_[i] = 1.0 - cumulativeDensity_[i-1];
135// excessProbability_[i] = excessProbability_[i-1]
136// - density_[i-1] * dx_[i-1];
137// cumulativeExcessProbability_[i]
138// = (excessProbability_[i-1] +
139// excessProbability_[i]) / 2 * dx_[i-1]
140// + cumulativeExcessProbability_[i-1];
141 cumulativeExcessProbability_[i]
142 = excessProbability_[i-1] * dx_[i-1]
143 + cumulativeExcessProbability_[i-1];
144 }
145 }
146
147 isNormalized_ = true;
148 }
149
150 //-------------------------------------------------------------------------
151 Real Distribution::confidenceLevel (Real quantil) {
152 //-------------------------------------------------------------------------
153 normalize();
154 for (int i = 0; i < size_; i++) {
155 if (cumulativeDensity_[i] > quantil)
156 return x_[i] + dx_[i];
157 }
158 return x_.back() + dx_.back();
159 }
160
161 //-------------------------------------------------------------------------
162 Real Distribution::expectedValue () {
163 //-------------------------------------------------------------------------
164 normalize();
165 Real expected = 0;
166 for (int i = 0; i < size_; i++) {
167 Real x = x_[i] + dx_[i]/2;
168 expected += x * dx_[i] * density_[i];
169 }
170 return expected;
171 }
172
173 //-------------------------------------------------------------------------
174 Real Distribution::trancheExpectedValue (Real a, Real d) {
175 //-------------------------------------------------------------------------
176 normalize();
177 Real expected = 0;
178 for (int i = 0; i < size_; i++) {
179 Real x = x_[i] + dx_[i]/2;
180 if (x < a)
181 continue;
182 if (x > d)
183 break;
184 expected += (x - a) * dx_[i] * density_[i];
185 }
186
187 expected += (d - a) * (1.0 - cumulativeDensity (x: d));
188
189 return expected;
190 }
191
192// Real Distribution::cumulativeExcessProbability (Real a, Real b) {
193// //normalize();
194// Real integral = 0.0;
195// for (int i = 0; i < size_; i++) {
196// if (x_[i] >= b) break;
197// if (x_[i] >= a)
198// integral += dx_[i] * excessProbability_[i];
199// }
200// return integral;
201// }
202
203 //-------------------------------------------------------------------------
204 Real Distribution::cumulativeExcessProbability (Real a, Real b) {
205 //-------------------------------------------------------------------------
206 normalize();
207 QL_REQUIRE (b <= xmax_,
208 "end of interval " << b << " out of range ["
209 << xmin_ << ", " << xmax_ << "]");
210 QL_REQUIRE (a >= xmin_,
211 "start of interval " << a << " out of range ["
212 << xmin_ << ", " << xmax_ << "]");
213
214 int i = locate (x: a);
215 int j = locate (x: b);
216 return cumulativeExcessProbability_[j]-cumulativeExcessProbability_[i];
217 }
218
219 //-------------------------------------------------------------------------
220 Real Distribution::cumulativeDensity (Real x) {
221 //-------------------------------------------------------------------------
222 Real tiny = dx_.back() * 1e-3;
223 QL_REQUIRE (x > 0, "x must be positive");
224 normalize();
225 for (int i = 0; i < size_; i++) {
226 if (x_[i] + dx_[i] + tiny >= x)
227 return ((x - x_[i]) * cumulativeDensity_[i]
228 + (x_[i] + dx_[i] - x) * cumulativeDensity_[i-1]) / dx_[i];
229 }
230 QL_FAIL ("x = " << x << " beyond distribution cutoff "
231 << x_.back() + dx_.back());
232 }
233
234 //-------------------------------------------------------------------------
235 // Dangerous to perform calls to members after this; transform and clone?
236 void Distribution::tranche (Real attachmentPoint, Real detachmentPoint) {
237 //-------------------------------------------------------------------------
238 QL_REQUIRE (attachmentPoint < detachmentPoint,
239 "attachment >= detachment point");
240 QL_REQUIRE (x_.back() > attachmentPoint &&
241 x_.back()+dx_.back() >= detachmentPoint,
242 "attachment or detachment too large");
243
244 normalize();
245
246 // shift
247 while (x_[0] < attachmentPoint) {
248 x_.erase(position: x_.begin());
249 dx_.erase(position: dx_.begin());
250 count_.erase(position: count_.begin());
251 density_.erase(position: density_.begin());
252 cumulativeDensity_.erase(position: cumulativeDensity_.begin());
253 excessProbability_.erase(position: excessProbability_.begin());
254 }
255
256 // remove losses over detachment point:
257 auto detachPosit = std::find_if(first: x_.begin(), last: x_.end(), pred: [=](Real x){ return x > detachmentPoint; });
258 if(detachPosit != x_.end())
259 x_.erase(first: detachPosit + 1, last: x_.end());
260
261 size_ = x_.size();
262 cumulativeDensity_.erase(first: cumulativeDensity_.begin() + size_,
263 last: cumulativeDensity_.end());
264 cumulativeDensity_.back() = 1.;
265 count_.erase(first: count_.begin() + size_, last: count_.end());
266 dx_.erase(first: dx_.begin() + size_, last: dx_.end());
267
268 // truncate
269 for (Real& i : x_) {
270 i = std::min(a: std::max(a: i - attachmentPoint, b: 0.), b: detachmentPoint - attachmentPoint);
271 }
272
273 density_.clear();
274 excessProbability_.clear();
275 cumulativeExcessProbability_.clear(); //? reuse?
276 density_.push_back(x: (cumulativeDensity_[0]-0.)/dx_[0]);
277 excessProbability_.push_back(x: 1.);
278 for(Integer i=1; i<size_-1; i++) {
279 excessProbability_.push_back(x: 1.-cumulativeDensity_[i-1]);
280 density_.push_back(x: (cumulativeDensity_[i]-
281 cumulativeDensity_[i-1])/dx_[i]);
282 }
283 excessProbability_.push_back(x: 1.-cumulativeDensity_.back());
284 density_.push_back(x: (1.-cumulativeDensity_.back())/dx_.back());
285 }
286
287 //-------------------------------------------------------------------------
288 Distribution ManipulateDistribution::convolve (const Distribution& d1,
289 const Distribution& d2) {
290 //-------------------------------------------------------------------------
291 // force equal constant bucket sizes
292 QL_REQUIRE (d1.dx_[0] == d2.dx_[0], "bucket sizes differ in d1 and d2");
293 for (Size i = 1; i < d1.size(); i++)
294 QL_REQUIRE (d1.dx_[i] == d1.dx_[i-1], "bucket size varies in d1");
295 for (Size i = 1; i < d2.size(); i++)
296 QL_REQUIRE (d2.dx_[i] == d2.dx_[i-1], "bucket size varies in d2");
297
298 // force offset 0
299 QL_REQUIRE (d1.xmin_ == 0.0 && d2.xmin_ == 0.0,
300 "distributions offset larger than 0");
301
302 Distribution dist(d1.size() + d2.size() - 1,
303 0.0, // assuming both distributions have xmin = 0
304 d1.xmax_ + d2.xmax_);
305
306 for (Size i1 = 0; i1 < d1.size(); i1++) {
307 Real dx = d1.dx_[i1];
308 for (Size i2 = 0; i2 < d2.size(); i2++)
309 dist.density_[i1+i2] = d1.density_[i1] * d2.density_[i2] * dx;
310 }
311
312 // update cumulated and excess
313 dist.excessProbability_[0] = 1.0;
314 for (Size i = 0; i < dist.size(); i++) {
315 dist.cumulativeDensity_[i] = dist.density_[i] * dist.dx_[i];
316 if (i > 0) {
317 dist.cumulativeDensity_[i] += dist.cumulativeDensity_[i-1];
318 dist.excessProbability_[i] = dist.excessProbability_[i-1]
319 - dist.density_[i-1] * dist.dx_[i-1];
320 }
321 }
322
323 return dist;
324 }
325
326
327 //-------------------------------------------------------------------------
328 Real Distribution::expectedShortfall (Real percValue) {
329 //-------------------------------------------------------------------------
330 QL_REQUIRE(percValue >= 0. && percValue <= 1.,
331 "Incorrect percentile");
332 normalize();
333 Real expected = 0;
334 Integer iVal = locate(x: confidenceLevel(quantil: percValue));
335
336 if(iVal == size_-1) return x_.back();
337
338 for (int i = iVal; i < size_; i++)
339 expected += x_[i] *
340 (cumulativeDensity_[i] - cumulativeDensity_[i-1]);
341 return expected/(1.-cumulativeDensity_.at(n: iVal));
342 }
343
344}
345

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