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 22 namespace QuantLib { 23 24 namespace { 25 26 Real epsilon = 1.0e-8; 27 sign(Real val)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 */ arctan(Real x,Real y)36 Real arctan(Real x, Real y) { 37 Real res = std::atan2(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 f_x(Real m,Real h,Real k,Real rho)42 Real f_x(Real m, Real h, Real k, Real rho) { 43 Real unCor = 1 - rho*rho; 44 Real sub = std::pow(h - rho * k, 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 P_n(Real h,Real k,Natural n,Real rho)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(n * M_PI); 56 Real xHK = f_x(n, h, k, rho); 57 Real xKH = f_x(n, k, h, rho); 58 Real divH = 1 + h*h / n; 59 Real divK = 1 + k*k / n; 60 Real sgnHK = sign(h - rho * k); 61 Real sgnKH = sign(k - rho * h); 62 63 if (n % 2 == 0) { // n is even, equation (10) 64 // first line of (10) 65 Real res = arctan(std::sqrt(unCor), -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(std::sqrt(xHK), std::sqrt(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(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(std::sqrt(xKH), std::sqrt(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(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(h*h - 2 * rho * hk + k*k + n * unCor); 119 Real res = arctan(std::sqrt(Real(n)) * (-(h + k) * hkcn - (hk - n) * sqrtExpr), 120 (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(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(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:: BivariateCumulativeStudentDistribution(Natural n,Real rho)163 BivariateCumulativeStudentDistribution(Natural n, 164 Real rho) 165 : n_(n), rho_(rho) {} 166 operator ()(Real x,Real y) const167 Real BivariateCumulativeStudentDistribution::operator()(Real x, 168 Real y) const { 169 return P_n(x, y, n_, rho_); 170 } 171 172 } 173