QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
hestonprocess.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) 2005, 2007, 2009, 2014 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
31#include <boost/math/distributions/non_central_chi_squared.hpp>
32#include <complex>
33#include <utility>
34
35namespace QuantLib {
36
38 Handle<YieldTermStructure> dividendYield,
40 Real v0,
41 Real kappa,
42 Real theta,
43 Real sigma,
44 Real rho,
47 riskFreeRate_(std::move(riskFreeRate)), dividendYield_(std::move(dividendYield)),
48 s0_(std::move(s0)), v0_(v0), kappa_(kappa), theta_(theta), sigma_(sigma), rho_(rho),
49 discretization_(d) {
50
54 }
55
57 return 2;
58 }
59
64 }
65
67 return { s0_->value(), v0_ };
68 }
69
71 const Real vol = (x[1] > 0.0) ? std::sqrt(x[1])
72 : (discretization_ == Reflection) ? Real(- std::sqrt(-x[1]))
73 : 0.0;
74
75 return {
76 riskFreeRate_->forwardRate(t, t, Continuous).rate()
77 - dividendYield_->forwardRate(t, t, Continuous).rate()
78 - 0.5 * vol * vol,
79 kappa_* (theta_-((discretization_==PartialTruncation) ? x[1] : vol*vol))
80 };
81 }
82
84 /* the correlation matrix is
85 | 1 rho |
86 | rho 1 |
87 whose square root (which is used here) is
88 | 1 0 |
89 | rho sqrt(1-rho^2) |
90 */
91 Matrix tmp(2,2);
92 const Real vol = (x[1] > 0.0) ? std::sqrt(x[1])
93 : (discretization_ == Reflection) ? Real(-std::sqrt(-x[1]))
94 : 1e-8; // set vol to (almost) zero but still
95 // expose some correlation information
96 const Real sigma2 = sigma_ * vol;
97 const Real sqrhov = std::sqrt(1.0 - rho_*rho_);
98
99 tmp[0][0] = vol; tmp[0][1] = 0.0;
100 tmp[1][0] = rho_*sigma2; tmp[1][1] = sqrhov*sigma2;
101 return tmp;
102 }
103
105 const Array& dx) const {
106 return {
107 x0[0] * std::exp(dx[0]),
108 x0[1] + dx[1]
109 };
110 }
111
112 namespace {
113 // This is the continuous version of a characteristic function
114 // for the exact sampling of the Heston process, s. page 8, formula 13,
115 // M. Broadie, O. Kaya, Exact Simulation of Stochastic Volatility and
116 // other Affine Jump Diffusion Processes
117 // http://finmath.stanford.edu/seminars/documents/Broadie.pdf
118 //
119 // This version does not need a branch correction procedure.
120 // For details please see:
121 // Roger Lord, "Efficient Pricing Algorithms for exotic Derivatives",
122 // http://repub.eur.nl/pub/13917/LordR-Thesis.pdf
123 std::complex<Real> Phi(const HestonProcess& process,
124 const std::complex<Real>& a,
125 Real nu_0, Real nu_t, Time dt) {
126 const Real theta = process.theta();
127 const Real kappa = process.kappa();
128 const Real sigma = process.sigma();
129
130 const Volatility sigma2 = sigma*sigma;
131 const std::complex<Real> ga = std::sqrt(
132 kappa*kappa - 2*sigma2*a*std::complex<Real>(0.0, 1.0));
133 const Real d = 4*theta*kappa/sigma2;
134
135 const Real nu = 0.5*d-1;
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)));
140
141 const std::complex<Real> alpha
142 = 4.0*ga*std::exp(-0.5*ga*dt)/(sigma2*(1.0-std::exp(-ga*dt)));
143 const std::complex<Real> beta = 4.0*kappa*std::exp(-0.5*kappa*dt)
144 /(sigma2*(1.0-std::exp(-kappa*dt)));
145
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 * (
149 kappa*(1.0+std::exp(-kappa*dt))/(1.0-std::exp(-kappa*dt))
150 - ga*(1.0+std::exp(-ga*dt))/(1.0-std::exp(-ga*dt))))
151 *std::exp(nu*log_z)/std::pow(z, nu)
152 *((nu_t > 1e-8)
154 nu, std::sqrt(nu_0*nu_t)*alpha)
156 nu, std::sqrt(nu_0*nu_t)*beta)
157 : std::pow(alpha/beta, nu)
158 );
159 }
160
161 Real ch(const HestonProcess& process,
162 Real x, Real u, Real nu_0, Real nu_t, Time dt) {
163 return M_2_PI*std::sin(u*x)/u
164 * Phi(process, u, nu_0, nu_t, dt).real();
165 }
166
167 Real ph(const HestonProcess& process,
168 Real x, Real u, Real nu_0, Real nu_t, Time dt) {
169 return M_2_PI*std::cos(u*x)*Phi(process, u, nu_0, nu_t, dt).real();
170 }
171
172 Real int_ph(const HestonProcess& process,
173 Real a, Real x, Real y, Real nu_0, Real nu_t, Time t) {
174 static const GaussLaguerreIntegration gaussLaguerreIntegration(128);
175
176 const Real rho = process.rho();
177 const Real kappa = process.kappa();
178 const Real sigma = process.sigma();
179 const Real x0 = std::log(process.s0()->value());
180
181 return gaussLaguerreIntegration(
182 [&](Real u){ return ph(process, y, u, nu_0, nu_t, t); })
183 / std::sqrt(2*M_PI*(1-rho*rho)*y)
184 * std::exp(-0.5*squared(x - x0 - a + y*(0.5-rho*kappa/sigma))
185 /(y*(1-rho*rho)));
186 }
187
188
189 Real pade(Real x, const Real* nominator, const Real* denominator, Size m) {
190 Real n=0.0, d=0.0;
191 for (Integer i=m-1; i >= 0; --i) {
192 n = (n+nominator[i])*x;
193 d = (d+denominator[i])*x;
194 }
195 return (1+n)/(1+d);
196 }
197
198 // For the definition of the Pade approximation please see e.g.
199 // http://wikipedia.org/wiki/Sine_integral#Sine_integral
200 Real Si(Real x) {
201 if (x <=4.0) {
202 const Real n[] =
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 };
207 const Real d[] =
208 { 1.01162145739225565e-2,4.99175116169755106e-5,
209 1.55654986308745614e-7,3.28067571055789734e-10,
210 4.5049097575386581e-13,3.21107051193712168e-16,
211 0.0 };
212
213 return x*pade(x*x, n, d, sizeof(n)/sizeof(Real));
214 }
215 else {
216 const Real y = 1/(x*x);
217 const Real fn[] =
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 };
223 const Real fd[] =
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;
230
231 const Real gn[] =
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 };
237 const Real gd[] =
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);
244
245 return M_PI_2 - f*std::cos(x)-g*std::sin(x);
246 }
247 }
248
249 Real cornishFisherEps(const HestonProcess& process,
250 Real nu_0, Real nu_t, Time dt, Real eps) {
251 // use moment generating function to get the
252 // first,second, third and fourth moment of the distribution
253 const Real d = 1e-2;
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();
264
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);
269
270 const Real m3 = (-0.5*pm2 + pm1 - p1 + 0.5*p2)/(d*d*d);
271 const Real skew
272 = (m3 - 3*var*avg - avg*avg*avg) / (var*stdDev);
273
274 const Real m4 = (pm2 - 4*pm1 + 6*p0 - 4*p1 + p2)/(d*d*d*d);
275 const Real kurt
276 = (m4 - 4*m3*avg + 6*m2*avg*avg - 3*avg*avg*avg*avg)
277 /(var*var);
278
279 // Cornish-Fisher relation to come up with an improved
280 // estimate of 1-F(u_\eps) < \eps
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;
284
285 return avg + w*stdDev;
286 }
287
288 Real cdf_nu_ds(const HestonProcess& process,
289 Real x, Real nu_0, Real nu_t, Time dt,
290 HestonProcess::Discretization discretization) {
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)));
294
295 switch (discretization) {
297 {
298 static const GaussLaguerreIntegration
299 gaussLaguerreIntegration(128);
300
301 // get the upper bound for the integration
302 Real upper = u_eps/2.0;
303 while (std::abs(Phi(process,upper,nu_0,nu_t,dt)/upper)
304 > eps) upper*=2.0;
305
306 return (x < 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); })))
310 : Real(1.0);
311 }
313 {
314 // get the upper bound for the integration
315 Real upper = u_eps/2.0;
316 while (std::abs(Phi(process, upper,nu_0,nu_t,dt)/upper)
317 > eps) upper*=2.0;
318
319 return (x < 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); },
323 QL_EPSILON, upper)))
324 : Real(1.0);
325 }
327 {
328 const Real h = 0.05;
329
330 Real si = Si(0.5*h*x);
331 Real s = M_2_PI*si;
332 std::complex<Real> f;
333 Size j = 0;
334 do {
335 ++j;
336 const Real u = h*j;
337 const Real si_n = Si(x*(u+0.5*h));
338
339 f = Phi(process, u, nu_0, nu_t, dt);
340 s+= M_2_PI*f.real()*(si_n-si);
341 si = si_n;
342 }
343 while (M_2_PI*std::abs(f)/j > eps);
344
345 return s;
346 }
347 default:
348 QL_FAIL("unknown integration method");
349 }
350 }
351 }
352
354 Real nu_t, Time dt,
355 HestonProcess::Discretization discretization,
356 Real x0) {
357 return cdf_nu_ds(process, x, nu_0, nu_t, dt, discretization) - x0;
358 }
359
361 const Real k = sigma_*sigma_*(1-std::exp(-kappa_*t))/(4*kappa_);
362 const Real a = std::log( dividendYield_->discount(t)
363 / riskFreeRate_->discount(t))
364 + rho_/sigma_*(v - v0_ - kappa_*theta_*t);
365
366 const Real x0 = std::log(s0()->value());
367 Real upper = std::max(0.1, -(x-x0-a)/(0.5-rho_*kappa_/sigma_)), f=0, df=1;
368
369 while (df > 0.0 || f > 0.1*eps) {
370 const Real f1 = x-x0-a+upper*(0.5-rho_*kappa_/sigma_);
371 const Real f2 = -0.5*f1*f1/(upper*(1-rho_*rho_));
372
373 df = 1/std::sqrt(2*M_PI*(1-rho_*rho_))
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
377 + 2/upper*f1*(0.5-rho_*kappa_/sigma_)));
378
379 f = std::exp(f2)/ std::sqrt(2*M_PI*(1-rho_*rho_)*upper);
380 upper*=1.5;
381 }
382
383 upper = 2.0*cornishFisherEps(*this, v0_, v, t,1e-3);
384
385 return SegmentIntegral(100)(
386 [&](Real xi){ return int_ph(*this, a, x, xi, v0_, v, t); },
387 QL_EPSILON, upper)
388 * boost::math::pdf(
389 boost::math::non_central_chi_squared_distribution<Real>(
391 4*kappa_*std::exp(-kappa_*t)
392 /((sigma_*sigma_)*(1-std::exp(-kappa_*t)))*v0_),
393 v/k) / k;
394 }
395
397 Time dt, const Array& dw) const {
398 Array retVal(2);
399 Real vol, vol2, mu, nu, dy;
400
401 const Real sdt = std::sqrt(dt);
402 const Real sqrhov = std::sqrt(1.0 - rho_*rho_);
403
404 switch (discretization_) {
405 // For the definition of PartialTruncation, FullTruncation
406 // and Reflection see Lord, R., R. Koekkoek and D. van Dijk (2006),
407 // "A Comparison of biased simulation schemes for
408 // stochastic volatility models",
409 // Working Paper, Tinbergen Institute
411 vol = (x0[1] > 0.0) ? std::sqrt(x0[1]) : Real(0.0);
412 vol2 = sigma_ * vol;
413 mu = riskFreeRate_->forwardRate(t0, t0+dt, Continuous).rate()
414 - dividendYield_->forwardRate(t0, t0+dt, Continuous).rate()
415 - 0.5 * vol * vol;
416 nu = kappa_*(theta_ - x0[1]);
417
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]);
420 break;
421 case FullTruncation:
422 vol = (x0[1] > 0.0) ? std::sqrt(x0[1]) : Real(0.0);
423 vol2 = sigma_ * vol;
424 mu = riskFreeRate_->forwardRate(t0, t0+dt, Continuous).rate()
425 - dividendYield_->forwardRate(t0, t0+dt, Continuous).rate()
426 - 0.5 * vol * vol;
427 nu = kappa_*(theta_ - vol*vol);
428
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]);
431 break;
432 case Reflection:
433 vol = std::sqrt(std::fabs(x0[1]));
434 vol2 = sigma_ * vol;
435 mu = riskFreeRate_->forwardRate(t0, t0+dt, Continuous).rate()
436 - dividendYield_->forwardRate(t0, t0+dt, Continuous).rate()
437 - 0.5 * vol*vol;
438 nu = kappa_*(theta_ - vol*vol);
439
440 retVal[0] = x0[0]*std::exp(mu*dt+vol*dw[0]*sdt);
441 retVal[1] = vol*vol
442 +nu*dt + vol2*sdt*(rho_*dw[0] + sqrhov*dw[1]);
443 break;
445 // use Alan Lewis trick to decorrelate the equity and the variance
446 // process by using y(t)=x(t)-\frac{rho}{sigma}\nu(t)
447 // and Ito's Lemma. Then use exact sampling for the variance
448 // process. For further details please read the Wilmott thread
449 // "QuantLib code is very high quality"
450 vol = (x0[1] > 0.0) ? std::sqrt(x0[1]) : Real(0.0);
451 mu = riskFreeRate_->forwardRate(t0, t0+dt, Continuous).rate()
452 - dividendYield_->forwardRate(t0, t0+dt, Continuous).rate()
453 - 0.5 * vol*vol;
454
455 retVal[1] = varianceDistribution(x0[1], dw[1], dt);
456 dy = (mu - rho_/sigma_*kappa_
457 *(theta_-vol*vol)) * dt + vol*sqrhov*dw[0]*sdt;
458
459 retVal[0] = x0[0]*std::exp(dy + rho_/sigma_*(retVal[1]-x0[1]));
460 break;
463 {
464 // for details of the quadratic exponential discretization scheme
465 // see Leif Andersen,
466 // Efficient Simulation of the Heston Stochastic Volatility Model
467 const Real ex = std::exp(-kappa_*dt);
468
469 const Real m = theta_+(x0[1]-theta_)*ex;
470 const Real s2 = x0[1]*sigma_*sigma_*ex/kappa_*(1-ex)
471 + theta_*sigma_*sigma_/(2*kappa_)*(1-ex)*(1-ex);
472 const Real psi = s2/(m*m);
473
474 const Real g1 = 0.5;
475 const Real g2 = 0.5;
476 Real k0 = -rho_*kappa_*theta_*dt/sigma_;
477 const Real k1 = g1*dt*(kappa_*rho_/sigma_-0.5)-rho_/sigma_;
478 const Real k2 = g2*dt*(kappa_*rho_/sigma_-0.5)+rho_/sigma_;
479 const Real k3 = g1*dt*(1-rho_*rho_);
480 const Real k4 = g2*dt*(1-rho_*rho_);
481 const Real A = k2+0.5*k4;
482
483 if (psi < 1.5) {
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);
487
489 // martingale correction
490 QL_REQUIRE(A < 1/(2*a), "illegal value");
491 k0 = -A*b2*a/(1-2*A*a)+0.5*std::log(1-2*A*a)
492 -(k1+0.5*k3)*x0[1];
493 }
494 retVal[1] = a*(b+dw[1])*(b+dw[1]);
495 }
496 else {
497 const Real p = (psi-1)/(psi+1);
498 const Real beta = (1-p)/m;
499
500 const Real u = CumulativeNormalDistribution()(dw[1]);
501
503 // martingale correction
504 QL_REQUIRE(A < beta, "illegal value");
505 k0 = -std::log(p+beta*(1-p)/(beta-A))-(k1+0.5*k3)*x0[1];
506 }
507 retVal[1] = ((u <= p) ? Real(0.0) : std::log((1-p)/(1-u))/beta);
508 }
509
510 mu = riskFreeRate_->forwardRate(t0, t0+dt, Continuous).rate()
511 - dividendYield_->forwardRate(t0, t0+dt, Continuous).rate();
512
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]);
515 }
516 break;
520 {
521 const Real nu_0 = x0[1];
522 const Real nu_t = varianceDistribution(nu_0, dw[1], dt);
523
524 const Real x = std::min(1.0-QL_EPSILON,
525 std::max(0.0, CumulativeNormalDistribution()(dw[2])));
526
527 const Real vds = Brent().solve(
528 [&](Real xi){ return cdf_nu_ds_minus_x(*this, xi, nu_0, nu_t, dt, discretization_, x); },
529 1e-5, theta_*dt, 0.1*theta_*dt);
530
531 const Real vdw
532 = (nu_t - nu_0 - kappa_*theta_*dt + kappa_*vds)/sigma_;
533
534 mu = ( riskFreeRate_->forwardRate(t0, t0+dt, Continuous).rate()
535 -dividendYield_->forwardRate(t0, t0+dt, Continuous).rate())*dt
536 - 0.5*vds + rho_*vdw;
537
538 const Volatility sig = std::sqrt((1-rho_*rho_)*vds);
539 const Real s = x0[0]*std::exp(mu + sig*dw[0]);
540
541 retVal[0] = s;
542 retVal[1] = nu_t;
543 }
544 break;
545 default:
546 QL_FAIL("unknown discretization schema");
547 }
548
549 return retVal;
550 }
551
553 return s0_;
554 }
555
557 return dividendYield_;
558 }
559
561 return riskFreeRate_;
562 }
563
565 return riskFreeRate_->dayCounter().yearFraction(
566 riskFreeRate_->referenceDate(), d);
567 }
568
570 const Real df = 4*theta_*kappa_/(sigma_*sigma_);
571 const Real ncp = 4*kappa_*std::exp(-kappa_*dt)
572 /(sigma_*sigma_*(1-std::exp(-kappa_*dt)))*v;
573
574 const Real p = std::min(1.0-QL_EPSILON,
575 std::max(0.0, CumulativeNormalDistribution()(dw)));
576
577 return sigma_*sigma_*(1-std::exp(-kappa_*dt))/(4*kappa_)
579 }
580}
const Real kappa_
Brent 1-D solver.
Chi-square (central and non-central) distributions.
1-D array used in linear algebra.
Definition: array.hpp:52
Brent 1-D solver
Definition: brent.hpp:37
Cumulative normal distribution function.
Concrete date class.
Definition: date.hpp:125
Euler discretization for stochastic processes.
Shared handle to an observable.
Definition: handle.hpp:41
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
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.
Definition: matrix.hpp:41
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
Integral of a one-dimensional function.
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
discretization of a stochastic process over a given time interval
multi-dimensional stochastic process class.
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Euler discretization for stochastic processes.
Date d
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
#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 Volatility
volatility
Definition: types.hpp:78
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
Heston stochastic process.
Real kappa
Real theta
Real v0
Real rho
Real sigma
functionals and combinators not included in the STL
#define M_PI_2
#define M_2_PI
#define M_PI
modified Bessel functions of first and second kind
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
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)
STL namespace.
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
ext::shared_ptr< BlackVolTermStructure > v
Real beta
Definition: sabr.cpp:200
Real nu
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
Integral of a one-dimensional function using segment algorithm.
simple quote class