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