QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
hestonrndcalculator.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) 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
29#include <complex>
30#include <utility>
31
32namespace QuantLib {
33
34namespace {
35 struct HestonParams {
37 };
38
39 HestonParams getHestonParams(
40 const ext::shared_ptr<HestonProcess>& process) {
41 const HestonParams p = { process->v0(), process->kappa(),
42 process->theta(), process->sigma(),
43 process->rho() };
44 return p;
45 }
46
47 std::complex<Real> gamma(const HestonParams& p, Real p_x) {
48 return std::complex<Real>(p.kappa, p.rho*p.sigma*p_x);
49 }
50
51 std::complex<Real> omega(const HestonParams& p, Real p_x) {
52 const std::complex<Real> g = gamma(p, p_x);
53 return std::sqrt(g*g
54 + p.sigma*p.sigma*std::complex<Real>(p_x*p_x, -p_x));
55 }
56
57 class CpxPv_Helper {
58 public:
59 CpxPv_Helper(const HestonParams& p, Real x, Time t)
60 : p_(p), t_(t), x_(x),
61 c_inf_(std::min(10.0, std::max(0.0001,
62 std::sqrt(1.0-squared(p_.rho))/p_.sigma))
63 *(p_.v0 + p_.kappa*p_.theta*t)) {}
64
65 Real operator()(Real x) const {
66 return std::real(transformPhi(x));
67 }
68
69 Real p0(Real p_x) const {
70 if (p_x < QL_EPSILON) {
71 return 0.0;
72 }
73
74 const Real u_x = std::max(QL_EPSILON, -std::log(p_x)/c_inf_);
75 return std::real(phi(u_x)
76 /((p_x*c_inf_)*std::complex<Real>(0.0, u_x)));
77 }
78
79 private:
80 std::complex<Real> transformPhi(Real x) const {
81 if (x < QL_EPSILON) {
82 return std::complex<Real>(0.0, 0.0);
83 }
84
85 const Real u_x = -std::log(x)/c_inf_;
86 return phi(u_x)/(x*c_inf_);
87 }
88
89 std::complex<Real> phi(Real p_x) const {
90 const Real sigma2 = p_.sigma*p_.sigma;
91 const std::complex<Real> g = gamma(p_, p_x);
92 const std::complex<Real> o = omega(p_, p_x);
93 const std::complex<Real> gamma = (g-o)/(g+o);
94
95 return 2.0*std::exp(std::complex<Real>(0.0, p_x*x_)
96 - p_.v0*std::complex<Real>(p_x*p_x, -p_x)
97 /(g+o*(1.0+std::exp(-o*t_))/(1.0-std::exp(-o*t_)))
98 +p_.kappa*p_.theta/sigma2*(
99 (g-o)*t_ - 2.0*std::log((1.0-gamma*std::exp(-o*t_))
100 /(1.0-gamma))));
101 }
102
103 const HestonParams p_;
104 const Time t_;
105 const Real x_, c_inf_;
106 };
107 }
108
109
110 HestonRNDCalculator::HestonRNDCalculator(ext::shared_ptr<HestonProcess> hestonProcess,
111 Real integrationEps,
112 Size maxIntegrationIterations)
113 : hestonProcess_(std::move(hestonProcess)), x0_(std::log(hestonProcess_->s0()->value())),
114 integrationEps_(integrationEps), maxIntegrationIterations_(maxIntegrationIterations) {}
115
117 const DiscountFactor dr = hestonProcess_->riskFreeRate()->discount(t);
118 const DiscountFactor dq = hestonProcess_->dividendYield()->discount(t);
119
120 return x - x0_ + std::log(dr/dq);
121 }
122
126 CpxPv_Helper(getHestonParams(hestonProcess_), x_t(x, t), t),
127 0.0, 1.0)/M_TWOPI;
128 }
129
131 CpxPv_Helper helper(getHestonParams(hestonProcess_), x_t(x, t), t);
132
134 [&](Real p_x){ return helper.p0(p_x); },
135 0.0, 1.0)/M_TWOPI + 0.5;
136 }
137
139 const Real v0 = hestonProcess_->v0();
140 const Real kappa = hestonProcess_->kappa();
141 const Real theta = hestonProcess_->theta();
142
143 const Volatility expVol
144 = std::sqrt(theta + (v0-theta)*(1-std::exp(-kappa*t))/(t*kappa));
145
146 const ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
147 ext::make_shared<BlackScholesMertonProcess>(
148 hestonProcess_->s0(),
149 hestonProcess_->dividendYield(),
150 hestonProcess_->riskFreeRate(),
152 ext::make_shared<BlackConstantVol>(
153 hestonProcess_->riskFreeRate()->referenceDate(),
154 NullCalendar(),
155 expVol,
156 hestonProcess_->riskFreeRate()->dayCounter()))));
157
158 const Real guess = BSMRNDCalculator(bsmProcess).invcdf(p, t);
159
162 .inverseCDF(p, t);
163 }
164}
Black constant volatility, no time dependence, no strike dependence.
Black-Scholes processes.
risk neutral terminal density calculator for the Black-Scholes-Merton model with constant volatility
Real invcdf(Real q, Time t) const override
Integral of a one-dimensional function.
Shared handle to an observable.
Definition: handle.hpp:41
HestonRNDCalculator(ext::shared_ptr< HestonProcess > hestonProcess, Real integrationEps=1e-6, Size maxIntegrationIterations=10000UL)
Real x_t(Real x, Time t) const
Real invcdf(Real q, Time t) const override
Real pdf(Real x, Time t) const override
Real cdf(Real x, Time t) const override
const ext::shared_ptr< HestonProcess > hestonProcess_
Calendar for reproducing theoretical calculations.
const DefaultType & t
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
#define QL_EPSILON
Definition: qldefines.hpp:178
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real DiscountFactor
discount factor between dates
Definition: types.hpp:66
Real Volatility
volatility
Definition: types.hpp:78
std::size_t Size
size of a container
Definition: types.hpp:58
Heston stochastic process.
Real kappa
Real theta
const Real c_inf_
Real v0
const HestonParams p_
Real rho
Real sigma
risk neutral terminal density calculator for the Heston stochastic volatility model
functionals and combinators not included in the STL
#define M_TWOPI
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
STL namespace.
Calendar for reproducing theoretical calculations.
std::uint64_t x_