QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
chisquaredistribution.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) 2002, 2003 Sadruddin Rejeb
5 Copyright (C) 2007 Klaus Spanderen
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
26
27namespace QuantLib {
28
30 return CumulativeGammaDistribution(0.5*df_)(0.5*x);
31 }
32
34 if (x <= 0.0)
35 return 0.0;
36
37 const Real errmax = 1e-12;
38 const Size itrmax = 10000;
39 Real lam = 0.5*ncp_;
40
41 Real u = std::exp(-lam);
42 Real v = u;
43 Real x2 = 0.5*x;
44 Real f2 = 0.5*df_;
45 Real f_x_2n = df_ - x;
46
47 Real t = 0.0;
48 if (f2*QL_EPSILON > 0.125 &&
49 std::fabs(x2-f2) < std::sqrt(QL_EPSILON)*f2) {
50 t = std::exp((1 - t) *
51 (2 - t/(f2+1)))/std::sqrt(2.0*M_PI*(f2 + 1.0));
52 }
53 else {
54 t = std::exp(f2*std::log(x2) - x2 -
55 GammaFunction().logValue(f2 + 1));
56 }
57
58 Real ans = v*t;
59
60 bool flag = false;
61 Size n = 1;
62 Real f_2n = df_ + 2.0;
63 f_x_2n += 2.0;
64
65 Real bound;
66 for (;;) {
67 if (f_x_2n > 0) {
68 flag = true;
69 goto L10;
70 }
71 for (;;) {
72 u *= lam / n;
73 v += u;
74 t *= x / f_2n;
75 ans += v*t;
76 n++;
77 f_2n += 2.0;
78 f_x_2n += 2.0;
79 if (!flag && n <= itrmax)
80 break;
81 L10:
82 bound = t * x / f_x_2n;
83 if (bound <= errmax || n > itrmax)
84 goto L_End;
85 }
86 }
87 L_End:
88 if (bound > errmax) QL_FAIL("didn't converge");
89 return (ans);
90
91 }
92
94
95 const Real h = 1-2*(df_+ncp_)*(df_+3*ncp_)/(3*squared(df_+2*ncp_));
96 const Real p = (df_+2*ncp_)/squared(df_+ncp_);
97 const Real m = (h-1)*(1-3*h);
98
99 const Real u= (std::pow(x/(df_+ncp_), h) - (1 + h*p*(h-1-0.5*(2-h)*m*p)))/
100 (h*std::sqrt(2*p)*(1+0.5*m*p));
101
103 }
104
107 Size maxEvaluations,
108 Real accuracy)
109 : nonCentralDist_(df, ncp),
110 guess_(df+ncp),
111 maxEvaluations_(maxEvaluations),
112 accuracy_(accuracy) {
113 }
114
116
117 // first find the right side of the interval
118 Real upper = guess_;
119 Size evaluations = maxEvaluations_;
120 while (nonCentralDist_(upper) < x && evaluations > 0) {
121 upper*=2.0;
122 --evaluations;
123 }
124
125 // use a Brent solver for the rest
126 Brent solver;
127 solver.setMaxEvaluations(evaluations);
128 return solver.solve([&](Real y) { return nonCentralDist_(y) - x; },
129 accuracy_, 0.75*upper,
130 (evaluations == maxEvaluations_)? 0.0: Real(0.5*upper),
131 upper);
132 }
133}
Brent 1-D solver.
Chi-square (central and non-central) distributions.
Brent 1-D solver
Definition: brent.hpp:37
Cumulative normal distribution function.
Gamma function class.
InverseNonCentralCumulativeChiSquareDistribution(Real df, Real ncp, Size maxEvaluations=10, Real accuracy=1e-8)
NonCentralCumulativeChiSquareDistribution nonCentralDist_
void setMaxEvaluations(Size evaluations)
Definition: solver1d.hpp:238
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
const DefaultType & t
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Gamma distribution.
#define QL_EPSILON
Definition: qldefines.hpp:178
QL_REAL Real
real number
Definition: types.hpp:50
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
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< BlackVolTermStructure > v