31#include <boost/math/distributions/non_central_chi_squared.hpp>
47 riskFreeRate_(
std::move(riskFreeRate)), dividendYield_(
std::move(dividendYield)),
67 return {
s0_->value(),
v0_ };
71 const Real vol = (x[1] > 0.0) ? std::sqrt(x[1])
92 const Real vol = (x[1] > 0.0) ? std::sqrt(x[1])
99 tmp[0][0] = vol; tmp[0][1] = 0.0;
100 tmp[1][0] =
rho_*sigma2; tmp[1][1] = sqrhov*sigma2;
105 const Array& dx)
const {
107 x0[0] * std::exp(dx[0]),
124 const std::complex<Real>& a,
131 const std::complex<Real> ga = std::sqrt(
132 kappa*
kappa - 2*sigma2*a*std::complex<Real>(0.0, 1.0));
136 const std::complex<Real> z
137 = ga*std::exp(-0.5*ga*dt)/(1.0-std::exp(-ga*dt));
138 const std::complex<Real> log_z
139 = -0.5*ga*dt + std::log(ga/(1.0-std::exp(-ga*dt)));
141 const std::complex<Real>
alpha
142 = 4.0*ga*std::exp(-0.5*ga*dt)/(sigma2*(1.0-std::exp(-ga*dt)));
144 /(sigma2*(1.0-std::exp(-
kappa*dt)));
146 return ga*std::exp(-0.5*(ga-
kappa)*dt)*(1-std::exp(-
kappa*dt))
147 / (
kappa*(1.0-std::exp(-ga*dt)))
148 *std::exp((nu_0+nu_t)/sigma2 * (
150 - ga*(1.0+std::exp(-ga*dt))/(1.0-std::exp(-ga*dt))))
151 *std::exp(
nu*log_z)/std::pow(z,
nu)
154 nu, std::sqrt(nu_0*nu_t)*
alpha)
156 nu, std::sqrt(nu_0*nu_t)*
beta)
161 Real ch(
const HestonProcess& process,
163 return M_2_PI*std::sin(u*x)/u
164 * Phi(process, u, nu_0, nu_t, dt).real();
167 Real ph(
const HestonProcess& process,
169 return M_2_PI*std::cos(u*x)*Phi(process, u, nu_0, nu_t, dt).real();
172 Real int_ph(
const HestonProcess& process,
174 static const GaussLaguerreIntegration gaussLaguerreIntegration(128);
176 const Real rho = process.rho();
179 const Real x0 = std::log(process.s0()->value());
181 return gaussLaguerreIntegration(
182 [&](
Real u){
return ph(process,
y, u, nu_0, nu_t,
t); })
191 for (
Integer i=m-1; i >= 0; --i) {
192 n = (
n+nominator[i])*x;
193 d = (
d+denominator[i])*x;
203 { -4.54393409816329991e-2,1.15457225751016682e-3,
204 -1.41018536821330254e-5,9.43280809438713025e-8,
205 -3.53201978997168357e-10,7.08240282274875911e-13,
206 -6.05338212010422477e-16 };
208 { 1.01162145739225565e-2,4.99175116169755106e-5,
209 1.55654986308745614e-7,3.28067571055789734e-10,
210 4.5049097575386581e-13,3.21107051193712168e-16,
213 return x*pade(x*x,
n,
d,
sizeof(
n)/
sizeof(
Real));
216 const Real y = 1/(x*x);
218 { 7.44437068161936700618e2,1.96396372895146869801e5,
219 2.37750310125431834034e7,1.43073403821274636888e9,
220 4.33736238870432522765e10,6.40533830574022022911e11,
221 4.20968180571076940208e12,1.00795182980368574617e13,
222 4.94816688199951963482e12,-4.94701168645415959931e11 };
224 { 7.46437068161927678031e2,1.97865247031583951450e5,
225 2.41535670165126845144e7,1.47478952192985464958e9,
226 4.58595115847765779830e10,7.08501308149515401563e11,
227 5.06084464593475076774e12,1.43468549171581016479e13,
228 1.11535493509914254097e13, 0.0 };
229 const Real f = pade(
y, fn, fd, 10)/x;
232 { 8.1359520115168615e2,2.35239181626478200e5,
233 3.12557570795778731e7,2.06297595146763354e9,
234 6.83052205423625007e10,1.09049528450362786e12,
235 7.57664583257834349e12,1.81004487464664575e13,
236 6.43291613143049485e12,-1.36517137670871689e12 };
238 { 8.19595201151451564e2,2.40036752835578777e5,
239 3.26026661647090822e7,2.23355543278099360e9,
240 7.87465017341829930e10,1.39866710696414565e12,
241 1.17164723371736605e13,4.01839087307656620e13,
242 3.99653257887490811e13, 0.0};
243 const Real g =
y*pade(
y, gn, gd, 10);
245 return M_PI_2 -
f*std::cos(x)-
g*std::sin(x);
249 Real cornishFisherEps(
const HestonProcess& process,
254 const Real p2 = Phi(process, std::complex<Real>(0, -2*
d),
255 nu_0, nu_t, dt).real();
256 const Real p1 = Phi(process, std::complex<Real>(0, -
d),
257 nu_0, nu_t, dt).real();
258 const Real p0 = Phi(process, std::complex<Real>(0, 0),
259 nu_0, nu_t, dt).real();
260 const Real pm1= Phi(process, std::complex<Real>(0,
d),
261 nu_0, nu_t, dt).real();
262 const Real pm2= Phi(process, std::complex<Real>(0, 2*
d),
263 nu_0, nu_t, dt).real();
265 const Real avg = (pm2-8*pm1+8*p1-p2)/(12*
d);
266 const Real m2 = (-pm2+16*pm1-30*p0+16*p1-p2)/(12*
d*
d);
267 const Real var = m2 - avg*avg;
268 const Real stdDev = std::sqrt(var);
270 const Real m3 = (-0.5*pm2 + pm1 - p1 + 0.5*p2)/(
d*
d*
d);
272 = (m3 - 3*var*avg - avg*avg*avg) / (var*stdDev);
274 const Real m4 = (pm2 - 4*pm1 + 6*p0 - 4*p1 + p2)/(
d*
d*
d*
d);
276 = (m4 - 4*m3*avg + 6*m2*avg*avg - 3*avg*avg*avg*avg)
281 const Real q = InverseCumulativeNormal()(1-eps);
282 const Real w =
q + (
q*
q-1)/6*skew + (
q*
q*
q-3*
q)/24*(kurt-3)
283 - (2*
q*
q*
q-5*
q)/36*skew*skew;
285 return avg + w*stdDev;
288 Real cdf_nu_ds(
const HestonProcess& process,
291 const Real eps = 1e-4;
292 const Real u_eps = std::min(100.0,
293 std::max(0.1, cornishFisherEps(process, nu_0, nu_t, dt, eps)));
295 switch (discretization) {
298 static const GaussLaguerreIntegration
299 gaussLaguerreIntegration(128);
302 Real upper = u_eps/2.0;
303 while (std::abs(Phi(process,upper,nu_0,nu_t,dt)/upper)
307 ? std::max(0.0, std::min(1.0,
308 gaussLaguerreIntegration(
309 [&](
Real u){
return ch(process, x, u, nu_0, nu_t, dt); })))
315 Real upper = u_eps/2.0;
316 while (std::abs(Phi(process, upper,nu_0,nu_t,dt)/upper)
320 ? std::max(0.0, std::min(1.0,
321 GaussLobattoIntegral(Null<Size>(), eps)(
322 [&](
Real xi){
return ch(process, x, xi, nu_0, nu_t, dt); },
332 std::complex<Real>
f;
337 const Real si_n =
Si(x*(u+0.5*h));
339 f = Phi(process, u, nu_0, nu_t, dt);
343 while (
M_2_PI*std::abs(
f)/j > eps);
348 QL_FAIL(
"unknown integration method");
357 return cdf_nu_ds(process, x, nu_0, nu_t, dt, discretization) - x0;
366 const Real x0 = std::log(
s0()->value());
369 while (df > 0.0 ||
f > 0.1*eps) {
374 * ( -0.5/(upper*std::sqrt(upper))*std::exp(f2)
375 + 1/std::sqrt(upper)*std::exp(f2)*(-0.5/(1-
rho_*
rho_))
376 *(-1/(upper*upper)*f1*f1
383 upper = 2.0*cornishFisherEps(*
this,
v0_,
v,
t,1e-3);
386 [&](
Real xi){
return int_ph(*
this, a, x, xi,
v0_,
v,
t); },
389 boost::math::non_central_chi_squared_distribution<Real>(
399 Real vol, vol2, mu,
nu, dy;
401 const Real sdt = std::sqrt(dt);
411 vol = (x0[1] > 0.0) ? std::sqrt(x0[1]) :
Real(0.0);
418 retVal[0] = x0[0] * std::exp(mu*dt+vol*dw[0]*sdt);
419 retVal[1] = x0[1] +
nu*dt + vol2*sdt*(
rho_*dw[0] + sqrhov*dw[1]);
422 vol = (x0[1] > 0.0) ? std::sqrt(x0[1]) :
Real(0.0);
429 retVal[0] = x0[0] * std::exp(mu*dt+vol*dw[0]*sdt);
430 retVal[1] = x0[1] +
nu*dt + vol2*sdt*(
rho_*dw[0] + sqrhov*dw[1]);
433 vol = std::sqrt(std::fabs(x0[1]));
440 retVal[0] = x0[0]*std::exp(mu*dt+vol*dw[0]*sdt);
442 +
nu*dt + vol2*sdt*(
rho_*dw[0] + sqrhov*dw[1]);
450 vol = (x0[1] > 0.0) ? std::sqrt(x0[1]) :
Real(0.0);
457 *(
theta_-vol*vol)) * dt + vol*sqrhov*dw[0]*sdt;
459 retVal[0] = x0[0]*std::exp(dy +
rho_/
sigma_*(retVal[1]-x0[1]));
472 const Real psi = s2/(m*m);
481 const Real A = k2+0.5*k4;
484 const Real b2 = 2/psi-1+std::sqrt(2/psi*(2/psi-1));
485 const Real b = std::sqrt(b2);
486 const Real a = m/(1+b2);
491 k0 = -A*b2*a/(1-2*A*a)+0.5*std::log(1-2*A*a)
494 retVal[1] = a*(
b+dw[1])*(
b+dw[1]);
497 const Real p = (psi-1)/(psi+1);
505 k0 = -std::log(p+
beta*(1-p)/(
beta-A))-(k1+0.5*k3)*x0[1];
507 retVal[1] = ((u <= p) ?
Real(0.0) : std::log((1-p)/(1-u))/
beta);
513 retVal[0] = x0[0]*std::exp(mu*dt + k0 + k1*x0[1] + k2*retVal[1]
514 +std::sqrt(k3*x0[1]+k4*retVal[1])*dw[0]);
521 const Real nu_0 = x0[1];
536 - 0.5*vds +
rho_*vdw;
539 const Real s = x0[0]*std::exp(mu + sig*dw[0]);
546 QL_FAIL(
"unknown discretization schema");
Chi-square (central and non-central) distributions.
1-D array used in linear algebra.
Cumulative normal distribution function.
Euler discretization for stochastic processes.
Shared handle to an observable.
Square-root stochastic-volatility Heston process.
HestonProcess(Handle< YieldTermStructure > riskFreeRate, Handle< YieldTermStructure > dividendYield, Handle< Quote > s0, Real v0, Real kappa, Real theta, Real sigma, Real rho, Discretization d=QuadraticExponentialMartingale)
Array drift(Time t, const Array &x) const override
returns the drift part of the equation, i.e.,
Size size() const override
returns the number of dimensions of the stochastic process
Array evolve(Time t0, const Array &x0, Time dt, const Array &dw) const override
Real pdf(Real x, Real v, Time t, Real eps=1e-3) const
@ BroadieKayaExactSchemeLobatto
@ NonCentralChiSquareVariance
@ BroadieKayaExactSchemeLaguerre
@ QuadraticExponentialMartingale
@ BroadieKayaExactSchemeTrapezoidal
Real varianceDistribution(Real v, Real dw, Time dt) const
Matrix diffusion(Time t, const Array &x) const override
returns the diffusion part of the equation, i.e.
const Handle< YieldTermStructure > & dividendYield() const
Time time(const Date &) const override
Size factors() const override
returns the number of independent factors of the process
Discretization discretization_
Array initialValues() const override
returns the initial values of the state variables
Handle< YieldTermStructure > dividendYield_
Array apply(const Array &x0, const Array &dx) const override
const Handle< Quote > & s0() const
const Handle< YieldTermStructure > & riskFreeRate() const
Handle< YieldTermStructure > riskFreeRate_
Matrix used in linear algebra.
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Integral of a one-dimensional function.
Real solve(const F &f, Real accuracy, Real guess, Real step) const
discretization of a stochastic process over a given time interval
multi-dimensional stochastic process class.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Euler discretization for stochastic processes.
ext::function< Real(Real)> b
Integral of a 1-dimensional function using the Gauss quadratures.
integral of a one-dimensional function using the adaptive Gauss-Lobatto integral
Real Time
continuous quantity with 1-year units
Real Volatility
volatility
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
Heston stochastic process.
functionals and combinators not included in the STL
modified Bessel functions of first and second kind
Real modifiedBesselFunction_i(Real nu, Real x)
Real cdf_nu_ds_minus_x(const HestonProcess &process, Real x, Real nu_0, Real nu_t, Time dt, HestonProcess::Discretization discretization, Real x0)
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< BlackVolTermStructure > v
Integral of a one-dimensional function using segment algorithm.