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/onefactorcopula.hpp>
21
22using namespace std;
23
24namespace QuantLib {
25
26 //-------------------------------------------------------------------------
27 Real OneFactorCopula::conditionalProbability(Real p, Real m) const {
28 //-------------------------------------------------------------------------
29 calculate ();
30 // FIXME
31 if (p < 1e-10) return 0;
32
33 Real c = correlation_->value();
34
35 Real res = cumulativeZ (z: (inverseCumulativeY (p) - sqrt(x: c) * m)
36 / sqrt (x: 1. - c));
37
38 QL_REQUIRE (res >= 0 && res <= 1,
39 "conditional probability " << res << "out of range");
40
41 return res;
42 }
43
44 //-------------------------------------------------------------------------
45 vector<Real> OneFactorCopula::conditionalProbability(
46 const vector<Real>& prob,
47 Real m) const {
48 //-------------------------------------------------------------------------
49 calculate ();
50 vector<Real> p (prob.size(), 0);
51 for (Size i = 0; i < p.size(); i++)
52 p[i] = conditionalProbability (p: prob[i], m);
53 return p;
54 }
55
56 //-------------------------------------------------------------------------
57 Real OneFactorCopula::cumulativeY (Real y) const {
58 //-------------------------------------------------------------------------
59 calculate ();
60
61 QL_REQUIRE(!y_.empty(), "cumulative Y not tabulated yet");
62
63 // linear interpolation on the tabulated cumulative distribution of Y
64 if (y < y_.front())
65 return cumulativeY_.front();
66
67 for (Size i = 0; i < y_.size(); i++) {
68 if (y_[i] > y)
69 return ( (y_[i] - y) * cumulativeY_[i-1]
70 + (y - y_[i-1]) * cumulativeY_[i] )
71 / (y_[i] - y_[i-1]);
72 }
73
74 return cumulativeY_.back();
75 }
76
77 //-------------------------------------------------------------------------
78 Real OneFactorCopula::inverseCumulativeY (Real x) const {
79 //-------------------------------------------------------------------------
80 calculate ();
81
82 QL_REQUIRE(!y_.empty(), "cumulative Y not tabulated yet");
83
84 // linear interpolation on the tabulated cumulative distribution of Y
85 if (x < cumulativeY_.front())
86 return y_.front();
87
88 for (Size i = 0; i < cumulativeY_.size(); i++) {
89 if (cumulativeY_[i] > x)
90 return ( (cumulativeY_[i] - x) * y_[i-1]
91 + (x - cumulativeY_[i-1]) * y_[i] )
92 / (cumulativeY_[i] - cumulativeY_[i-1]);
93 }
94
95 return y_.back();
96 }
97
98 //-------------------------------------------------------------------------
99 int OneFactorCopula::checkMoments (Real tolerance) const {
100 //-------------------------------------------------------------------------
101 calculate ();
102
103 Real norm = 0, mean = 0, var = 0;
104 for (Size i = 0; i < steps(); i++) {
105 norm += densitydm (i);
106 mean += m(i) * densitydm (i);
107 var += pow (x: m(i), y: 2) * densitydm (i);
108 }
109
110 QL_REQUIRE (fabs (norm - 1.0) < tolerance, "norm out of tolerance range");
111 QL_REQUIRE (fabs (mean) < tolerance, "mean out of tolerance range");
112 QL_REQUIRE (fabs (var - 1.0) < tolerance, "variance out of tolerance range");
113
114 // FIXME: define range for Y via cutoff quantil?
115 Real zMin = -10;
116 Real zMax = +10;
117 Size zSteps = 200;
118 norm = 0;
119 mean = 0;
120 var = 0;
121 for (Size i = 1; i < zSteps; i++) {
122 Real z1 = zMin + (zMax - zMin) / zSteps * (i - 1);
123 Real z2 = zMin + (zMax - zMin) / zSteps * i;
124 Real z = (z1 + z2) / 2;
125 Real densitydz = cumulativeZ (z: z2) - cumulativeZ (z: z1);
126 norm += densitydz;
127 mean += z * densitydz;
128 var += pow (x: z, y: 2) * densitydz;
129 }
130
131 QL_REQUIRE (fabs (norm - 1.0) < tolerance, "norm out of tolerance range");
132 QL_REQUIRE (fabs (mean) < tolerance, "mean out of tolerance range");
133 QL_REQUIRE (fabs (var - 1.0) < tolerance, "variance out of tolerance range");
134
135 // FIXME: define range for Y via cutoff quantil?
136 Real yMin = -10;
137 Real yMax = +10;
138 Size ySteps = 200;
139 norm = 0;
140 mean = 0;
141 var = 0;
142 for (Size i = 1; i < ySteps; i++) {
143 Real y1 = yMin + (yMax - yMin) / ySteps * (i - 1);
144 Real y2 = yMin + (yMax - yMin) / ySteps * i;
145 Real y = (y1 + y2) / 2;
146 Real densitydy = cumulativeY (y: y2) - cumulativeY (y: y1);
147 norm += densitydy;
148 mean += y * densitydy;
149 var += y * y * densitydy;
150 }
151
152 QL_REQUIRE (fabs (norm - 1.0) < tolerance, "norm out of tolerance range");
153 QL_REQUIRE (fabs (mean) < tolerance, "mean out of tolerance range");
154 QL_REQUIRE (fabs (var - 1.0) < tolerance, "variance out of tolerance range");
155
156 return 0;
157 }
158
159}
160
161

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