QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
convolvedstudentt.cpp
Go to the documentation of this file.
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
21#include <ql/errors.hpp>
22#include <ql/math/factorial.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(polynCharactT((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(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_)
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(static_cast<Real>(degreesFreedom_[i]))
68 * std::abs(factors_[i]);
69 a2_ = a_ * a_;
70 }
71
73 Natural nu = 2 * n +1;
74 std::vector<Real> low(1,1.), high(1,1.);
75 high.push_back(std::sqrt(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(0.); // 1 coef
82 recursionFactor.push_back(nu/((2.*k+1.)*(2.*k-1.))); // 2 coef
83 std::vector<Real> lowUp =
84 convolveVectorPolynomials(recursionFactor, 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
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>(0, polyOrdr-longer.size()+1);
109 i<=std::min(polyOrdr, shorter.size()-1); i++)
110 resultB[polyOrdr] += shorter[i]*longer[polyOrdr-i];
111 }
112 return resultB;
113 }
114
116 // 1st & 0th terms with the table integration
117 Real integral = polyConvolved_[0] * std::atan(x/a_);
118 Real squared = a2_ + x*x;
119 Real rootsqr = std::sqrt(squared);
120 Real atan2xa = std::atan2(-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(exponent-1) * std::sin((exponent)*atan2xa)
127 /std::pow(rootsqr, static_cast<Real>(exponent));
128 }
129 return .5 + integral / M_PI;
130 }
131
134 Real squared = a2_ + x*x;
135 Real integral = polyConvolved_[0] * a_ / squared;
136 Real rootsqr = std::sqrt(squared);
137 Real atan2xa = std::atan2(-x,a_);
138 for(Size exponent=1; exponent <polyConvolved_.size(); exponent++) {
139 integral += polyConvolved_[exponent] *
140 Factorial::get(exponent) * std::cos((exponent+1)*atan2xa)
141 /std::pow(rootsqr, static_cast<Real>(exponent+1) );
142 }
143 return integral / M_PI;
144 }
145
146
147
149 const std::vector<Integer>& degreesFreedom,
150 const std::vector<Real>& factors,
151 Real accuracy)
152 : normSqr_(std::inner_product(factors.begin(), factors.end(),
153 factors.begin(), Real(0.))),
154 accuracy_(accuracy), distrib_(degreesFreedom, factors) { }
155
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 =
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([&](Real x) -> Real { return distrib_(x) - effectiveq; },
176 accuracy_, (xMin+xMax)/2., xMin, xMax);
177 }
178
179}
Brent 1-D solver.
Brent 1-D solver
Definition: brent.hpp:37
std::vector< std::vector< Real > > polynCharFnc_
const std::vector< Real > & factors() const
Factors in the linear combination.
Probability density(Real x) const
Returns the probability density of the resulting distribution.
Probability operator()(Real x) const
Returns the cumulative probability of the resulting distribution.
CumulativeBehrensFisher(const std::vector< Integer > &degreesFreedom=std::vector< Integer >(), const std::vector< Real > &factors=std::vector< Real >())
std::vector< Real > convolveVectorPolynomials(const std::vector< Real > &v1, const std::vector< Real > &v2) const
std::vector< Integer > degreesFreedom_
std::vector< Real > polynCharactT(Natural n) const
Student t characteristic polynomials.
static Real get(Natural n)
Definition: factorial.cpp:50
Real operator()(Probability q) const
Returns the cumulative inverse value.
InverseCumulativeBehrensFisher(const std::vector< Integer > &degreesFreedom=std::vector< Integer >(), const std::vector< Real > &factors=std::vector< Real >(), Real accuracy=1.e-6)
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Factorial numbers calculator.
QL_REAL Real
real number
Definition: types.hpp:50
unsigned QL_INTEGER Natural
positive integer
Definition: types.hpp:43
Real Probability
probability
Definition: types.hpp:82
std::size_t Size
size of a container
Definition: types.hpp:58
functionals and combinators not included in the STL
#define M_PI
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
STL namespace.
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
Real nu
Definition: sabr.cpp:200