1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2014 Jose Aparicio
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/math/convolvedstudentt.hpp>
21#include <ql/errors.hpp>
22#include <ql/math/factorial.hpp>
23#include <ql/math/distributions/normaldistribution.hpp>
24#include <ql/math/solvers1d/brent.hpp>
25#include <ql/math/functional.hpp>
26#include <boost/math/distributions/students_t.hpp>
27
28namespace QuantLib {
29
30 CumulativeBehrensFisher::CumulativeBehrensFisher(const std::vector<Integer>& degreesFreedom,
31 const std::vector<Real>& factors)
32 : degreesFreedom_(degreesFreedom), factors_(factors), polyConvolved_(std::vector<Real>(1, 1.))
33
34 {
35 QL_REQUIRE(degreesFreedom.size() == factors.size(),
36 "Incompatible sizes in convolution.");
37 for (int i : degreesFreedom) {
38 QL_REQUIRE(i % 2 != 0, "Even degree of freedom not allowed");
39 QL_REQUIRE(i >= 0, "Negative degree of freedom not allowed");
40 }
41 for(Size i=0; i<degreesFreedom_.size(); i++)
42 polynCharFnc_.push_back(x: polynCharactT(n: (degreesFreedom[i]-1)/2));
43 // adjust the polynomial coefficients by the factors in the linear
44 // combination:
45 for(Size i=0; i<degreesFreedom_.size(); i++) {
46 Real multiplier = 1.;
47 for(Size k=1; k<polynCharFnc_[i].size(); k++) {
48 multiplier *= std::abs(x: factors_[i]);
49 polynCharFnc_[i][k] *= multiplier;
50 }
51 }
52 //convolution, here it is a product of polynomials and exponentials
53 for (auto& i : polynCharFnc_)
54 polyConvolved_ = convolveVectorPolynomials(v1: polyConvolved_, v2: i);
55 // trim possible zeros that might have arised:
56 auto it = polyConvolved_.rbegin();
57 while (it != polyConvolved_.rend()) {
58 if (*it == 0.) {
59 polyConvolved_.pop_back();
60 it = polyConvolved_.rbegin();
61 }else{
62 break;
63 }
64 }
65 // cache 'a' value (the exponent)
66 for(Size i=0; i<degreesFreedom_.size(); i++)
67 a_ += std::sqrt(x: static_cast<Real>(degreesFreedom_[i]))
68 * std::abs(x: factors_[i]);
69 a2_ = a_ * a_;
70 }
71
72 std::vector<Real> CumulativeBehrensFisher::polynCharactT(Natural n) const {
73 Natural nu = 2 * n +1;
74 std::vector<Real> low(1,1.), high(1,1.);
75 high.push_back(x: std::sqrt(x: static_cast<Real>(nu)));
76 if(n==0) return low;
77 if(n==1) return high;
78
79 for(Size k=1; k<n; k++) {
80 std::vector<Real> recursionFactor(1,0.); // 0 coef
81 recursionFactor.push_back(x: 0.); // 1 coef
82 recursionFactor.push_back(x: nu/((2.*k+1.)*(2.*k-1.))); // 2 coef
83 std::vector<Real> lowUp =
84 convolveVectorPolynomials(v1: recursionFactor, v2: low);
85 //add them up:
86 for(Size i=0; i<high.size(); i++)
87 lowUp[i] += high[i];
88 low = high;
89 high = lowUp;
90 }
91 return high;
92 }
93
94 std::vector<Real> CumulativeBehrensFisher::convolveVectorPolynomials(
95 const std::vector<Real>& v1,
96 const std::vector<Real>& v2) const {
97 #if defined(QL_EXTRA_SAFETY_CHECKS)
98 QL_REQUIRE(!v1.empty() && !v2.empty(),
99 "Incorrect vectors in polynomial.");
100 #endif
101
102 const std::vector<Real>& shorter = v1.size() < v2.size() ? v1 : v2;
103 const std::vector<Real>& longer = (v1 == shorter) ? v2 : v1;
104
105 Size newDegree = v1.size()+v2.size()-2;
106 std::vector<Real> resultB(newDegree+1, 0.);
107 for(Size polyOrdr=0; polyOrdr<resultB.size(); polyOrdr++) {
108 for(Size i=std::max<Integer>(a: 0, b: polyOrdr-longer.size()+1);
109 i<=std::min(a: polyOrdr, b: shorter.size()-1); i++)
110 resultB[polyOrdr] += shorter[i]*longer[polyOrdr-i];
111 }
112 return resultB;
113 }
114
115 Probability CumulativeBehrensFisher::operator()(const Real x) const {
116 // 1st & 0th terms with the table integration
117 Real integral = polyConvolved_[0] * std::atan(x: x/a_);
118 Real squared = a2_ + x*x;
119 Real rootsqr = std::sqrt(x: squared);
120 Real atan2xa = std::atan2(y: -x,x: a_);
121 if(polyConvolved_.size()>1)
122 integral += polyConvolved_[1] * x/squared;
123
124 for(Size exponent = 2; exponent <polyConvolved_.size(); exponent++) {
125 integral -= polyConvolved_[exponent] *
126 Factorial::get(n: exponent-1) * std::sin(x: (exponent)*atan2xa)
127 /std::pow(x: rootsqr, y: static_cast<Real>(exponent));
128 }
129 return .5 + integral / M_PI;
130 }
131
132 Probability
133 CumulativeBehrensFisher::density(const Real x) const {
134 Real squared = a2_ + x*x;
135 Real integral = polyConvolved_[0] * a_ / squared;
136 Real rootsqr = std::sqrt(x: squared);
137 Real atan2xa = std::atan2(y: -x,x: a_);
138 for(Size exponent=1; exponent <polyConvolved_.size(); exponent++) {
139 integral += polyConvolved_[exponent] *
140 Factorial::get(n: exponent) * std::cos(x: (exponent+1)*atan2xa)
141 /std::pow(x: rootsqr, y: static_cast<Real>(exponent+1) );
142 }
143 return integral / M_PI;
144 }
145
146
147
148 InverseCumulativeBehrensFisher::InverseCumulativeBehrensFisher(
149 const std::vector<Integer>& degreesFreedom,
150 const std::vector<Real>& factors,
151 Real accuracy)
152 : normSqr_(std::inner_product(first1: factors.begin(), last1: factors.end(),
153 first2: factors.begin(), init: Real(0.))),
154 accuracy_(accuracy), distrib_(degreesFreedom, factors) { }
155
156 Real InverseCumulativeBehrensFisher::operator()(const Probability q) const {
157 Probability effectiveq;
158 Real sign;
159 // since the distrib is symmetric solve only on the right side:
160 if(q==0.5) {
161 return 0.;
162 }else if(q < 0.5) {
163 sign = -1.;
164 effectiveq = 1.-q;
165 }else{
166 sign = 1.;
167 effectiveq = q;
168 }
169 Real xMin =
170 InverseCumulativeNormal::standard_value(x: effectiveq) * normSqr_;
171 // inversion will fail at the Brent's bounds-check if this is not enough
172 // (q is very close to 1.), in a bad combination fails around 1.-1.e-7
173 Real xMax = 1.e6;
174 return sign *
175 Brent().solve(f: [&](Real x) -> Real { return distrib_(x) - effectiveq; },
176 accuracy: accuracy_, guess: (xMin+xMax)/2., xMin, xMax);
177 }
178
179}
180

source code of quantlib/ql/experimental/math/convolvedstudentt.cpp