1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Roland Lichters
5 Copyright (C) 2009, 2014 Jose Aparicio
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#ifndef quantlib_randomdefault_latent_model_hpp
22#define quantlib_randomdefault_latent_model_hpp
23
24#include <ql/experimental/credit/basket.hpp>
25#include <ql/experimental/credit/constantlosslatentmodel.hpp>
26#include <ql/experimental/credit/defaultlossmodel.hpp>
27#include <ql/experimental/math/gaussiancopulapolicy.hpp>
28#include <ql/experimental/math/latentmodel.hpp>
29#include <ql/experimental/math/tcopulapolicy.hpp>
30#include <ql/math/beta.hpp>
31#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
32#include <ql/math/randomnumbers/sobolrsg.hpp>
33#include <ql/math/solvers1d/brent.hpp>
34#include <ql/math/statistics/histogram.hpp>
35#include <ql/math/statistics/riskstatistics.hpp>
36#include <ql/tuple.hpp>
37#include <utility>
38
39/* Intended to replace
40 ql\experimental\credit\randomdefaultmodel.Xpp
41*/
42
43namespace QuantLib {
44
45 /*! Simulation event trait class template forward declaration.
46 Each latent model will be modelling different entities according to the
47 meaning of the model function which depends on the random \$ Y_i\$
48 variables. Still the generation of the factors and variables it is common to
49 any model. Only within a given transformation function the model and event
50 generated becomes concrete.
51
52 However here these simulations are already made specific to a default event.
53 Yet other variables contingent to default can be modelled (recovery,
54 market...) So we are placed in a less generic stage where default is
55 modelled possibly jointly with other unespecified magnitudes.
56
57 Another role of this trait class is to compact in memory the simulation
58 data. The statistic post processing needs to have the results stored in
59 memory and simulations can not be consumed at generation time, typically
60 because some statistics are conditional on others (e.g. ESF) or/and
61 parametric (percentile, etc...)
62
63 Simulation events do not derive from each other, and they are specialized
64 for each type; duck typing applies for variable names (see the statistic
65 methods)
66 */
67 // replaces class Loss
68 template <class simEventOwner> struct simEvent;
69
70
71 /*! Base class for latent model monte carlo simulation. Independent of the
72 copula type and the generator.
73 Generates the factors and variable samples and determines event threshold
74 but it is not responsible for actual event specification; thats the derived
75 classes responsibility according to what they model.
76 Derived classes need mainly to implement nextSample (Worker::nextSample in
77 the multithreaded version) to compute the simulation event generated, if
78 any, from the latent variables sample. They also have the accompanying
79 event trait to specify.
80 */
81 /* CRTP used for performance to avoid virtual table resolution in the Monte
82 Carlo. Not only in sample generation but access; quite an amount of time can
83 go into statistics computation, for a portfolio of tens of thousands
84 positions that part of the problem will be starting to overtake the
85 simulation costs.
86
87 \todo: someone with sound experience on cache misses look into this, the
88 statistics will be getting memory in and out of the cpu heavily and it
89 might be possible to get performance out of that.
90 \todo: parallelize the statistics computation, things like Var/ESF splits
91 are very expensive.
92 \todo: consider another design, taking the statistics outside the models.
93 */
94 template<template <class, class> class derivedRandomLM, class copulaPolicy,
95 class USNG = SobolRsg>
96 class RandomLM : public virtual LazyObject,
97 public virtual DefaultLossModel {
98 private:
99 // Takes the parents type, all children have the same type, the
100 // random generation is performed in this class only.
101 typedef typename LatentModel<copulaPolicy>::template FactorSampler<USNG>
102 copulaRNG_type;
103 protected:
104 RandomLM(Size numFactors, Size numLMVars, copulaPolicy copula, Size nSims, BigNatural seed)
105 : seed_(seed), numFactors_(numFactors), numLMVars_(numLMVars), nSims_(nSims),
106 copula_(std::move(copula)) {}
107
108 void update() override {
109 simsBuffer_.clear();
110 // tell basket to notify instruments, etc, we are invalid
111 if (!basket_.empty())
112 basket_->notifyObservers();
113 LazyObject::update();
114 }
115
116 void performCalculations() const override {
117 static_cast<const derivedRandomLM<copulaPolicy, USNG>* >(
118 this)->initDates();//in update?
119 copulasRng_ = ext::make_shared<copulaRNG_type>(copula_, seed_);
120 performSimulations();
121 }
122
123 void performSimulations() const {
124 // Next sequence should determine the event and push it into buffer
125 for (Size i = nSims_; i != 0U; i--) {
126 const std::vector<Real>& sample =
127 copulasRng_->nextSequence().value;
128 static_cast<const derivedRandomLM<copulaPolicy, USNG>* >(
129 this)->nextSample(sample);
130 // alternatively make call an explicit local method...
131 }
132 }
133
134 /* Method to access simulation results and avoiding a copy of
135 each thread results buffer. PerformCalculations should have been called.
136 Here in the monothread version this method is redundant/trivial but
137 serves to detach the statistics access to the way the simulations are
138 stored.
139 */
140 const std::vector<simEvent<derivedRandomLM<copulaPolicy, USNG> > >&
141 getSim(const Size iSim) const { return simsBuffer_[iSim]; }
142
143 /* Allows statistics to be written generically for fixed and random
144 recovery rates. */
145 Real getEventRecovery(
146 const simEvent<derivedRandomLM<copulaPolicy, USNG> >& evt) const
147 {
148 return static_cast<const derivedRandomLM<copulaPolicy, USNG>* >(
149 this)->getEventRecovery(evt);
150 }
151
152 //! \name Statistics, DefaultLossModel interface.
153 // These are virtual and allow for children-specific optimization and
154 // variance reduction. The virtual table is ok, they are not part
155 // of the simulation.
156 //@{
157 /*! Returns the probaility of having a given or larger number of
158 defaults in the basket portfolio at a given time.
159 */
160 Probability probAtLeastNEvents(Size n, const Date& d) const override;
161 /*! Order of results refers to the simulated (super)pool not the
162 basket's pool.
163 Notice that this statistic suffers from heavy dispersion. To see
164 techniques to improve it (not implemented here) see:
165 Joshi, M., D. Kainth. 2004. Rapid and accurate development of prices
166 and Greeks for nth to default credit swaps in the Li model. Quantitative
167 Finance, Vol. 4. Institute of Physics Publishing, London, UK, 266-275
168 and:
169 Chen, Z., Glasserman, P. 'Fast pricing of basket default swaps' in
170 Operations Research Vol. 56, No. 2, March/April 2008, pp. 286-303
171 */
172 std::vector<Probability> probsBeingNthEvent(Size n, const Date& d) const override;
173 //! Pearsons' default probability correlation.
174 Real defaultCorrelation(const Date& d, Size iName, Size jName) const override;
175 Real expectedTrancheLoss(const Date& d) const override;
176 virtual std::pair<Real, Real> expectedTrancheLossInterval(const Date& d,
177 Probability confidencePerc) const;
178 std::map<Real, Probability> lossDistribution(const Date& d) const override;
179 virtual Histogram computeHistogram(const Date& d) const;
180 Real expectedShortfall(const Date& d, Real percent) const override;
181 Real percentile(const Date& d, Real percentile) const override;
182 /*! Returns the VaR value for a given percentile and the 95 confidence
183 interval of that value. */
184 virtual ext::tuple<Real, Real, Real> percentileAndInterval(
185 const Date& d, Real percentile) const;
186 /*! Distributes the total VaR amount along the portfolio counterparties.
187 The passed loss amount is in loss units.
188 */
189 std::vector<Real> splitVaRLevel(const Date& date, Real loss) const override;
190 /*! Distributes the total VaR amount along the portfolio
191 counterparties.
192
193 Provides confidence interval for split so that portfolio
194 optimization can be performed outside those limits.
195
196 The passed loss amount is in loss units.
197 */
198 virtual std::vector<std::vector<Real> > splitVaRAndError(
199 const Date& date, Real loss, Probability confInterval) const;
200 //@}
201 public:
202 ~RandomLM() override = default;
203
204 private:
205 BigNatural seed_;
206 protected:
207 const Size numFactors_;
208 const Size numLMVars_;
209
210 const Size nSims_;
211
212 mutable std::vector<std::vector<simEvent<derivedRandomLM<copulaPolicy,
213 USNG > > > > simsBuffer_;
214
215 mutable copulaPolicy copula_;
216 mutable ext::shared_ptr<copulaRNG_type> copulasRng_;
217
218 // Maximum time inversion horizon
219 static const Size maxHorizon_ = 4050; // over 11 years
220 // Inversion probability limits are computed by children in initdates()
221 };
222
223
224 /* ---- Statistics --------------------------------------------------- */
225
226 template<template <class, class> class D, class C, class URNG>
227 Probability RandomLM<D, C, URNG>::probAtLeastNEvents(Size n,
228 const Date& d) const
229 {
230 calculate();
231 Date today = Settings::instance().evaluationDate();
232
233 QL_REQUIRE(d>today, "Date for statistic must be in the future.");
234 // casted to natural to avoid warning, we have just checked the sign
235 Natural val = d.serialNumber() - today.serialNumber();
236
237 if(n==0) return 1.;
238
239 Real counts = 0.;
240 for(Size iSim=0; iSim < nSims_; iSim++) {
241 Size simCount = 0;
242 const std::vector<simEvent<D<C, URNG> > >& events =
243 getSim(iSim);
244 for(Size iEvt=0; iEvt < events.size(); iEvt++)
245 // duck type on the members:
246 if(val > events[iEvt].dayFromRef) simCount++;
247 if(simCount >= n) counts++;
248 }
249 return counts/nSims_;
250 // \todo Provide confidence interval
251 }
252
253 template<template <class, class> class D, class C, class URNG>
254 std::vector<Probability> RandomLM<D, C, URNG>::probsBeingNthEvent(Size n,
255 const Date& d) const
256 {
257 calculate();
258 Size basketSize = basket_->size();
259
260 QL_REQUIRE(n>0 && n<=basketSize, "Impossible number of defaults.");
261 Date today = Settings::instance().evaluationDate();
262
263 QL_REQUIRE(d>today, "Date for statistic must be in the future.");
264 // casted to natural to avoid warning, we have just checked the sign
265 Natural val = d.serialNumber() - today.serialNumber();
266
267 std::vector<Probability> hitsByDate(basketSize, 0.);
268 for(Size iSim=0; iSim < nSims_; iSim++) {
269 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
270 std::map<unsigned short, unsigned short> namesDefaulting;
271 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
272 // if event is within time horizon...
273 if(val > events[iEvt].dayFromRef)
274 //...count it. notice insertion sorts by date.
275 namesDefaulting.insert(std::make_pair<unsigned short,
276 unsigned short>(events[iEvt].dayFromRef,
277 events[iEvt].nameIdx));
278 }
279 if(namesDefaulting.size() >= n) {
280 std::map<unsigned short, unsigned short>::const_iterator
281 itdefs = namesDefaulting.begin();
282 // locate nth default in time:
283 std::advance(i&: itdefs, n: n-1);
284 // update statistic:
285 hitsByDate[itdefs->second]++;
286 }
287 }
288 std::transform(hitsByDate.begin(), hitsByDate.end(),
289 hitsByDate.begin(),
290 [this](Real x){ return x/nSims_; });
291 return hitsByDate;
292 // \todo Provide confidence interval
293 }
294
295
296 template<template <class, class> class D, class C, class URNG>
297 Real RandomLM<D, C, URNG>::defaultCorrelation(const Date& d,
298 Size iName, Size jName) const
299 {
300 // a control variate with the probabilities is possible
301 calculate();
302 Date today = Settings::instance().evaluationDate();
303
304 QL_REQUIRE(d>today, "Date for statistic must be in the future.");
305 // casted to natural to avoid warning, we have just checked the sign
306 Natural val = d.serialNumber() - today.serialNumber();
307
308 Real expectedDefiDefj = 0.;// E[1_i 1_j]
309 // the rest of magnitudes have known values (probabilities) but that
310 // would distort the simulation results.
311 Real expectedDefi = 0.;
312 Real expectedDefj = 0.;
313 for(Size iSim=0; iSim < nSims_; iSim++) {
314 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
315 Real imatch = 0., jmatch = 0.;
316 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
317 if((val > events[iEvt].dayFromRef) &&
318 (events[iEvt].nameIdx == iName)) imatch = 1.;
319 if((val > events[iEvt].dayFromRef) &&
320 (events[iEvt].nameIdx == jName)) jmatch = 1.;
321 }
322 expectedDefiDefj += imatch * jmatch;
323 expectedDefi += imatch;
324 expectedDefj += jmatch;
325 }
326 expectedDefiDefj = expectedDefiDefj / (nSims_-1);// unbiased
327 expectedDefi = expectedDefi / nSims_;
328 expectedDefj = expectedDefj / nSims_;
329
330 return (expectedDefiDefj - expectedDefi*expectedDefj) /
331 std::sqrt(x: (expectedDefi*expectedDefj*(1.-expectedDefi)
332 *(1.-expectedDefj)));
333 // \todo Provide confidence interval
334 }
335
336
337 template<template <class, class> class D, class C, class URNG>
338 Real RandomLM<D, C, URNG>::expectedTrancheLoss(
339 const Date& d) const {
340 return expectedTrancheLossInterval(d, confidencePerc: 0.95).first;
341 }
342
343
344 template<template <class, class> class D, class C, class URNG>
345 std::pair<Real, Real> RandomLM<D, C, URNG>::expectedTrancheLossInterval(
346 const Date& d, Probability confidencePerc) const
347 {
348 calculate();
349 Date today = Settings::instance().evaluationDate();
350 Date::serial_type val = d.serialNumber() - today.serialNumber();
351
352 Real attachAmount = basket_->attachmentAmount();
353 Real detachAmount = basket_->detachmentAmount();
354
355 // Real trancheLoss= 0.;
356 GeneralStatistics lossStats;
357 for(Size iSim=0; iSim < nSims_; iSim++) {
358 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
359
360 Real portfSimLoss=0.;
361 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
362 // if event is within time horizon...
363 if(val > static_cast<Date::serial_type>(
364 events[iEvt].dayFromRef)) {
365 Size iName = events[iEvt].nameIdx;
366 // ...and is contained in the basket.
367 portfSimLoss +=
368 basket_->exposure(name: basket_->names()[iName],
369 Date(events[iEvt].dayFromRef +
370 today.serialNumber())) *
371 (1.-getEventRecovery(evt: events[iEvt]));
372 }
373 }
374 lossStats.add(// d ates? current losses? realized defaults, not yet
375 value: std::min(a: std::max(a: portfSimLoss - attachAmount, b: 0.),
376 b: detachAmount - attachAmount) );
377 }
378 return std::make_pair(x: lossStats.mean(), y: lossStats.errorEstimate() *
379 InverseCumulativeNormal::standard_value(x: 0.5*(1.+confidencePerc)));
380 }
381
382
383 template<template <class, class> class D, class C, class URNG>
384 std::map<Real, Probability> RandomLM<D, C, URNG>::lossDistribution(const Date& d) const {
385
386 Histogram hist = computeHistogram(d);
387 std::map<Real, Probability> distrib;
388
389 // prob of losses less or equal to
390 Real suma = hist.frequency(i: 0);
391 distrib.insert(x: std::make_pair(x: 0., y&: suma));
392 for(Size i=1; i<hist.bins(); i++) {
393 suma += hist.frequency(i);
394 distrib.insert(x: std::make_pair( x: hist.breaks()[i-1], y&: suma ));
395 }
396 return distrib;
397 }
398
399
400 template<template <class, class> class D, class C, class URNG>
401 Histogram RandomLM<D, C, URNG>::computeHistogram(const Date& d) const {
402 std::vector<Real> data;
403 std::set<Real> keys;// attainable loss values
404 keys.insert(x: 0.);
405 Date today = Settings::instance().evaluationDate();
406 Date::serial_type val = d.serialNumber() - today.serialNumber();
407 // redundant test? should have been tested by the basket caller?
408 QL_REQUIRE(d >= today,
409 "Requested percentile date must lie after computation date.");
410 calculate();
411
412 Real attachAmount = basket_->attachmentAmount();
413 Real detachAmount = basket_->detachmentAmount();
414
415 for(Size iSim=0; iSim < nSims_; iSim++) {
416 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
417
418 Real portfSimLoss=0.;
419 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
420 if(val > static_cast<Date::serial_type>(
421 events[iEvt].dayFromRef)) {
422 Size iName = events[iEvt].nameIdx;
423 // test needed (here and the others) to reuse simulations:
424 // if(basket_->pool()->has(copula_->pool()->names()[iName]))
425 portfSimLoss +=
426 basket_->exposure(name: basket_->names()[iName],
427 Date(events[iEvt].dayFromRef +
428 today.serialNumber())) *
429 (1.-getEventRecovery(evt: events[iEvt]));
430 }
431 }
432 data.push_back(x: std::min(a: std::max(a: portfSimLoss - attachAmount, b: 0.),
433 b: detachAmount - attachAmount));
434 keys.insert(x: data.back());
435 }
436 // avoid using as many points as in the simulation.
437 Size nPts = std::min<Size>(a: data.size(), b: 150);// fix
438 return Histogram(data.begin(), data.end(), nPts);
439 }
440
441
442 template<template <class, class> class D, class C, class URNG>
443 Real RandomLM<D, C, URNG>::expectedShortfall(const Date& d,
444 Real percent) const {
445
446 const Date today = Settings::instance().evaluationDate();
447 QL_REQUIRE(d >= today,
448 "Requested percentile date must lie after computation date.");
449 calculate();
450
451 Real attachAmount = basket_->attachmentAmount();
452 Real detachAmount = basket_->detachmentAmount();
453
454 Date::serial_type val = d.serialNumber() - today.serialNumber();
455 if(val <= 0) return 0.;// plus basket realized losses
456
457 //GenericRiskStatistics<GeneralStatistics> statsX;
458 std::vector<Real> losses;
459 for(Size iSim=0; iSim < nSims_; iSim++) {
460 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
461 Real portfSimLoss=0.;
462 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
463 if(val > static_cast<Date::serial_type>(
464 events[iEvt].dayFromRef)) {
465 Size iName = events[iEvt].nameIdx;
466 // ...and is contained in the basket.
467 //if(basket_->pool()->has(copula_->pool()->names()[iName]))
468 portfSimLoss +=
469 basket_->exposure(name: basket_->names()[iName],
470 Date(events[iEvt].dayFromRef +
471 today.serialNumber())) *
472 (1.-getEventRecovery(evt: events[iEvt]));
473 }
474 }
475 portfSimLoss = std::min(a: std::max(a: portfSimLoss - attachAmount, b: 0.),
476 b: detachAmount - attachAmount);
477 losses.push_back(x: portfSimLoss);
478 }
479
480 std::sort(first: losses.begin(), last: losses.end());
481 Real posit = std::ceil(x: percent * nSims_);
482 posit = posit >= 0. ? posit : 0.;
483 Size position = static_cast<Size>(posit);
484 Real perctlInf = losses[position];//q_{\alpha}
485
486 // the prob of values strictly larger than the quantile value.
487 Probability probOverQ =
488 static_cast<Real>(std::distance(first: losses.begin() + position,
489 last: losses.end())) / static_cast<Real>(nSims_);
490
491 return ( perctlInf * (1.-percent-probOverQ) +//<-correction term
492 std::accumulate(first: losses.begin() + position, last: losses.end(),
493 init: Real(0.))/nSims_
494 )/(1.-percent);
495
496 /* Alternative ESF definition; find the first loss larger than the
497 one of the percentile. Notice the choice here, the expected shortfall
498 is understood in the sense that we are looking for the average given
499 than losses are above a certain value rather than above a certain
500 probability:
501 (Unlikely to be the algorithm of choice)*/
502 /*
503 std::vector<Real>::iterator itPastPerc =
504 std::find_if(losses.begin() + position, losses.end(),
505 [=](Real x){ return x >= perctlInf; });
506 // notice if the sample is flat at the end this might be zero
507 Size pointsOverVal = nSims_ - std::distance(itPastPerc, losses.end());
508 return pointsOverVal == 0 ? 0. :
509 std::accumulate(itPastPerc, losses.end(), 0.) / pointsOverVal;
510 */
511
512 /* For the definition of ESF see for instance: 'Quantitative Risk
513 Management' by A.J. McNeil, R.Frey and P.Embrechts, princeton series in
514 finance, 2005; equations on page 39 sect 2.12:
515 $q_{\alpha}(F) = inf{x \in R : F(x) \le \alpha}$
516 and equation 2.25 on p. 45:
517 $ESF_{\alpha} = \frac{1}{1-\alpha} [E(L; L \ge q_{\alpha} ) +
518 q_{\alpha} (1-\alpha-P(L \ge q_{\alpha})) ]$
519 The second term accounts for non continuous distributions.
520 */
521 }
522
523
524 template<template <class, class> class D, class C, class URNG>
525 Real RandomLM<D, C, URNG>::percentile(const Date& d, Real perc) const {
526 // need to specify return type in tuples' get is parametric
527 return ext::get<0>(percentileAndInterval(d, percentile: perc));
528 }
529
530
531 /* See Appendix-A of "Evaluating value-at-risk methodologies: Accuracy
532 versus computational time.", M. Pritsker, Wharton FIC, November 1996
533 Strictly speaking this gives the interval with a 95% probability of
534 the true value being within the interval; which is different to the error
535 of the stimator just computed. See the reference for a discussion.
536 */
537 template<template <class, class> class D, class C, class URNG>
538 ext::tuple<Real, Real, Real> RandomLM<D, C, URNG>::percentileAndInterval(const Date& d,
539 Real percentile) const {
540
541 QL_REQUIRE(percentile >= 0. && percentile <= 1.,
542 "Incorrect percentile");
543 calculate();
544
545 Real attachAmount = basket_->attachmentAmount();
546 Real detachAmount = basket_->detachmentAmount();
547
548 std::vector<Real> rankLosses;
549 Date today = Settings::instance().evaluationDate();
550 Date::serial_type val = d.serialNumber() - today.serialNumber();
551 for(Size iSim=0; iSim < nSims_; iSim++) {
552 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
553 Real portfSimLoss=0.;
554 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
555 if(val > static_cast<Date::serial_type>(
556 events[iEvt].dayFromRef)) {
557 Size iName = events[iEvt].nameIdx;
558 // if(basket_->pool()->has(copula_->pool()->names()[iName]))
559 portfSimLoss +=
560 basket_->exposure(name: basket_->names()[iName],
561 Date(events[iEvt].dayFromRef +
562 today.serialNumber())) *
563 (1.-getEventRecovery(evt: events[iEvt]));
564 }
565 }
566 portfSimLoss = std::min(a: std::max(a: portfSimLoss - attachAmount, b: 0.),
567 b: detachAmount - attachAmount);
568 // update dataset for rank stat:
569 rankLosses.push_back(x: portfSimLoss);
570 }
571
572 std::sort(first: rankLosses.begin(), last: rankLosses.end());
573 Size quantilePosition = static_cast<Size>(floor(x: nSims_*percentile));
574 Real quantileValue = rankLosses[quantilePosition];
575
576 // compute confidence interval:
577 const Probability confInterval = 0.95;// as an argument?
578 Real lowerPercentile, upperPercentile;
579 Size r = quantilePosition - 1;
580 Size s = quantilePosition + 1;
581 bool rLocked = false,
582 sLocked = false;
583 // Size rfinal = 0,
584 // sfinal = 0;
585 for(Size delta=1; delta < quantilePosition; delta++) {
586 Real cached =
587 incompleteBetaFunction(a: Real(s), b: Real(nSims_+1-s),
588 x: percentile, accuracy: 1.e-8, maxIteration: 500);
589 Real pMinus =
590 /* There was a fix in the repository on the gammadistribution. It
591 might impact these, it might be neccesary to multiply these values
592 by '-1'*/
593 incompleteBetaFunction(a: Real(r+1), b: Real(nSims_-r),
594 x: percentile, accuracy: 1.e-8, maxIteration: 500)
595 - cached;
596 Real pPlus =
597 incompleteBetaFunction(a: Real(r), b: Real(nSims_-r+1),
598 x: percentile, accuracy: 1.e-8, maxIteration: 500)
599 - cached;
600 if((pMinus > confInterval) && !rLocked ) {
601 // rfinal = r + 1;
602 rLocked = true;
603 }
604 if((pPlus >= confInterval) && !sLocked) {
605 // sfinal = s;
606 sLocked = true;
607 }
608 if(rLocked && sLocked) break;
609 r--;
610 s++;
611 s = std::min(a: nSims_-1, b: s);
612 }
613 lowerPercentile = rankLosses[r];
614 upperPercentile = rankLosses[s];
615
616 return ext::make_tuple(args&: quantileValue, args&: lowerPercentile, args&: upperPercentile);
617 }
618
619
620 template<template <class, class> class D, class C, class URNG>
621 std::vector<Real> RandomLM<D, C, URNG>::splitVaRLevel(
622 const Date& date, Real loss) const
623 {
624 std::vector<Real> varLevels = splitVaRAndError(date, loss, confInterval: 0.95)[0];
625 // turn relative units into absolute:
626 std::transform(varLevels.begin(), varLevels.end(), varLevels.begin(),
627 [=](Real x) -> Real { return x * loss; });
628 return varLevels;
629 }
630
631
632 // parallelize this one(if possible), it is really expensive
633 template<template <class, class> class D, class C, class URNG>
634 /* FIX ME: some trouble on limit cases, like zero loss or no losses over the
635 requested level.*/
636 std::vector<std::vector<Real> > RandomLM<D, C, URNG>::splitVaRAndError(const Date& date, Real loss,
637 Probability confInterval) const
638 {
639 /* Check 'loss' value integrity: i.e. is within tranche limits? (should
640 have been done basket...)*/
641 calculate();
642
643 Real attachAmount = basket_->attachmentAmount();
644 Real detachAmount = basket_->detachmentAmount();
645 Size numLiveNames = basket_->remainingSize();
646
647 std::vector<Real> split(numLiveNames, 0.);
648 std::vector<GeneralStatistics> splitStats(numLiveNames,
649 GeneralStatistics());
650 Date today = Settings::instance().evaluationDate();
651 Date::serial_type val = date.serialNumber() - today.serialNumber();
652
653 for(Size iSim=0; iSim < nSims_; iSim++) {
654 const std::vector<simEvent<D<C, URNG> > >& events = getSim(iSim);
655 Real portfSimLoss=0.;
656 //std::vector<Real> splitBuffer(numLiveNames_, 0.);
657 std::vector<simEvent<D<C, URNG> > > splitEventsBuffer;
658
659 for(Size iEvt=0; iEvt < events.size(); iEvt++) {
660 if(val > static_cast<Date::serial_type>(
661 events[iEvt].dayFromRef)) {
662 Size iName = events[iEvt].nameIdx;
663 // if(basket_->pool()->has(copula_->pool()->names()[iName])) {
664 portfSimLoss +=
665 basket_->exposure(name: basket_->names()[iName],
666 Date(events[iEvt].dayFromRef +
667 today.serialNumber())) *
668 (1.-getEventRecovery(evt: events[iEvt]));
669 //and will sort later if buffer applies:
670 splitEventsBuffer.push_back(events[iEvt]);
671 }
672 }
673 portfSimLoss = std::min(a: std::max(a: portfSimLoss - attachAmount, b: 0.),
674 b: detachAmount - attachAmount);
675
676 /* second pass; split is conditional to total losses within target
677 losses/percentile: */
678 Real ptflCumulLoss = 0.;
679 if(portfSimLoss > loss) {
680 std::sort(splitEventsBuffer.begin(), splitEventsBuffer.end());
681 //NOW THIS:
682 split.assign(n: numLiveNames, val: 0.);
683 /* if the name triggered a loss in the portf limits assign
684 this loss to that name.. */
685 for(Size i=0; i<splitEventsBuffer.size(); i++) {
686 Size iName = splitEventsBuffer[i].nameIdx;
687 Real lossName =
688 // allows amortizing (others should be like this)
689 // basket_->remainingNotionals(Date(simsBuffer_[i].dayFromRef +
690 // today.serialNumber()))[iName] *
691 basket_->exposure(name: basket_->names()[iName],
692 Date(splitEventsBuffer[i].dayFromRef +
693 today.serialNumber())) *
694 (1.-getEventRecovery(evt: splitEventsBuffer[i]));
695
696 Real tranchedLossBefore =
697 std::min(a: std::max(a: ptflCumulLoss - attachAmount, b: 0.),
698 b: detachAmount - attachAmount);
699 ptflCumulLoss += lossName;
700 Real tranchedLossAfter =
701 std::min(a: std::max(a: ptflCumulLoss - attachAmount, b: 0.),
702 b: detachAmount - attachAmount);
703 // assign new losses:
704 split[iName] += tranchedLossAfter - tranchedLossBefore;
705 }
706 for(Size iName=0; iName<numLiveNames; iName++) {
707 splitStats[iName].add(value: split[iName] /
708 std::min(a: std::max(a: ptflCumulLoss - attachAmount, b: 0.),
709 b: detachAmount - attachAmount) );
710 }
711 }
712 }
713
714 // Compute error in VaR split
715 std::vector<Real> means, rangeUp, rangeDown;
716 Real confidFactor = InverseCumulativeNormal()(0.5+confInterval/2.);
717 for(Size iName=0; iName<numLiveNames; iName++) {
718 means.push_back(x: splitStats[iName].mean());
719 Real error = confidFactor * splitStats[iName].errorEstimate() ;
720 rangeDown.push_back(x: means.back() - error);
721 rangeUp.push_back(x: means.back() + error);
722 }
723
724 std::vector<std::vector<Real> > results;
725 results.push_back(x: means);
726 results.push_back(x: rangeDown);
727 results.push_back(x: rangeUp);
728
729 return results;
730 }
731
732
733
734
735 // --------- Time inversion solver target function: -----------------------
736
737 /* It could be argued that this concept is part of the copula (more generic)
738 In general when the modelled magnitude is parametric one can solve for
739 inversion to get the parameter value for a given magnitude value (provided
740 the modelled variable dependence in invertible). In this particular problem
741 the parameter is Time and it is solved here where we are alredy in the
742 context of default
743 See default transition models for another instance of this inversion.
744 Alternatively use the faster trick (flat HR) mentioned in the code or make
745 the algorithm parametric on the type of interpolation in the default TS.
746 */
747 namespace detail {// not template dependent .....move it
748 //! Utility for the numerical time solver
749 class Root {
750 public:
751 /* See a faster algorithm (neeeds to locate the points) in
752 D.O'KANE p.249 sect 13.5 */
753 Root(const Handle<DefaultProbabilityTermStructure>& dts, Real pd)
754 : dts_(dts), pd_(pd), curveRef_(dts->referenceDate()) {}
755 /* The cast I am forcing here comes from the requirement of 1D
756 solvers to take in a target (cost) function of Real domain. It could
757 be possible to change the template arg F in the 1D solvers to a
758 boost function and then use the (template arg) domain argument type
759 of the function for use with the 'guess' and operator() ?
760 */
761 Real operator()(Real t) const {
762 QL_REQUIRE (t >= 0.0, "t < 0");
763 /* As long as this doesnt involve modifying a mutable member
764 it should be thread safe (they are const methods and access is
765 read only)
766 */
767 return dts_->defaultProbability(d: curveRef_ +
768 Period(static_cast<Integer>(t), Days), extrapolate: true) - pd_;
769 }
770 private:
771 const Handle<DefaultProbabilityTermStructure> dts_;
772 Real pd_;
773 const Date curveRef_;
774 };
775 }
776
777 /*
778 ---------------------------------------------------------------------------
779 ---------------------------------------------------------------------------
780 */
781
782 // move this one to a separte file?
783 /*! Random default with deterministic recovery event type.\par
784 Stores sims results in a bitfield buffer for lean memory storage.
785 Although strictly speaking this is not guaranteed by the compiler it
786 amounts to reducing the memory storage by half.
787 Some computations, like conditional statistics, precise that all sims
788 results be available.
789 */
790 template<class , class > class RandomDefaultLM;
791 template<class copulaPolicy, class USNG>
792 struct simEvent<RandomDefaultLM<copulaPolicy, USNG> > {
793 simEvent(unsigned int n, unsigned int d)
794 : nameIdx(n), dayFromRef(d){}
795 unsigned int nameIdx : 16; // can index up to 65535 names
796 unsigned int dayFromRef : 16; //indexes up to 65535 days ~179 years
797 bool operator<(const simEvent& evt) const {
798 return dayFromRef < evt.dayFromRef;
799 }
800 };
801
802 /*! Default only latent model simulation with trivially fixed recovery
803 amounts.
804 */
805 template<class copulaPolicy, class USNG = SobolRsg>
806 class RandomDefaultLM : public RandomLM<RandomDefaultLM, copulaPolicy, USNG>
807 {
808 private:
809 typedef simEvent<RandomDefaultLM> defaultSimEvent;
810
811 // \todo Consider this to be only a ConstantLossLM instead
812 const ext::shared_ptr<DefaultLatentModel<copulaPolicy> > model_;
813 const std::vector<Real> recoveries_;
814 // for time inversion:
815 Real accuracy_;
816 public:
817 // \todo: Allow a constructor building its own default latent model.
818 explicit RandomDefaultLM(const ext::shared_ptr<DefaultLatentModel<copulaPolicy> >& model,
819 const std::vector<Real>& recoveries = std::vector<Real>(),
820 Size nSims = 0, // stats will crash on div by zero, FIX ME.
821 Real accuracy = 1.e-6,
822 BigNatural seed = 2863311530UL)
823 : RandomLM< ::QuantLib::RandomDefaultLM, copulaPolicy, USNG>(
824 model->numFactors(), model->size(), model->copula(), nSims, seed),
825 model_(model),
826 recoveries_(recoveries.empty() ? std::vector<Real>(model->size(), 0.) : recoveries),
827 accuracy_(accuracy) {
828 // redundant through basket?
829 this->registerWith(Settings::instance().evaluationDate());
830 this->registerWith(model_);
831 }
832 explicit RandomDefaultLM(
833 const ext::shared_ptr<ConstantLossLatentmodel<copulaPolicy> >& model,
834 Size nSims = 0,// stats will crash on div by zero, FIX ME.
835 Real accuracy = 1.e-6,
836 BigNatural seed = 2863311530UL)
837 : RandomLM< ::QuantLib::RandomDefaultLM, copulaPolicy, USNG>
838 (model->numFactors(), model->size(), model->copula(),
839 nSims, seed ),
840 model_(model),
841 recoveries_(model->recoveries()),
842 accuracy_(accuracy)
843 {
844 // redundant through basket?
845 this->registerWith(Settings::instance().evaluationDate());
846 this->registerWith(model_);
847 }
848
849 // grant access to static polymorphism:
850 /* While this works on g++, VC9 refuses to compile it.
851 Not completely sure whos right; individually making friends of the
852 calling members or writting explicitly the derived class T parameters
853 throws the same errors.
854 The access is then open to the member fucntions.
855 Another solution is to use this http://accu.org/index.php/journals/296
856
857 It might well be that gcc is allowing some c11 features silently, which
858 wont pass on a lower gcc version.
859 */
860 friend class RandomLM< ::QuantLib::RandomDefaultLM, copulaPolicy, USNG>;
861 protected:
862 void nextSample(const std::vector<Real>& values) const;
863 void initDates() const {
864 /* Precalculate horizon time default probabilities (used to
865 determine if the default took place and subsequently compute its
866 event time)
867 */
868 Date today = Settings::instance().evaluationDate();
869 Date maxHorizonDate = today + Period(this->maxHorizon_, Days);
870
871 const ext::shared_ptr<Pool>& pool = this->basket_->pool();
872 for(Size iName=0; iName < this->basket_->size(); ++iName)//use'live'
873 horizonDefaultPs_.push_back(pool->get(name: pool->names()[iName]).
874 defaultProbability(key: this->basket_->defaultKeys()[iName])
875 ->defaultProbability(maxHorizonDate, true));
876 }
877 Real getEventRecovery(const defaultSimEvent& evt) const {
878 return recoveries_[evt.nameIdx];
879 }
880 Real expectedRecovery(const Date&, Size iName, const DefaultProbKey&) const override {
881 // deterministic
882 return recoveries_[iName];
883 }
884
885 Real latentVarValue(const std::vector<Real>& factorsSample,
886 Size iVar) const {
887 return model_->latentVarValue(factorsSample, iVar);
888 }
889 //allows statistics to know the portfolio size (could be moved to parent
890 //invoking duck typing on the variable name or a handle to the basket)
891 Size basketSize() const { return model_->size(); }
892 private:
893 void resetModel() override /*const*/ {
894 /* Explore: might save recalculation if the basket is the same
895 (some situations, like BC or control variates) in that case do not
896 update, only reset the model's basket.
897 */
898 model_->resetBasket(this->basket_.currentLink());
899
900 QL_REQUIRE(this->basket_->size() == model_->size(),
901 "Incompatible basket and model sizes.");
902 QL_REQUIRE(recoveries_.size() == this->basket_->size(),
903 "Incompatible basket and recovery sizes.");
904 // invalidate current calculations if any and notify observers
905 // NOLINTNEXTLINE(bugprone-parent-virtual-call)
906 LazyObject::update();
907 }
908 // This one and the buffer might be moved to the parent, only some
909 // dates might be specific to a particular model.
910 // Default probabilities for each name at the time of the maximun
911 // horizon date. Cached for perf.
912 mutable std::vector<Probability> horizonDefaultPs_;
913 };
914
915
916
917
918
919 template<class C, class URNG>
920 void RandomDefaultLM<C, URNG>::nextSample(
921 const std::vector<Real>& values) const
922 {
923 const ext::shared_ptr<Pool>& pool = this->basket_->pool();
924 // starts with no events
925 this->simsBuffer_.push_back(std::vector<defaultSimEvent> ());
926
927 for(Size iName=0; iName<model_->size(); iName++) {
928 Real latentVarSample =
929 model_->latentVarValue(values, iName);
930 Probability simDefaultProb =
931 model_->cumulativeY(latentVarSample, iName);
932 // If the default simulated lies before the max date:
933 if (horizonDefaultPs_[iName] >= simDefaultProb) {
934 const Handle<DefaultProbabilityTermStructure>& dfts =
935 pool->get(name: pool->names()[iName]).// use 'live' names
936 defaultProbability(key: this->basket_->defaultKeys()[iName]);
937 // compute and store default time with respect to the
938 // curve ref date:
939 Size dateSTride =
940 static_cast<Size>(Brent().solve(// casted from Real:
941 f: detail::Root(dfts, simDefaultProb),
942 accuracy: accuracy_,guess: 0.,step: 1.));
943 /*
944 // value if one approximates to a flat HR;
945 // faster (>x2) but it introduces an error:..
946 // \todo: see how to include this 'polymorphically'.
947 // While not the case in pricing in risk metrics/real
948 // probabilities the curves are often flat
949 static_cast<Size>(ceil(maxHorizon_ *
950 std::log(1.-simDefaultProb)
951 /std::log(1.-data_.horizonDefaultPs_[iName])));
952 */
953 this->simsBuffer_.back().push_back(defaultSimEvent(iName,
954 dateSTride));
955 //emplace_back
956 }
957 /* Used to remove sims with no events. Uses less memory, faster
958 post-statistics. But only if all names in the portfolio have low
959 default probability, otherwise is more expensive and sim access has
960 to be modified. However low probability is also an indicator that
961 variance reduction is needed. */
962 }
963 }
964
965
966
967
968 // Common usage typedefs (notice they vary in the multithread version)
969 // ---------- Gaussian default generators options ------------------------
970 /* Uses copula direct normal inversion and MT generator
971 typedef RandomDefaultLM<GaussianCopulaPolicy,
972 RandomSequenceGenerator<MersenneTwisterUniformRng> >
973 GaussianRandomDefaultLM;
974 */
975 /* Uses BoxMuller for gaussian generation, bypassing copula inversions
976 typedef RandomDefaultLM<GaussianCopulaPolicy, RandomSequenceGenerator<
977 BoxMullerGaussianRng<MersenneTwisterUniformRng> > >
978 GaussianRandomDefaultLM;
979 */
980 /* Default case, uses the copula inversion directly and sobol sequence */
981 typedef RandomDefaultLM<GaussianCopulaPolicy> GaussianRandomDefaultLM;
982
983 // ---------- T default generators options ----------------------------
984 /* Uses copula inversion and MT base generation
985 typedef RandomDefaultLM<TCopulaPolicy,
986 RandomSequenceGenerator<MersenneTwisterUniformRng> > TRandomDefaultLM;
987 */
988 /* Uses MT and polar direct strudent-T generation
989 typedef RandomDefaultLM<TCopulaPolicy,
990 RandomSequenceGenerator<PolarStudentTRng<MersenneTwisterUniformRng> > >
991 TRandomDefaultLM;
992 */
993 /* Default case, uses sobol sequence and copula inversion */
994 typedef RandomDefaultLM<TCopulaPolicy> TRandomDefaultLM;
995
996}
997
998#endif
999

source code of quantlib/ql/experimental/credit/randomdefaultlatentmodel.hpp