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 
21 #include <ql/math/solvers1d/brent.hpp>
22 #include <ql/math/functional.hpp>
23 #include <ql/math/distributions/chisquaredistribution.hpp>
24 #include <ql/math/distributions/gammadistribution.hpp>
25 #include <ql/math/distributions/normaldistribution.hpp>
26 
27 namespace QuantLib {
28 
operator ()(Real x) const29     Real CumulativeChiSquareDistribution::operator()(Real x) const {
30         return CumulativeGammaDistribution(0.5*df_)(0.5*x);
31     }
32 
operator ()(Real x) const33     Real NonCentralCumulativeChiSquareDistribution::operator()(Real x) const {
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 
operator ()(Real x) const93     Real NonCentralCumulativeChiSquareSankaranApprox::operator()(Real x) const {
94 
95         const Real h = 1-2*(df_+ncp_)*(df_+3*ncp_)/(3*square<Real>()(df_+2*ncp_));
96         const Real p = (df_+2*ncp_)/square<Real>()(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 
102         return CumulativeNormalDistribution()(u);
103     }
104 
105     InverseNonCentralCumulativeChiSquareDistribution::
InverseNonCentralCumulativeChiSquareDistribution(Real df,Real ncp,Size maxEvaluations,Real accuracy)106       InverseNonCentralCumulativeChiSquareDistribution(Real df, Real ncp,
107                                              Size maxEvaluations,
108                                              Real accuracy)
109     : nonCentralDist_(df, ncp),
110       guess_(df+ncp),
111       maxEvaluations_(maxEvaluations),
112       accuracy_(accuracy) {
113     }
114 
operator ()(Real x) const115     Real InverseNonCentralCumulativeChiSquareDistribution::operator()(Real x) const {
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(compose(subtract<Real>(x),
129                                     nonCentralDist_),
130                             accuracy_, 0.75*upper,
131                             (evaluations == maxEvaluations_)? 0.0: 0.5*upper,
132                             upper);
133     }
134 }
135