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
20#include <ql/math/distributions/chisquaredistribution.hpp>
21#include <ql/math/distributions/normaldistribution.hpp>
22#include <ql/math/functional.hpp>
23#include <ql/math/integrals/gaussianquadratures.hpp>
24#include <ql/math/integrals/gausslobattointegral.hpp>
25#include <ql/math/integrals/segmentintegral.hpp>
26#include <ql/math/modifiedbessel.hpp>
27#include <ql/math/solvers1d/brent.hpp>
28#include <ql/processes/eulerdiscretization.hpp>
29#include <ql/processes/hestonprocess.hpp>
30#include <ql/quotes/simplequote.hpp>
31#include <boost/math/distributions/non_central_chi_squared.hpp>
32#include <complex>
33#include <utility>
34
35namespace QuantLib {
36
37 HestonProcess::HestonProcess(Handle<YieldTermStructure> riskFreeRate,
38 Handle<YieldTermStructure> dividendYield,
39 Handle<Quote> s0,
40 Real v0,
41 Real kappa,
42 Real theta,
43 Real sigma,
44 Real rho,
45 Discretization d)
46 : StochasticProcess(ext::shared_ptr<discretization>(new EulerDiscretization)),
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
51 registerWith(h: riskFreeRate_);
52 registerWith(h: dividendYield_);
53 registerWith(h: s0_);
54 }
55
56 Size HestonProcess::size() const {
57 return 2;
58 }
59
60 Size HestonProcess::factors() const {
61 return ( discretization_ == BroadieKayaExactSchemeLobatto
62 || discretization_ == BroadieKayaExactSchemeTrapezoidal
63 || discretization_ == BroadieKayaExactSchemeLaguerre) ? 3 : 2;
64 }
65
66 Array HestonProcess::initialValues() const {
67 return { s0_->value(), v0_ };
68 }
69
70 Array HestonProcess::drift(Time t, const Array& x) const {
71 const Real vol = (x[1] > 0.0) ? std::sqrt(x: x[1])
72 : (discretization_ == Reflection) ? Real(- std::sqrt(x: -x[1]))
73 : 0.0;
74
75 return {
76 riskFreeRate_->forwardRate(t1: t, t2: t, comp: Continuous).rate()
77 - dividendYield_->forwardRate(t1: t, t2: t, comp: Continuous).rate()
78 - 0.5 * vol * vol,
79 kappa_* (theta_-((discretization_==PartialTruncation) ? x[1] : vol*vol))
80 };
81 }
82
83 Matrix HestonProcess::diffusion(Time, const Array& x) const {
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: x[1])
93 : (discretization_ == Reflection) ? Real(-std::sqrt(x: -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(x: 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
104 Array HestonProcess::apply(const Array& x0,
105 const Array& dx) const {
106 return {
107 x0[0] * std::exp(x: 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 z: 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(z: -0.5*ga*dt)/(1.0-std::exp(z: -ga*dt));
138 const std::complex<Real> log_z
139 = -0.5*ga*dt + std::log(z: ga/(1.0-std::exp(z: -ga*dt)));
140
141 const std::complex<Real> alpha
142 = 4.0*ga*std::exp(z: -0.5*ga*dt)/(sigma2*(1.0-std::exp(z: -ga*dt)));
143 const std::complex<Real> beta = 4.0*kappa*std::exp(x: -0.5*kappa*dt)
144 /(sigma2*(1.0-std::exp(x: -kappa*dt)));
145
146 return ga*std::exp(z: -0.5*(ga-kappa)*dt)*(1-std::exp(x: -kappa*dt))
147 / (kappa*(1.0-std::exp(z: -ga*dt)))
148 *std::exp(z: (nu_0+nu_t)/sigma2 * (
149 kappa*(1.0+std::exp(x: -kappa*dt))/(1.0-std::exp(x: -kappa*dt))
150 - ga*(1.0+std::exp(z: -ga*dt))/(1.0-std::exp(z: -ga*dt))))
151 *std::exp(z: nu*log_z)/std::pow(x: z, y: nu)
152 *((nu_t > 1e-8)
153 ? modifiedBesselFunction_i(
154 nu, z: std::sqrt(x: nu_0*nu_t)*alpha)
155 / modifiedBesselFunction_i(
156 nu, z: std::sqrt(x: nu_0*nu_t)*beta)
157 : std::pow(x: alpha/beta, y: 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(x: u*x)/u
164 * Phi(process, a: 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(x: u*x)*Phi(process, a: 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(x: process.s0()->value());
180
181 return gaussLaguerreIntegration(
182 [&](Real u){ return ph(process, x: y, u, nu_0, nu_t, dt: t); })
183 / std::sqrt(x: 2*M_PI*(1-rho*rho)*y)
184 * std::exp(x: -0.5*squared(x: 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*x, nominator: n, denominator: d, m: 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(x: y, nominator: fn, denominator: fd, m: 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(x: y, nominator: gn, denominator: gd, m: 10);
244
245 return M_PI_2 - f*std::cos(x: x)-g*std::sin(x: 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, a: std::complex<Real>(0, -2*d),
255 nu_0, nu_t, dt).real();
256 const Real p1 = Phi(process, a: std::complex<Real>(0, -d),
257 nu_0, nu_t, dt).real();
258 const Real p0 = Phi(process, a: std::complex<Real>(0, 0),
259 nu_0, nu_t, dt).real();
260 const Real pm1= Phi(process, a: std::complex<Real>(0, d),
261 nu_0, nu_t, dt).real();
262 const Real pm2= Phi(process, a: 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(x: 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(a: 100.0,
293 b: std::max(a: 0.1, b: cornishFisherEps(process, nu_0, nu_t, dt, eps)));
294
295 switch (discretization) {
296 case HestonProcess::BroadieKayaExactSchemeLaguerre:
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(z: Phi(process,a: upper,nu_0,nu_t,dt)/upper)
304 > eps) upper*=2.0;
305
306 return (x < upper)
307 ? std::max(a: 0.0, b: std::min(a: 1.0,
308 b: gaussLaguerreIntegration(
309 [&](Real u){ return ch(process, x, u, nu_0, nu_t, dt); })))
310 : Real(1.0);
311 }
312 case HestonProcess::BroadieKayaExactSchemeLobatto:
313 {
314 // get the upper bound for the integration
315 Real upper = u_eps/2.0;
316 while (std::abs(z: Phi(process, a: upper,nu_0,nu_t,dt)/upper)
317 > eps) upper*=2.0;
318
319 return (x < upper)
320 ? std::max(a: 0.0, b: std::min(a: 1.0,
321 b: GaussLobattoIntegral(Null<Size>(), eps)(
322 [&](Real xi){ return ch(process, x, u: xi, nu_0, nu_t, dt); },
323 QL_EPSILON, upper)))
324 : Real(1.0);
325 }
326 case HestonProcess::BroadieKayaExactSchemeTrapezoidal:
327 {
328 const Real h = 0.05;
329
330 Real si = Si(x: 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: x*(u+0.5*h));
338
339 f = Phi(process, a: 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(z: f)/j > eps);
344
345 return s;
346 }
347 default:
348 QL_FAIL("unknown integration method");
349 }
350 }
351 }
352
353 Real cdf_nu_ds_minus_x(const HestonProcess &process, Real x, Real nu_0,
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
360 Real HestonProcess::pdf(Real x, Real v, Time t, Real eps) const {
361 const Real k = sigma_*sigma_*(1-std::exp(x: -kappa_*t))/(4*kappa_);
362 const Real a = std::log( x: dividendYield_->discount(t)
363 / riskFreeRate_->discount(t))
364 + rho_/sigma_*(v - v0_ - kappa_*theta_*t);
365
366 const Real x0 = std::log(x: s0()->value());
367 Real upper = std::max(a: 0.1, b: -(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(x: 2*M_PI*(1-rho_*rho_))
374 * ( -0.5/(upper*std::sqrt(x: upper))*std::exp(x: f2)
375 + 1/std::sqrt(x: upper)*std::exp(x: 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(x: f2)/ std::sqrt(x: 2*M_PI*(1-rho_*rho_)*upper);
380 upper*=1.5;
381 }
382
383 upper = 2.0*cornishFisherEps(process: *this, nu_0: v0_, nu_t: v, dt: t,eps: 1e-3);
384
385 return SegmentIntegral(100)(
386 [&](Real xi){ return int_ph(process: *this, a, x, y: xi, nu_0: v0_, nu_t: v, t); },
387 QL_EPSILON, upper)
388 * boost::math::pdf(
389 dist: boost::math::non_central_chi_squared_distribution<Real>(
390 4*theta_*kappa_/(sigma_*sigma_),
391 4*kappa_*std::exp(x: -kappa_*t)
392 /((sigma_*sigma_)*(1-std::exp(x: -kappa_*t)))*v0_),
393 x: v/k) / k;
394 }
395
396 Array HestonProcess::evolve(Time t0, const Array& x0,
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(x: dt);
402 const Real sqrhov = std::sqrt(x: 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
410 case PartialTruncation:
411 vol = (x0[1] > 0.0) ? std::sqrt(x: x0[1]) : Real(0.0);
412 vol2 = sigma_ * vol;
413 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
414 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
415 - 0.5 * vol * vol;
416 nu = kappa_*(theta_ - x0[1]);
417
418 retVal[0] = x0[0] * std::exp(x: 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(x: x0[1]) : Real(0.0);
423 vol2 = sigma_ * vol;
424 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
425 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
426 - 0.5 * vol * vol;
427 nu = kappa_*(theta_ - vol*vol);
428
429 retVal[0] = x0[0] * std::exp(x: 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(x: std::fabs(x: x0[1]));
434 vol2 = sigma_ * vol;
435 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
436 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
437 - 0.5 * vol*vol;
438 nu = kappa_*(theta_ - vol*vol);
439
440 retVal[0] = x0[0]*std::exp(x: mu*dt+vol*dw[0]*sdt);
441 retVal[1] = vol*vol
442 +nu*dt + vol2*sdt*(rho_*dw[0] + sqrhov*dw[1]);
443 break;
444 case NonCentralChiSquareVariance:
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(x: x0[1]) : Real(0.0);
451 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
452 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
453 - 0.5 * vol*vol;
454
455 retVal[1] = varianceDistribution(v: x0[1], dw: 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(x: dy + rho_/sigma_*(retVal[1]-x0[1]));
460 break;
461 case QuadraticExponential:
462 case QuadraticExponentialMartingale:
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(x: -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(x: 2/psi*(2/psi-1));
485 const Real b = std::sqrt(x: b2);
486 const Real a = m/(1+b2);
487
488 if (discretization_ == QuadraticExponentialMartingale) {
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(x: 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
502 if (discretization_ == QuadraticExponentialMartingale) {
503 // martingale correction
504 QL_REQUIRE(A < beta, "illegal value");
505 k0 = -std::log(x: p+beta*(1-p)/(beta-A))-(k1+0.5*k3)*x0[1];
506 }
507 retVal[1] = ((u <= p) ? Real(0.0) : std::log(x: (1-p)/(1-u))/beta);
508 }
509
510 mu = riskFreeRate_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate()
511 - dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate();
512
513 retVal[0] = x0[0]*std::exp(x: mu*dt + k0 + k1*x0[1] + k2*retVal[1]
514 +std::sqrt(x: k3*x0[1]+k4*retVal[1])*dw[0]);
515 }
516 break;
517 case BroadieKayaExactSchemeLobatto:
518 case BroadieKayaExactSchemeLaguerre:
519 case BroadieKayaExactSchemeTrapezoidal:
520 {
521 const Real nu_0 = x0[1];
522 const Real nu_t = varianceDistribution(v: nu_0, dw: dw[1], dt);
523
524 const Real x = std::min(a: 1.0-QL_EPSILON,
525 b: std::max(a: 0.0, b: CumulativeNormalDistribution()(dw[2])));
526
527 const Real vds = Brent().solve(
528 f: [&](Real xi){ return cdf_nu_ds_minus_x(process: *this, x: xi, nu_0, nu_t, dt, discretization: discretization_, x0: x); },
529 accuracy: 1e-5, guess: theta_*dt, step: 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(t1: t0, t2: t0+dt, comp: Continuous).rate()
535 -dividendYield_->forwardRate(t1: t0, t2: t0+dt, comp: Continuous).rate())*dt
536 - 0.5*vds + rho_*vdw;
537
538 const Volatility sig = std::sqrt(x: (1-rho_*rho_)*vds);
539 const Real s = x0[0]*std::exp(x: 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
552 const Handle<Quote>& HestonProcess::s0() const {
553 return s0_;
554 }
555
556 const Handle<YieldTermStructure>& HestonProcess::dividendYield() const {
557 return dividendYield_;
558 }
559
560 const Handle<YieldTermStructure>& HestonProcess::riskFreeRate() const {
561 return riskFreeRate_;
562 }
563
564 Time HestonProcess::time(const Date& d) const {
565 return riskFreeRate_->dayCounter().yearFraction(
566 d1: riskFreeRate_->referenceDate(), d2: d);
567 }
568
569 Real HestonProcess::varianceDistribution(Real v, Real dw, Time dt) const {
570 const Real df = 4*theta_*kappa_/(sigma_*sigma_);
571 const Real ncp = 4*kappa_*std::exp(x: -kappa_*dt)
572 /(sigma_*sigma_*(1-std::exp(x: -kappa_*dt)))*v;
573
574 const Real p = std::min(a: 1.0-QL_EPSILON,
575 b: std::max(a: 0.0, b: CumulativeNormalDistribution()(dw)));
576
577 return sigma_*sigma_*(1-std::exp(x: -kappa_*dt))/(4*kappa_)
578 *InverseNonCentralCumulativeChiSquareDistribution(df, ncp, 100)(p);
579 }
580}
581

source code of quantlib/ql/processes/hestonprocess.cpp

Morty Proxy This is a proxified and sanitized view of the page, visit original site.