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/onefactorstudentcopula.hpp>
21
22namespace QuantLib {
23
24 //-------------------------------------------------------------------------
25 OneFactorStudentCopula::OneFactorStudentCopula (
26 const Handle<Quote>& correlation,
27 int nz, int nm,
28 Real maximum,
29 Size integrationSteps)
30 : OneFactorCopula (correlation, maximum, integrationSteps),
31 density_ (nm), cumulative_ (nz), nz_(nz), nm_(nm) {
32 //-------------------------------------------------------------------------
33
34 QL_REQUIRE (nz > 2 && nm > 2, "degrees of freedom must be > 2");
35
36 scaleM_ = std::sqrt (x: Real (nm_ - 2) / nm_);
37 scaleZ_ = std::sqrt (x: Real (nz_ - 2) / nz_);
38
39 calculate ();
40 }
41
42 //-------------------------------------------------------------------------
43 void OneFactorStudentCopula::performCalculations () const {
44 //-------------------------------------------------------------------------
45 y_.clear();
46 cumulativeY_.clear();
47
48 // FIXME:
49 // compute F(ymin) and F(ymax) for the fattest case nm = nz = 2
50 // set a desired confidence and work out ymin, ymax
51 Real ymin = -10;
52 Real ymax = +10;
53 Size steps = 200;
54 for (Size i = 0; i <= steps; i++) {
55 Real y = ymin + (ymax - ymin) * i / steps;
56 Real c = cumulativeYintegral (y);
57 y_.push_back (x: y);
58 cumulativeY_.push_back (x: c);
59 }
60 }
61
62 //-------------------------------------------------------------------------
63 Real OneFactorStudentCopula::cumulativeYintegral (Real y) const {
64 //-------------------------------------------------------------------------
65 Real c = correlation_->value();
66
67 if (c == 0)
68 return CumulativeStudentDistribution(nz_)(y / scaleZ_);
69
70 if (c == 1)
71 return CumulativeStudentDistribution(nm_)(y / scaleM_);
72
73 StudentDistribution dz (nz_);
74 StudentDistribution dm (nm_);
75
76 // FIXME:
77 // Find a sensitive way of setting these parameters,
78 // e.g. depending on nm and nz, and the desired table range
79 Real minimum = -10; // -15
80 Real maximum = +10; // +15
81 int steps = 400;
82
83 Real delta = (maximum - minimum) / steps;
84 Real cumulated = 0;
85
86 if (c < 0.5) {
87 // outer integral -> 1 for c -> 0
88 // inner integral -> cumulativeStudent(nz)(y) for c-> 0
89 for (Real m = minimum + delta/2; m < maximum; m += delta)
90 for (Real z = minimum + delta/2;
91 z < (y - std::sqrt(x: c) * m) / std::sqrt (x: 1. - c); z += delta)
92 cumulated += dm (m / scaleM_) / scaleM_
93 * dz (z / scaleZ_) / scaleZ_;
94 }
95 else {
96 // outer integral -> 1 for c -> 1
97 // inner integral -> cumulativeStudent(nm)(y) for c-> 1
98 for (Real z = minimum + delta/2; z < maximum; z += delta)
99 for (Real m = minimum + delta/2;
100 m < (y - std::sqrt(x: 1.0 - c) * z) / std::sqrt(x: c); m += delta)
101 cumulated += dm (m / scaleM_) / scaleM_
102 * dz (z / scaleZ_) / scaleZ_;
103 }
104
105 return cumulated * delta * delta;
106 }
107
108 //-------------------------------------------------------------------------
109 OneFactorGaussianStudentCopula::OneFactorGaussianStudentCopula (
110 const Handle<Quote>& correlation,
111 int nz, Real maximum,
112 Size integrationSteps)
113 : OneFactorCopula (correlation, maximum, integrationSteps),
114 cumulative_(nz), nz_(nz) {
115 //-------------------------------------------------------------------------
116
117 QL_REQUIRE (nz > 2, "degrees of freedom must be > 2");
118
119 scaleZ_ = std::sqrt (x: Real (nz_ - 2) / nz_);
120
121 calculate ();
122 }
123
124 //-------------------------------------------------------------------------
125 void OneFactorGaussianStudentCopula::performCalculations () const {
126 //-------------------------------------------------------------------------
127 y_.clear();
128 cumulativeY_.clear();
129
130 // FIXME:
131 // compute F(ymin) and F(ymax) for the fattest case nm = nz = 2
132 // set a desired confidence and work out ymin, ymax
133 Real ymin = -10;
134 Real ymax = +10;
135 Size steps = 200;
136 for (Size i = 0; i <= steps; i++) {
137 Real y = ymin + (ymax - ymin) * i / steps;
138 Real c = cumulativeYintegral (y);
139 y_.push_back (x: y);
140 cumulativeY_.push_back (x: c);
141 }
142 }
143
144 //-------------------------------------------------------------------------
145 Real OneFactorGaussianStudentCopula::cumulativeYintegral (Real y) const {
146 //-------------------------------------------------------------------------
147 Real c = correlation_->value();
148
149 if (c == 0)
150 return CumulativeStudentDistribution(nz_)(y / scaleZ_);
151
152 if (c == 1)
153 return CumulativeNormalDistribution()(y);
154
155 StudentDistribution dz (nz_);
156 NormalDistribution dm;
157
158 // FIXME:
159 // Find a sensitive way of setting these parameters,
160 // e.g. depending on nm and nz, and the desired table range
161 Real minimum = -10;
162 Real maximum = +10;
163 int steps = 400;
164
165 Real delta = (maximum - minimum) / steps;
166 Real cumulated = 0;
167
168 if (c < 0.5) {
169 // outer integral -> 1 for c -> 0
170 // inner integral -> cumulativeStudent(nz)(y) for c-> 0
171 for (Real m = minimum + delta/2; m < maximum; m += delta)
172 for (Real z = minimum + delta/2;
173 z < (y - std::sqrt(x: c) * m) / std::sqrt (x: 1. - c);
174 z += delta)
175 cumulated += dm (m) * dz (z / scaleZ_) / scaleZ_;
176 }
177 else {
178 // outer integral -> 1 for c -> 1
179 // inner integral -> cumulativeNormal(y) for c-> 1
180 for (Real z = minimum + delta/2; z < maximum; z += delta)
181 for (Real m = minimum + delta/2;
182 m < (y - std::sqrt(x: 1.0 - c) * z) / std::sqrt(x: c);
183 m += delta)
184 cumulated += dm (m) * dz (z / scaleZ_) / scaleZ_;
185 }
186
187 return cumulated * delta * delta;
188 }
189
190 //-------------------------------------------------------------------------
191 OneFactorStudentGaussianCopula::OneFactorStudentGaussianCopula (
192 const Handle<Quote>& correlation,
193 int nm, Real maximum,
194 Size integrationSteps)
195 : OneFactorCopula (correlation, maximum, integrationSteps),
196 density_ (nm), nm_(nm) {
197 //-------------------------------------------------------------------------
198
199 QL_REQUIRE (nm > 2, "degrees of freedom must be > 2");
200
201 scaleM_ = std::sqrt (x: Real (nm_ - 2) / nm_);
202
203 calculate ();
204 }
205
206 //-------------------------------------------------------------------------
207 void OneFactorStudentGaussianCopula::performCalculations () const {
208 //-------------------------------------------------------------------------
209 y_.clear();
210 cumulativeY_.clear();
211
212 // FIXME:
213 // compute F(ymin) and F(ymax) for the fattest case nm = nz = 2
214 // set a desired confidence and work out ymin, ymax
215 Real ymin = -10;
216 Real ymax = +10;
217 Size steps = 200;
218 for (Size i = 0; i <= steps; i++) {
219 Real y = ymin + (ymax - ymin) * i / steps;
220 Real c = cumulativeYintegral (y);
221 y_.push_back (x: y);
222 cumulativeY_.push_back (x: c);
223 }
224 }
225
226 //-------------------------------------------------------------------------
227 Real OneFactorStudentGaussianCopula::cumulativeYintegral (Real y) const {
228 //-------------------------------------------------------------------------
229 Real c = correlation_->value();
230
231 if (c == 0)
232 return CumulativeNormalDistribution()(y);
233
234 if (c == 1)
235 return CumulativeStudentDistribution(nm_)(y / scaleM_);
236
237
238 StudentDistribution dm (nm_);
239 NormalDistribution dz;
240
241 // FIXME:
242 // Find a sensitive way of setting these parameters,
243 // e.g. depending on nm and nz, and the desired table range
244 Real minimum = -10;
245 Real maximum = +10;
246 int steps = 400;
247
248 Real delta = (maximum - minimum) / steps;
249 Real cumulated = 0;
250
251 if (c < 0.5) {
252 // outer integral -> 1 for c -> 0
253 // inner integral -> cumulativeNormal(y) for c-> 0
254 for (Real m = minimum + delta/2; m < maximum; m += delta)
255 for (Real z = minimum + delta/2;
256 z < (y - std::sqrt(x: c) * m) / std::sqrt (x: 1. - c);
257 z += delta)
258 cumulated += dm (m / scaleM_) / scaleM_ * dz (z);
259 }
260 else {
261 // outer integral -> 1 for c -> 1
262 // inner integral -> cumulativeStudent(nm)(y) for c-> 1
263 for (Real z = minimum + delta/2; z < maximum; z += delta)
264 for (Real m = minimum + delta/2;
265 m < (y - std::sqrt(x: 1.0 - c) * z) / std::sqrt(x: c);
266 m += delta)
267 cumulated += dm (m / scaleM_) / scaleM_ * dz (z);
268 }
269
270 return cumulated * delta * delta;
271 }
272
273}
274

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