1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2014 Michal Kaut
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/math/distributions/bivariatestudenttdistribution.hpp>
21
22namespace QuantLib {
23
24 namespace {
25
26 Real epsilon = 1.0e-8;
27
28 Real sign(Real val) {
29 return val == 0.0 ? 0.0
30 : (val < 0.0 ? -1.0 : 1.0);
31 }
32
33 /* unlike the atan2 function in C++ that gives results in
34 [-pi,pi], this returns a value in [0, 2*pi]
35 */
36 Real arctan(Real x, Real y) {
37 Real res = std::atan2(y: x, x: y);
38 return res >= 0.0 ? res : res + 2 * M_PI;
39 }
40
41 // function x(m,h,k) defined on top of page 155
42 Real f_x(Real m, Real h, Real k, Real rho) {
43 Real unCor = 1 - rho*rho;
44 Real sub = std::pow(x: h - rho * k, y: 2);
45 Real denom = sub + unCor * (m + k*k);
46 if (denom < epsilon)
47 return 0.0; // limit case for rho = +/-1.0
48 return sub / (sub + unCor * (m + k*k));
49 }
50
51 // this calculates the cdf
52 Real P_n(Real h, Real k, Natural n, Real rho) {
53 Real unCor = 1.0 - rho*rho;
54
55 Real div = 4 * std::sqrt(x: n * M_PI);
56 Real xHK = f_x(m: n, h, k, rho);
57 Real xKH = f_x(m: n, h: k, k: h, rho);
58 Real divH = 1 + h*h / n;
59 Real divK = 1 + k*k / n;
60 Real sgnHK = sign(val: h - rho * k);
61 Real sgnKH = sign(val: k - rho * h);
62
63 if (n % 2 == 0) { // n is even, equation (10)
64 // first line of (10)
65 Real res = arctan(x: std::sqrt(x: unCor), y: -rho) / M_TWOPI;
66
67 // second line of (10)
68 Real dgM = 2 * (1 - xHK); // multiplier for dgj
69 Real gjM = sgnHK * 2 / M_PI; // multiplier for g_j
70 // initializations for j = 1:
71 Real f_j = std::sqrt(M_PI / divK);
72 Real g_j = 1 + gjM * arctan(x: std::sqrt(x: xHK), y: std::sqrt(x: 1 - xHK));
73 Real sum = f_j * g_j;
74 if (n >= 4) {
75 // different formulas for j = 2:
76 f_j *= 0.5 / divK; // (2 - 1.5) / (Real) (2 - 1) / divK;
77 Real dgj = gjM * std::sqrt(x: xHK * (1 - xHK));
78 g_j += dgj;
79 sum += f_j * g_j;
80 // and then the loop for the rest of the j's:
81 for (Natural j = 3; j <= n / 2; ++j) {
82 f_j *= (j - 1.5) / (Real) (j - 1) / divK;
83 dgj *= (Real) (j - 2) / (2 * j - 3) * dgM;
84 g_j += dgj;
85 sum += f_j * g_j;
86 }
87 }
88 res += k / div * sum;
89
90 // third line of (10)
91 dgM = 2 * (1 - xKH);
92 gjM = sgnKH * 2 / M_PI;
93 // initializations for j = 1:
94 f_j = std::sqrt(M_PI / divH);
95 g_j = 1 + gjM * arctan(x: std::sqrt(x: xKH), y: std::sqrt(x: 1 - xKH));
96 sum = f_j * g_j;
97 if (n >= 4) {
98 // different formulas for j = 2:
99 f_j *= 0.5 / divH; // (2 - 1.5) / (Real) (2 - 1) / divK;
100 Real dgj = gjM * std::sqrt(x: xKH * (1 - xKH));
101 g_j += dgj;
102 sum += f_j * g_j;
103 // and then the loop for the rest of the j's:
104 for (Natural j = 3; j <= n / 2; ++j) {
105 f_j *= (j - 1.5) / (Real) (j - 1) / divH;
106 dgj *= (Real) (j - 2) / (2 * j - 3) * dgM;
107 g_j += dgj;
108 sum += f_j * g_j;
109 }
110 }
111 res += h / div * sum;
112 return res;
113
114 } else { // n is odd, equation (11)
115 // first line of (11)
116 Real hk = h * k;
117 Real hkcn = hk + rho * n;
118 Real sqrtExpr = std::sqrt(x: h*h - 2 * rho * hk + k*k + n * unCor);
119 Real res = arctan(x: std::sqrt(x: Real(n)) * (-(h + k) * hkcn - (hk - n) * sqrtExpr),
120 y: (hk - n) * hkcn - n * (h + k) * sqrtExpr ) / M_TWOPI;
121
122 if (n > 1) {
123 // second line of (11)
124 Real mult = (1 - xHK) / 2;
125 // initializations for j = 1:
126 Real f_j = 2 / std::sqrt(M_PI) / divK;
127 Real dgj = sgnHK * std::sqrt(x: xHK);
128 Real g_j = 1 + dgj;
129 Real sum = f_j * g_j;
130 // and then the loop for the rest of the j's:
131 for (Natural j = 2; j <= (n - 1) / 2; ++j) {
132 f_j *= (Real) (j - 1) / (j - 0.5) / divK;
133 dgj *= (Real) (2 * j - 3) / (j - 1) * mult;
134 g_j += dgj;
135 sum += f_j * g_j;
136 }
137 res += k / div * sum;
138
139 // third line of (11)
140 mult = (1 - xKH) / 2;
141 // initializations for j = 1:
142 f_j = 2 / std::sqrt(M_PI) / divH;
143 dgj = sgnKH * std::sqrt(x: xKH);
144 g_j = 1 + dgj;
145 sum = f_j * g_j;
146 // and then the loop for the rest of the j's:
147 for (Natural j = 2; j <= (n - 1) / 2; ++j) {
148 f_j *= (Real) (j - 1) / (j - 0.5) / divH;
149 dgj *= (Real) (2 * j - 3) / (j - 1) * mult;
150 g_j += dgj;
151 sum += f_j * g_j;
152 }
153 res += h / div * sum;
154 }
155 return res;
156 }
157 }
158
159 }
160
161
162 BivariateCumulativeStudentDistribution::
163 BivariateCumulativeStudentDistribution(Natural n,
164 Real rho)
165 : n_(n), rho_(rho) {}
166
167 Real BivariateCumulativeStudentDistribution::operator()(Real x,
168 Real y) const {
169 return P_n(h: x, k: y, n: n_, rho: rho_);
170 }
171
172}
173

source code of quantlib/ql/math/distributions/bivariatestudenttdistribution.cpp