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#ifndef convolved_student_t_hpp
21#define convolved_student_t_hpp
22
23#include <ql/types.hpp>
24#include <vector>
25#include <numeric>
26#include <functional>
27
28namespace QuantLib {
29
30 /*! \brief Cumulative (generalized) BehrensFisher distribution.
31
32 Exact analitical computation of the cumulative probability distribution of
33 the linear combination of an arbitrary number (not just two) of T random
34 variables of odd integer order. Adapted from the algorithm in:\par
35 V. Witkovsky, Journal of Statistical Planning and Inference 94
36 (2001) 1-13\par
37 see also:\par
38 On the distribution of a linear combination of t-distributed
39 variables; Glenn Alan Walker, Ph.D.thesis University of Florida 1977\par
40 'Convolutions of the T Distribution'; S. Nadarajah, D. K. Dey in
41 Computers and Mathematics with Applications 49 (2005) 715-721\par
42 The last reference provides direct expressions for some of the densities
43 when the linear combination of only two Ts is just an addition. It can be
44 used for testing the results here.\par
45 Another available test on this algorithm stems from the realization that a
46 linear convex (\f$ \sum a_i=1\f$) combination of Ts of order one is stable
47 in the distribution sense (but this result is often of no practical use
48 because of its non-finite variance).\par
49 This implementation is for two or more T variables in the linear
50 combination albeit these must be of odd order. The case of exactly two T of
51 odd order is known to be a finite mixture of Ts but that result is not used
52 here. On this line see 'Linearization coefficients of Bessel polynomials'
53 C.Berg, C.Vignat; February 2008; arXiv:math/0506458
54
55 \todo Implement the series expansion solution for the addition of
56 two Ts of even order described in: 'On the density of the sum of two
57 independent Student t-random vectors' C.Berg, C.Vignat; June 2009;
58 eprint arXiv:0906.3037
59 */
60 class CumulativeBehrensFisher { // ODD orders only by now, rename?
61 public:
62 /*! \deprecated Use `auto` or `decltype` instead.
63 Deprecated in version 1.29.
64 */
65 QL_DEPRECATED
66 typedef Probability result_type;
67
68 /*! \deprecated Use `auto` or `decltype` instead.
69 Deprecated in version 1.29.
70 */
71 QL_DEPRECATED
72 typedef Real argument_type;
73 /*!
74 @param degreesFreedom Degrees of freedom of the Ts convolved. The
75 algorithm is limited to odd orders only.
76 @param factors Factors in the linear combination of the Ts.
77 */
78 CumulativeBehrensFisher(
79 const std::vector<Integer>& degreesFreedom = std::vector<Integer>(),
80 const std::vector<Real>& factors = std::vector<Real>());
81
82 //! Degrees of freedom of the Ts involved in the convolution.
83 const std::vector<Integer>& degreeFreedom() const {
84 return degreesFreedom_;
85 }
86 //! Factors in the linear combination.
87 const std::vector<Real>& factors() const {
88 return factors_;
89 }
90 private:
91 /*! \brief Student t characteristic polynomials.
92
93 Generates the polynomial coefficients defining the characteristic
94 function of a T distribution \f$T_\nu\f$ of odd order; \f$\nu=2n+1\f$.
95 In general the characteristic function is given by:
96 \f[
97 \phi_{\nu}(t) = \varphi_{n}(t) \exp{-\nu^{1/2}|t|} ;\,where\,\nu = 2n+1
98 \f]
99 where \f$ \varphi \f$ are polynomials that are computed recursively.
100
101 The convolved characteristic function is the product of the two previous
102 characteristic functions and the problem is then the convolution (a
103 product) of two polynomials.
104
105 @param n Natural number defining the order of the T for which
106 the characteristic function is to be computed. The order of the
107 T is then \f$ \nu=2n+1 \f$
108 */
109 // move outside of the class, as a separate problem?
110 std::vector<Real> polynCharactT(Natural n) const;
111
112 std::vector<Real> convolveVectorPolynomials(
113 const std::vector<Real>& v1,
114 const std::vector<Real>& v2) const ;
115 public:
116 /*! \brief Returns the cumulative probability of the resulting
117 distribution.\par
118 To obtain the cumulative probability the Gil-Pelaez theorem
119 is applied:\par
120 First compute the characteristic function of the linear combination
121 variable by multiplying the individual characteristic functions.
122 Then transform back integrating the characteristic function
123 according to the GP theorem; this is done here analytically feeding
124 in the expression of the total characteristic
125 function this:
126 \f[ \int_0^{\infty}x^n e^{-ax}sin(bx)dx =
127 (-1)^n \Gamma(n+1) \frac{sin((n+1)arctg2(-b/a))}
128 {(\sqrt{a^2+b^2})^{n+1}}; for\,a>0,\,b>0
129 \f]
130 and for the first term I use:
131 \f[
132 \int_0^{\infty} \frac{e^{-ax}sin(bx)}{x} dx = arctg2(b/a)
133 \f]
134 The GP complex integration is simplified thanks to the symetry of
135 the distribution.
136 */
137 Probability operator()(Real x) const;
138
139 /*! \brief Returns the probability density of the resulting
140 distribution.\par
141 Similarly to the cumulative probability, Gil-Pelaez theorem is
142 applied, the integration is similar.
143
144 \todo Implement in a separate class? given the name of this class..
145 */
146 Probability density(Real x) const;
147
148 private:
149 mutable std::vector<Integer> degreesFreedom_;
150 mutable std::vector<Real> factors_;
151
152 mutable std::vector<std::vector<Real> > polynCharFnc_;
153 mutable std::vector<Real> polyConvolved_;
154
155 // cached factor in the exponential of the characteristic function
156 mutable Real a_ = 0., a2_;
157 };
158
159
160
161 /*! \brief Inverse of the cumulative of the convolution of odd-T
162 distributions
163
164 Finds the inverse through a root solver. To find limits for the solver
165 domain use is made of the property that the convolved distribution is
166 bounded above by the normalized gaussian. If the coeffiecient in the linear
167 combination add up to a number below one the T of order one can be used as
168 a limit below but in general this is not necessarily the case and a constant
169 is used.
170 Also the fact that the combination is symmetric is used.
171 */
172 class InverseCumulativeBehrensFisher {
173 public:
174 /*! \deprecated Use `auto` or `decltype` instead.
175 Deprecated in version 1.29.
176 */
177 QL_DEPRECATED
178 typedef Real result_type;
179
180 /*! \deprecated Use `auto` or `decltype` instead.
181 Deprecated in version 1.29.
182 */
183 QL_DEPRECATED
184 typedef Probability argument_type;
185 /*!
186 @param degreesFreedom Degrees of freedom of the Ts convolved. The
187 algorithm is limited to odd orders only.
188 @param factors Factors in the linear combination of the Ts.
189 @param accuracy The accuracy of the root-solving process.
190 */
191 InverseCumulativeBehrensFisher(
192 const std::vector<Integer>& degreesFreedom = std::vector<Integer>(),
193 const std::vector<Real>& factors = std::vector<Real>(),
194 Real accuracy = 1.e-6);
195 //! Returns the cumulative inverse value.
196 Real operator()(Probability q) const;
197
198 private:
199 mutable Real normSqr_, accuracy_;
200 mutable CumulativeBehrensFisher distrib_;
201 };
202
203}
204
205#endif
206

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