1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4  Copyright (C) 2015 Johannes Göttker-Schnetmann
5  Copyright (C) 2015 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 
22 #include <ql/methods/finitedifferences/utilities/squarerootprocessrndcalculator.hpp>
23 
24 #include <boost/math/distributions/non_central_chi_squared.hpp>
25 
26 namespace QuantLib {
27 
SquareRootProcessRNDCalculator(Real v0,Real kappa,Real theta,Real sigma)28     SquareRootProcessRNDCalculator::SquareRootProcessRNDCalculator(
29         Real v0, Real kappa, Real theta, Real sigma)
30     : v0_(v0), kappa_(kappa), theta_(theta),
31       d_(4*kappa/(sigma*sigma)), df_(d_*theta) {    }
32 
33 
pdf(Real v,Time t) const34     Real SquareRootProcessRNDCalculator::pdf(Real v, Time t) const {
35         const Real e   = std::exp(-kappa_*t);
36         const Real k   = d_/(1-e);
37         const Real ncp = k*v0_*e;
38 
39         const boost::math::non_central_chi_squared_distribution<Real>
40             dist(df_, ncp);
41 
42         return boost::math::pdf(dist, v*k) * k;
43     }
44 
cdf(Real v,Time t) const45     Real SquareRootProcessRNDCalculator::cdf(Real v, Time t) const {
46         const Real e   = std::exp(-kappa_*t);
47         const Real k   = d_/(1-e);
48         const Real ncp = k*v0_*e;
49 
50         const boost::math::non_central_chi_squared_distribution<Real>
51             dist(df_, ncp);
52 
53         return boost::math::cdf(dist, v*k);
54     }
55 
invcdf(Real q,Time t) const56     Real SquareRootProcessRNDCalculator::invcdf(Real q, Time t) const {
57         const Real e   = std::exp(-kappa_*t);
58         const Real k   = d_/(1-e);
59         const Real ncp = k*v0_*e;
60 
61         const boost::math::non_central_chi_squared_distribution<Real>
62             dist(df_, ncp);
63 
64         return boost::math::quantile(dist, q) / k;
65     }
66 
stationary_pdf(Real v) const67     Real SquareRootProcessRNDCalculator::stationary_pdf(Real v) const {
68         const Real alpha = 0.5*df_;
69         const Real beta = alpha/theta_;
70 
71         return std::pow(beta, alpha)*std::pow(v, alpha-1)
72                 *std::exp(-beta*v-boost::math::lgamma(alpha));
73     }
74 
stationary_cdf(Real v) const75     Real SquareRootProcessRNDCalculator::stationary_cdf(Real v) const {
76         const Real alpha = 0.5*df_;
77         const Real beta = alpha/theta_;
78 
79         return boost::math::gamma_p(alpha, beta*v);
80     }
81 
stationary_invcdf(Real q) const82     Real SquareRootProcessRNDCalculator::stationary_invcdf(Real q) const {
83         const Real alpha = 0.5*df_;
84         const Real beta = alpha/theta_;
85 
86         return boost::math::gamma_p_inv(alpha, q)/beta;
87     }
88 }
89