1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2004, 2005, 2008 Klaus Spanderen 5 Copyright (C) 2007 StatPro Italia srl 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 /*! \file hestonmodel.hpp 22 \brief analytic pricing engine for a heston option 23 based on fourier transformation 24 */ 25 26 #include <ql/functional.hpp> 27 #include <ql/math/solvers1d/brent.hpp> 28 #include <ql/math/functional.hpp> 29 #include <ql/math/integrals/simpsonintegral.hpp> 30 #include <ql/math/integrals/kronrodintegral.hpp> 31 #include <ql/math/integrals/trapezoidintegral.hpp> 32 #include <ql/math/integrals/discreteintegrals.hpp> 33 #include <ql/math/integrals/gausslobattointegral.hpp> 34 #include <ql/math/integrals/exponentialintegrals.hpp> 35 36 #include <ql/instruments/payoffs.hpp> 37 #include <ql/pricingengines/blackcalculator.hpp> 38 #include <ql/pricingengines/vanilla/analytichestonengine.hpp> 39 40 #if defined(QL_PATCH_MSVC) 41 #pragma warning(disable: 4180) 42 #endif 43 44 namespace QuantLib { 45 46 namespace { 47 48 class integrand1 { 49 private: 50 const Real c_inf_; 51 const ext::function<Real(Real)> f_; 52 public: integrand1(Real c_inf,const ext::function<Real (Real)> & f)53 integrand1(Real c_inf, const ext::function<Real(Real)>& f) 54 : c_inf_(c_inf), f_(f) {} operator ()(Real x) const55 Real operator()(Real x) const { 56 if ((1.0-x)*c_inf_ > QL_EPSILON) 57 return f_(-std::log(0.5-0.5*x)/c_inf_)/((1.0-x)*c_inf_); 58 else 59 return 0.0; 60 } 61 }; 62 63 class integrand2 { 64 private: 65 const Real c_inf_; 66 const ext::function<Real(Real)> f_; 67 public: integrand2(Real c_inf,const ext::function<Real (Real)> & f)68 integrand2(Real c_inf, const ext::function<Real(Real)>& f) 69 : c_inf_(c_inf), f_(f) {} operator ()(Real x) const70 Real operator()(Real x) const { 71 if (x*c_inf_ > QL_EPSILON) { 72 return f_(-std::log(x)/c_inf_)/(x*c_inf_); 73 } else { 74 return 0.0; 75 } 76 } 77 }; 78 79 class integrand3 { 80 private: 81 const integrand2 int_; 82 public: integrand3(Real c_inf,const ext::function<Real (Real)> & f)83 integrand3(Real c_inf, const ext::function<Real(Real)>& f) 84 : int_(c_inf, f) {} 85 operator ()(Real x) const86 Real operator()(Real x) const { return int_(1.0-x); } 87 }; 88 89 class u_Max { 90 public: u_Max(Real c_inf,Real epsilon)91 u_Max(Real c_inf, Real epsilon) 92 : c_inf_(c_inf), logEpsilon_(std::log(epsilon)), 93 evaluations_(0) {} 94 operator ()(Real u) const95 Real operator()(Real u) const { 96 ++evaluations_; 97 return c_inf_*u + std::log(u) + logEpsilon_; 98 } 99 evaluations() const100 Size evaluations() const { return evaluations_; } 101 102 private: 103 const Real c_inf_, logEpsilon_; 104 mutable Size evaluations_; 105 }; 106 107 108 class uHat_Max { 109 public: uHat_Max(Real v0T2,Real epsilon)110 uHat_Max(Real v0T2, Real epsilon) 111 : v0T2_(v0T2), logEpsilon_(std::log(epsilon)), 112 evaluations_(0) {} 113 operator ()(Real u) const114 Real operator()(Real u) const { 115 ++evaluations_; 116 return v0T2_*u*u + std::log(u) + logEpsilon_; 117 } 118 evaluations() const119 Size evaluations() const { return evaluations_; } 120 121 private: 122 const Real v0T2_, logEpsilon_; 123 mutable Size evaluations_; 124 }; 125 } 126 127 // helper class for integration 128 class AnalyticHestonEngine::Fj_Helper { 129 public: 130 Fj_Helper(const VanillaOption::arguments& arguments, 131 const ext::shared_ptr<HestonModel>& model, 132 const AnalyticHestonEngine* engine, 133 ComplexLogFormula cpxLog, 134 Time term, 135 Real ratio, 136 Size j); 137 138 Fj_Helper(Real kappa, 139 Real theta, 140 Real sigma, 141 Real v0, 142 Real s0, 143 Real rho, 144 const AnalyticHestonEngine* engine, 145 ComplexLogFormula cpxLog, 146 Time term, 147 Real strike, 148 Real ratio, 149 Size j); 150 151 Fj_Helper(Real kappa, 152 Real theta, 153 Real sigma, 154 Real v0, 155 Real s0, 156 Real rho, 157 ComplexLogFormula cpxLog, 158 Time term, 159 Real strike, 160 Real ratio, 161 Size j); 162 163 Real operator()(Real phi) const; 164 165 private: 166 const Size j_; 167 // const VanillaOption::arguments& arg_; 168 const Real kappa_, theta_, sigma_, v0_; 169 const ComplexLogFormula cpxLog_; 170 171 // helper variables 172 const Time term_; 173 const Real x_, sx_, dd_; 174 const Real sigma2_, rsigma_; 175 const Real t0_; 176 177 // log branch counter 178 mutable int b_; // log branch counter 179 mutable Real g_km1_; // imag part of last log value 180 181 const AnalyticHestonEngine* const engine_; 182 }; 183 184 Fj_Helper(const VanillaOption::arguments & arguments,const ext::shared_ptr<HestonModel> & model,const AnalyticHestonEngine * const engine,ComplexLogFormula cpxLog,Time term,Real ratio,Size j)185 AnalyticHestonEngine::Fj_Helper::Fj_Helper( 186 const VanillaOption::arguments& arguments, 187 const ext::shared_ptr<HestonModel>& model, 188 const AnalyticHestonEngine* const engine, 189 ComplexLogFormula cpxLog, 190 Time term, Real ratio, Size j) 191 : j_ (j), //arg_(arguments), 192 kappa_(model->kappa()), theta_(model->theta()), 193 sigma_(model->sigma()), v0_(model->v0()), 194 cpxLog_(cpxLog), term_(term), 195 x_(std::log(model->process()->s0()->value())), 196 sx_(std::log(ext::dynamic_pointer_cast<StrikedTypePayoff> 197 (arguments.payoff)->strike())), 198 dd_(x_-std::log(ratio)), 199 sigma2_(sigma_*sigma_), 200 rsigma_(model->rho()*sigma_), 201 t0_(kappa_ - ((j_== 1)? model->rho()*sigma_ : 0)), 202 b_(0), g_km1_(0), 203 engine_(engine) 204 { 205 } 206 Fj_Helper(Real kappa,Real theta,Real sigma,Real v0,Real s0,Real rho,const AnalyticHestonEngine * const engine,ComplexLogFormula cpxLog,Time term,Real strike,Real ratio,Size j)207 AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa, Real theta, 208 Real sigma, Real v0, Real s0, Real rho, 209 const AnalyticHestonEngine* const engine, 210 ComplexLogFormula cpxLog, 211 Time term, 212 Real strike, 213 Real ratio, 214 Size j) 215 : 216 j_(j), 217 kappa_(kappa), 218 theta_(theta), 219 sigma_(sigma), 220 v0_(v0), 221 cpxLog_(cpxLog), 222 term_(term), 223 x_(std::log(s0)), 224 sx_(std::log(strike)), 225 dd_(x_-std::log(ratio)), 226 sigma2_(sigma_*sigma_), 227 rsigma_(rho*sigma_), 228 t0_(kappa - ((j== 1)? rho*sigma : 0)), 229 b_(0), 230 g_km1_(0), 231 engine_(engine) 232 { 233 } 234 Fj_Helper(Real kappa,Real theta,Real sigma,Real v0,Real s0,Real rho,ComplexLogFormula cpxLog,Time term,Real strike,Real ratio,Size j)235 AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa, Real theta, 236 Real sigma, Real v0, Real s0, Real rho, 237 ComplexLogFormula cpxLog, 238 Time term, 239 Real strike, 240 Real ratio, 241 Size j) 242 : 243 j_(j), 244 kappa_(kappa), 245 theta_(theta), 246 sigma_(sigma), 247 v0_(v0), 248 cpxLog_(cpxLog), 249 term_(term), 250 x_(std::log(s0)), 251 sx_(std::log(strike)), 252 dd_(x_-std::log(ratio)), 253 sigma2_(sigma_*sigma_), 254 rsigma_(rho*sigma_), 255 t0_(kappa - ((j== 1)? rho*sigma : 0)), 256 b_(0), 257 g_km1_(0), 258 engine_(0) 259 { 260 } 261 262 operator ()(Real phi) const263 Real AnalyticHestonEngine::Fj_Helper::operator()(Real phi) const 264 { 265 const Real rpsig(rsigma_*phi); 266 267 const std::complex<Real> t1 = t0_+std::complex<Real>(0, -rpsig); 268 const std::complex<Real> d = 269 std::sqrt(t1*t1 - sigma2_*phi 270 *std::complex<Real>(-phi, (j_== 1)? 1 : -1)); 271 const std::complex<Real> ex = std::exp(-d*term_); 272 const std::complex<Real> addOnTerm = 273 engine_ != 0 ? engine_->addOnTerm(phi, term_, j_) : Real(0.0); 274 275 if (cpxLog_ == Gatheral) { 276 if (phi != 0.0) { 277 if (sigma_ > 1e-5) { 278 const std::complex<Real> p = (t1-d)/(t1+d); 279 const std::complex<Real> g 280 = std::log((1.0 - p*ex)/(1.0 - p)); 281 282 return 283 std::exp(v0_*(t1-d)*(1.0-ex)/(sigma2_*(1.0-ex*p)) 284 + (kappa_*theta_)/sigma2_*((t1-d)*term_-2.0*g) 285 + std::complex<Real>(0.0, phi*(dd_-sx_)) 286 + addOnTerm 287 ).imag()/phi; 288 } 289 else { 290 const std::complex<Real> td = phi/(2.0*t1) 291 *std::complex<Real>(-phi, (j_== 1)? 1 : -1); 292 const std::complex<Real> p = td*sigma2_/(t1+d); 293 const std::complex<Real> g = p*(1.0-ex); 294 295 return 296 std::exp(v0_*td*(1.0-ex)/(1.0-p*ex) 297 + (kappa_*theta_)*(td*term_-2.0*g/sigma2_) 298 + std::complex<Real>(0.0, phi*(dd_-sx_)) 299 + addOnTerm 300 ).imag()/phi; 301 } 302 } 303 else { 304 // use l'Hospital's rule to get lim_{phi->0} 305 if (j_ == 1) { 306 const Real kmr = rsigma_-kappa_; 307 if (std::fabs(kmr) > 1e-7) { 308 return dd_-sx_ 309 + (std::exp(kmr*term_)*kappa_*theta_ 310 -kappa_*theta_*(kmr*term_+1.0) ) / (2*kmr*kmr) 311 - v0_*(1.0-std::exp(kmr*term_)) / (2.0*kmr); 312 } 313 else 314 // \kappa = \rho * \sigma 315 return dd_-sx_ + 0.25*kappa_*theta_*term_*term_ 316 + 0.5*v0_*term_; 317 } 318 else { 319 return dd_-sx_ 320 - (std::exp(-kappa_*term_)*kappa_*theta_ 321 +kappa_*theta_*(kappa_*term_-1.0))/(2*kappa_*kappa_) 322 - v0_*(1.0-std::exp(-kappa_*term_))/(2*kappa_); 323 } 324 } 325 } 326 else if (cpxLog_ == BranchCorrection) { 327 const std::complex<Real> p = (t1+d)/(t1-d); 328 329 // next term: g = std::log((1.0 - p*std::exp(d*term_))/(1.0 - p)) 330 std::complex<Real> g; 331 332 // the exp of the following expression is needed. 333 const std::complex<Real> e = std::log(p)+d*term_; 334 335 // does it fit to the machine precision? 336 if (std::exp(-e.real()) > QL_EPSILON) { 337 g = std::log((1.0 - p/ex)/(1.0 - p)); 338 } else { 339 // use a "big phi" approximation 340 g = d*term_ + std::log(p/(p - 1.0)); 341 342 if (g.imag() > M_PI || g.imag() <= -M_PI) { 343 // get back to principal branch of the complex logarithm 344 Real im = std::fmod(g.imag(), 2*M_PI); 345 if (im > M_PI) 346 im -= 2*M_PI; 347 else if (im <= -M_PI) 348 im += 2*M_PI; 349 350 g = std::complex<Real>(g.real(), im); 351 } 352 } 353 354 // be careful here as we have to use a log branch correction 355 // to deal with the discontinuities of the complex logarithm. 356 // the principal branch is not always the correct one. 357 // (s. A. Sepp, chapter 4) 358 // remark: there is still the change that we miss a branch 359 // if the order of the integration is not high enough. 360 const Real tmp = g.imag() - g_km1_; 361 if (tmp <= -M_PI) 362 ++b_; 363 else if (tmp > M_PI) 364 --b_; 365 366 g_km1_ = g.imag(); 367 g += std::complex<Real>(0, 2*b_*M_PI); 368 369 return std::exp(v0_*(t1+d)*(ex-1.0)/(sigma2_*(ex-p)) 370 + (kappa_*theta_)/sigma2_*((t1+d)*term_-2.0*g) 371 + std::complex<Real>(0,phi*(dd_-sx_)) 372 + addOnTerm 373 ).imag()/phi; 374 } 375 else { 376 QL_FAIL("unknown complex logarithm formula"); 377 } 378 } 379 380 AP_Helper(Time term,Real fwd,Real strike,ComplexLogFormula cpxLog,const AnalyticHestonEngine * const enginePtr)381 AnalyticHestonEngine::AP_Helper::AP_Helper( 382 Time term, Real fwd, Real strike, ComplexLogFormula cpxLog, 383 const AnalyticHestonEngine* const enginePtr) 384 : term_(term), 385 fwd_(fwd), 386 strike_(strike), 387 freq_(std::log(fwd/strike)), 388 cpxLog_(cpxLog), 389 enginePtr_(enginePtr) { 390 QL_REQUIRE(enginePtr != 0, "pricing engine required"); 391 392 const Real v0 = enginePtr->model_->v0(); 393 const Real kappa = enginePtr->model_->kappa(); 394 const Real theta = enginePtr->model_->theta(); 395 const Real sigma = enginePtr->model_->sigma(); 396 const Real rho = enginePtr->model_->rho(); 397 398 switch(cpxLog_) { 399 case AndersenPiterbarg: 400 vAvg_ = (1-std::exp(-kappa*term))*(v0 - theta) 401 /(kappa*term) + theta; 402 break; 403 case AndersenPiterbargOptCV: 404 vAvg_ = -8.0*std::log(enginePtr->chF( 405 std::complex<Real>(0, -0.5), term).real())/term; 406 break; 407 case AsymptoticChF: 408 phi_ = -(v0+term*kappa*theta)/sigma 409 * std::complex<Real>(std::sqrt(1-rho*rho), rho); 410 411 psi_ = std::complex<Real>( 412 (kappa- 0.5*rho*sigma)*(v0 + term*kappa*theta) 413 + kappa*theta*std::log(4*(1-rho*rho)), 414 - ((0.5*rho*rho*sigma - kappa*rho)/std::sqrt(1-rho*rho) 415 *(v0 + kappa*theta*term) 416 - 2*kappa*theta*std::atan(rho/std::sqrt(1-rho*rho)))) 417 /(sigma*sigma); 418 break; 419 default: 420 QL_FAIL("unknown control variate"); 421 } 422 } 423 operator ()(Real u) const424 Real AnalyticHestonEngine::AP_Helper::operator()(Real u) const { 425 QL_REQUIRE( enginePtr_->addOnTerm(u, term_, 1) 426 == std::complex<Real>(0.0) 427 && enginePtr_->addOnTerm(u, term_, 2) 428 == std::complex<Real>(0.0), 429 "only Heston model is supported"); 430 431 const std::complex<Real> z(u, -0.5); 432 433 std::complex<Real> phiBS; 434 435 switch (cpxLog_) { 436 case AndersenPiterbarg: 437 case AndersenPiterbargOptCV: 438 phiBS = std::exp( 439 -0.5*vAvg_*term_*(z*z + std::complex<Real>(-z.imag(), z.real()))); 440 break; 441 case AsymptoticChF: 442 phiBS = std::exp(u*phi_ + psi_); 443 break; 444 default: 445 QL_FAIL("unknown control variate"); 446 } 447 448 return (std::exp(std::complex<Real>(0.0, u*freq_)) 449 * (phiBS - enginePtr_->chF(z, term_)) / (u*u + 0.25)).real(); 450 } 451 controlVariateValue() const452 Real AnalyticHestonEngine::AP_Helper::controlVariateValue() const { 453 if (cpxLog_ == AndersenPiterbarg || cpxLog_ == AndersenPiterbargOptCV) { 454 return BlackCalculator( 455 Option::Call, strike_, fwd_, std::sqrt(vAvg_*term_)) 456 .value(); 457 } 458 else if (cpxLog_ == AsymptoticChF) { 459 const std::complex<Real> phiFreq(phi_.real(), phi_.imag() + freq_); 460 461 using namespace ExponentialIntegral; 462 return fwd_ - std::sqrt(strike_*fwd_)/M_PI* 463 (std::exp(psi_)*( 464 -2.0*Ci(-0.5*phiFreq)*std::sin(0.5*phiFreq) 465 +std::cos(0.5*phiFreq)*(M_PI+2.0*Si(0.5*phiFreq)))).real(); 466 } 467 else 468 QL_FAIL("unknown control variate"); 469 } 470 chF(const std::complex<Real> & z,Time t) const471 std::complex<Real> AnalyticHestonEngine::chF( 472 const std::complex<Real>& z, Time t) const { 473 474 const Real kappa = model_->kappa(); 475 const Real sigma = model_->sigma(); 476 const Real theta = model_->theta(); 477 const Real rho = model_->rho(); 478 const Real v0 = model_->v0(); 479 480 const Real sigma2 = sigma*sigma; 481 482 if (sigma > 1e-4) { 483 const std::complex<Real> g 484 = kappa + rho*sigma*std::complex<Real>(z.imag(), -z.real()); 485 486 const std::complex<Real> D = std::sqrt( 487 g*g + (z*z + std::complex<Real>(-z.imag(), z.real()))*sigma2); 488 489 const std::complex<Real> G = (g-D)/(g+D); 490 491 return std::exp(v0/sigma2*(1.0-std::exp(-D*t))/(1.0-G*std::exp(-D*t)) 492 *(g-D) + kappa*theta/sigma2*((g-D)*t 493 -2.0*std::log((1.0-G*std::exp(-D*t))/(1.0-G)))); 494 } 495 else { 496 const Real kt = kappa*t; 497 const Real ekt = std::exp(kt); 498 const Real e2kt = std::exp(2*kt); 499 const Real rho2 = rho*rho; 500 const std::complex<Real> zpi = z + std::complex<Real>(0.0, 1.0); 501 502 return std::exp(-(((theta - v0 + ekt*((-1 + kt)*theta + v0)) 503 *z*zpi)/ekt)/(2.*kappa)) 504 + (std::exp(-(kt) - ((theta - v0 + ekt 505 *((-1 + kt)*theta + v0))*z*zpi) 506 /(2.*ekt*kappa))*rho*(2*theta + kt*theta - 507 v0 - kt*v0 + ekt*((-2 + kt)*theta + v0)) 508 *(1.0 - std::complex<Real>(-z.imag(),z.real()))*z*z) 509 /(2.*kappa*kappa)*sigma 510 + (std::exp(-2*kt - ((theta - v0 + ekt 511 *((-1 + kt)*theta + v0))*z*zpi)/(2.*ekt*kappa))*z*z*zpi 512 *(-2*rho2*square<Real>()(2*theta + kt*theta - v0 - 513 kt*v0 + ekt*((-2 + kt)*theta + v0)) 514 *z*z*zpi + 2*kappa*v0*(-zpi 515 + e2kt*(zpi + 4*rho2*z) - 2*ekt*(2*rho2*z 516 + kt*(zpi + rho2*(2 + kt)*z))) + kappa*theta*(zpi + e2kt 517 *(-5.0*zpi - 24*rho2*z+ 2*kt*(zpi + 4*rho2*z)) + 518 4*ekt*(zpi + 6*rho2*z + kt*(zpi + rho2*(4 + kt)*z))))) 519 /(16.*square<Real>()(square<Real>()(kappa)))*sigma2; 520 } 521 } 522 lnChF(const std::complex<Real> & z,Time T) const523 std::complex<Real> AnalyticHestonEngine::lnChF( 524 const std::complex<Real>& z, Time T) const { 525 return std::log(chF(z, T)); 526 } 527 AnalyticHestonEngine(const ext::shared_ptr<HestonModel> & model,Size integrationOrder)528 AnalyticHestonEngine::AnalyticHestonEngine( 529 const ext::shared_ptr<HestonModel>& model, 530 Size integrationOrder) 531 : GenericModelEngine<HestonModel, 532 VanillaOption::arguments, 533 VanillaOption::results>(model), 534 evaluations_(0), 535 cpxLog_ (Gatheral), 536 integration_(new Integration( 537 Integration::gaussLaguerre(integrationOrder))), 538 andersenPiterbargEpsilon_(Null<Real>()) { 539 } 540 AnalyticHestonEngine(const ext::shared_ptr<HestonModel> & model,Real relTolerance,Size maxEvaluations)541 AnalyticHestonEngine::AnalyticHestonEngine( 542 const ext::shared_ptr<HestonModel>& model, 543 Real relTolerance, Size maxEvaluations) 544 : GenericModelEngine<HestonModel, 545 VanillaOption::arguments, 546 VanillaOption::results>(model), 547 evaluations_(0), 548 cpxLog_(Gatheral), 549 integration_(new Integration(Integration::gaussLobatto( 550 relTolerance, Null<Real>(), maxEvaluations))), 551 andersenPiterbargEpsilon_(Null<Real>()) { 552 } 553 AnalyticHestonEngine(const ext::shared_ptr<HestonModel> & model,ComplexLogFormula cpxLog,const Integration & integration,const Real andersenPiterbargEpsilon)554 AnalyticHestonEngine::AnalyticHestonEngine( 555 const ext::shared_ptr<HestonModel>& model, 556 ComplexLogFormula cpxLog, 557 const Integration& integration, 558 const Real andersenPiterbargEpsilon) 559 : GenericModelEngine<HestonModel, 560 VanillaOption::arguments, 561 VanillaOption::results>(model), 562 evaluations_(0), 563 cpxLog_(cpxLog), 564 integration_(new Integration(integration)), 565 andersenPiterbargEpsilon_(andersenPiterbargEpsilon) { 566 QL_REQUIRE( cpxLog_ != BranchCorrection 567 || !integration.isAdaptiveIntegration(), 568 "Branch correction does not work in conjunction " 569 "with adaptive integration methods"); 570 } 571 572 AnalyticHestonEngine::ComplexLogFormula optimalControlVariate(Time t,Real v0,Real kappa,Real theta,Real sigma,Real rho)573 AnalyticHestonEngine::optimalControlVariate( 574 Time t, Real v0, Real kappa, Real theta, Real sigma, Real rho) { 575 576 if (t > 0.1 && (v0+t*kappa*theta)/sigma*std::sqrt(1-rho*rho) < 0.055) { 577 return AsymptoticChF; 578 } 579 else { 580 return AndersenPiterbargOptCV; 581 } 582 } 583 584 numberOfEvaluations() const585 Size AnalyticHestonEngine::numberOfEvaluations() const { 586 return evaluations_; 587 } 588 doCalculation(Real riskFreeDiscount,Real dividendDiscount,Real spotPrice,Real strikePrice,Real term,Real kappa,Real theta,Real sigma,Real v0,Real rho,const TypePayoff & type,const Integration & integration,const ComplexLogFormula cpxLog,const AnalyticHestonEngine * const enginePtr,Real & value,Size & evaluations)589 void AnalyticHestonEngine::doCalculation(Real riskFreeDiscount, 590 Real dividendDiscount, 591 Real spotPrice, 592 Real strikePrice, 593 Real term, 594 Real kappa, Real theta, Real sigma, Real v0, Real rho, 595 const TypePayoff& type, 596 const Integration& integration, 597 const ComplexLogFormula cpxLog, 598 const AnalyticHestonEngine* const enginePtr, 599 Real& value, 600 Size& evaluations) { 601 602 const Real ratio = riskFreeDiscount/dividendDiscount; 603 604 evaluations = 0; 605 606 switch(cpxLog) { 607 case Gatheral: 608 case BranchCorrection: { 609 const Real c_inf = std::min(0.2, std::max(0.0001, 610 std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*term); 611 612 const Real p1 = integration.calculate(c_inf, 613 Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr, 614 cpxLog, term, strikePrice, ratio, 1))/M_PI; 615 evaluations += integration.numberOfEvaluations(); 616 617 const Real p2 = integration.calculate(c_inf, 618 Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr, 619 cpxLog, term, strikePrice, ratio, 2))/M_PI; 620 evaluations += integration.numberOfEvaluations(); 621 622 switch (type.optionType()) 623 { 624 case Option::Call: 625 value = spotPrice*dividendDiscount*(p1+0.5) 626 - strikePrice*riskFreeDiscount*(p2+0.5); 627 break; 628 case Option::Put: 629 value = spotPrice*dividendDiscount*(p1-0.5) 630 - strikePrice*riskFreeDiscount*(p2-0.5); 631 break; 632 default: 633 QL_FAIL("unknown option type"); 634 } 635 } 636 break; 637 case AndersenPiterbarg: 638 case AndersenPiterbargOptCV: 639 case AsymptoticChF: 640 case OptimalCV: { 641 const Real c_inf = 642 std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*term)/sigma; 643 644 const Real fwdPrice = spotPrice / ratio; 645 646 const Real epsilon = enginePtr->andersenPiterbargEpsilon_ 647 *M_PI/(std::sqrt(strikePrice*fwdPrice)*riskFreeDiscount); 648 649 const ext::function<Real()> uM = ext::bind( 650 Integration::andersenPiterbargIntegrationLimit, 651 c_inf, epsilon, v0, term); 652 653 AP_Helper cvHelper(term, fwdPrice, strikePrice, 654 (cpxLog == OptimalCV) 655 ? optimalControlVariate(term, v0, kappa, theta, sigma, rho) 656 : cpxLog, 657 enginePtr 658 ); 659 660 const Real cvValue = cvHelper.controlVariateValue(); 661 662 const Real h_cv = integration.calculate(c_inf, cvHelper, uM) 663 * std::sqrt(strikePrice * fwdPrice)/M_PI; 664 evaluations += integration.numberOfEvaluations(); 665 666 switch (type.optionType()) 667 { 668 case Option::Call: 669 value = (cvValue + h_cv)*riskFreeDiscount; 670 break; 671 case Option::Put: 672 value = (cvValue + h_cv - (fwdPrice - strikePrice))*riskFreeDiscount; 673 break; 674 default: 675 QL_FAIL("unknown option type"); 676 } 677 } 678 break; 679 680 default: 681 QL_FAIL("unknown complex log formula"); 682 } 683 } 684 calculate() const685 void AnalyticHestonEngine::calculate() const 686 { 687 // this is a european option pricer 688 QL_REQUIRE(arguments_.exercise->type() == Exercise::European, 689 "not an European option"); 690 691 // plain vanilla 692 ext::shared_ptr<PlainVanillaPayoff> payoff = 693 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 694 QL_REQUIRE(payoff, "non plain vanilla payoff given"); 695 696 const ext::shared_ptr<HestonProcess>& process = model_->process(); 697 698 const Real riskFreeDiscount = process->riskFreeRate()->discount( 699 arguments_.exercise->lastDate()); 700 const Real dividendDiscount = process->dividendYield()->discount( 701 arguments_.exercise->lastDate()); 702 703 const Real spotPrice = process->s0()->value(); 704 QL_REQUIRE(spotPrice > 0.0, "negative or null underlying given"); 705 706 const Real strikePrice = payoff->strike(); 707 const Real term = process->time(arguments_.exercise->lastDate()); 708 709 doCalculation(riskFreeDiscount, 710 dividendDiscount, 711 spotPrice, 712 strikePrice, 713 term, 714 model_->kappa(), 715 model_->theta(), 716 model_->sigma(), 717 model_->v0(), 718 model_->rho(), 719 *payoff, 720 *integration_, 721 cpxLog_, 722 this, 723 results_.value, 724 evaluations_); 725 } 726 727 Integration(Algorithm intAlgo,const ext::shared_ptr<Integrator> & integrator)728 AnalyticHestonEngine::Integration::Integration( 729 Algorithm intAlgo, 730 const ext::shared_ptr<Integrator>& integrator) 731 : intAlgo_(intAlgo), 732 integrator_(integrator) { } 733 Integration(Algorithm intAlgo,const ext::shared_ptr<GaussianQuadrature> & gaussianQuadrature)734 AnalyticHestonEngine::Integration::Integration( 735 Algorithm intAlgo, 736 const ext::shared_ptr<GaussianQuadrature>& gaussianQuadrature) 737 : intAlgo_(intAlgo), 738 gaussianQuadrature_(gaussianQuadrature) { } 739 740 AnalyticHestonEngine::Integration gaussLobatto(Real relTolerance,Real absTolerance,Size maxEvaluations)741 AnalyticHestonEngine::Integration::gaussLobatto(Real relTolerance, 742 Real absTolerance, 743 Size maxEvaluations) { 744 return Integration(GaussLobatto, 745 ext::shared_ptr<Integrator>( 746 new GaussLobattoIntegral(maxEvaluations, 747 absTolerance, 748 relTolerance, 749 false))); 750 } 751 752 AnalyticHestonEngine::Integration gaussKronrod(Real absTolerance,Size maxEvaluations)753 AnalyticHestonEngine::Integration::gaussKronrod(Real absTolerance, 754 Size maxEvaluations) { 755 return Integration(GaussKronrod, 756 ext::shared_ptr<Integrator>( 757 new GaussKronrodAdaptive(absTolerance, 758 maxEvaluations))); 759 } 760 761 AnalyticHestonEngine::Integration simpson(Real absTolerance,Size maxEvaluations)762 AnalyticHestonEngine::Integration::simpson(Real absTolerance, 763 Size maxEvaluations) { 764 return Integration(Simpson, 765 ext::shared_ptr<Integrator>( 766 new SimpsonIntegral(absTolerance, 767 maxEvaluations))); 768 } 769 770 AnalyticHestonEngine::Integration trapezoid(Real absTolerance,Size maxEvaluations)771 AnalyticHestonEngine::Integration::trapezoid(Real absTolerance, 772 Size maxEvaluations) { 773 return Integration(Trapezoid, 774 ext::shared_ptr<Integrator>( 775 new TrapezoidIntegral<Default>(absTolerance, 776 maxEvaluations))); 777 } 778 779 AnalyticHestonEngine::Integration gaussLaguerre(Size intOrder)780 AnalyticHestonEngine::Integration::gaussLaguerre(Size intOrder) { 781 QL_REQUIRE(intOrder <= 192, "maximum integraton order (192) exceeded"); 782 return Integration(GaussLaguerre, 783 ext::shared_ptr<GaussianQuadrature>( 784 new GaussLaguerreIntegration(intOrder))); 785 } 786 787 AnalyticHestonEngine::Integration gaussLegendre(Size intOrder)788 AnalyticHestonEngine::Integration::gaussLegendre(Size intOrder) { 789 return Integration(GaussLegendre, 790 ext::shared_ptr<GaussianQuadrature>( 791 new GaussLegendreIntegration(intOrder))); 792 } 793 794 AnalyticHestonEngine::Integration gaussChebyshev(Size intOrder)795 AnalyticHestonEngine::Integration::gaussChebyshev(Size intOrder) { 796 return Integration(GaussChebyshev, 797 ext::shared_ptr<GaussianQuadrature>( 798 new GaussChebyshevIntegration(intOrder))); 799 } 800 801 AnalyticHestonEngine::Integration gaussChebyshev2nd(Size intOrder)802 AnalyticHestonEngine::Integration::gaussChebyshev2nd(Size intOrder) { 803 return Integration(GaussChebyshev2nd, 804 ext::shared_ptr<GaussianQuadrature>( 805 new GaussChebyshev2ndIntegration(intOrder))); 806 } 807 808 AnalyticHestonEngine::Integration discreteSimpson(Size evaluations)809 AnalyticHestonEngine::Integration::discreteSimpson(Size evaluations) { 810 return Integration( 811 DiscreteSimpson, ext::shared_ptr<Integrator>( 812 new DiscreteSimpsonIntegrator(evaluations))); 813 } 814 815 AnalyticHestonEngine::Integration discreteTrapezoid(Size evaluations)816 AnalyticHestonEngine::Integration::discreteTrapezoid(Size evaluations) { 817 return Integration( 818 DiscreteTrapezoid, ext::shared_ptr<Integrator>( 819 new DiscreteTrapezoidIntegrator(evaluations))); 820 } 821 numberOfEvaluations() const822 Size AnalyticHestonEngine::Integration::numberOfEvaluations() const { 823 if (integrator_ != 0) { 824 return integrator_->numberOfEvaluations(); 825 } else if (gaussianQuadrature_ != 0) { 826 return gaussianQuadrature_->order(); 827 } 828 else { 829 QL_FAIL("neither Integrator nor GaussianQuadrature given"); 830 } 831 } 832 isAdaptiveIntegration() const833 bool AnalyticHestonEngine::Integration::isAdaptiveIntegration() const { 834 return intAlgo_ == GaussLobatto 835 || intAlgo_ == GaussKronrod 836 || intAlgo_ == Simpson 837 || intAlgo_ == Trapezoid; 838 } 839 calculate(Real c_inf,const ext::function<Real (Real)> & f,const ext::function<Real ()> & maxBound) const840 Real AnalyticHestonEngine::Integration::calculate( 841 Real c_inf, 842 const ext::function<Real(Real)>& f, 843 const ext::function<Real()>& maxBound) const { 844 845 Real retVal; 846 847 switch(intAlgo_) { 848 case GaussLaguerre: 849 retVal = (*gaussianQuadrature_)(f); 850 break; 851 case GaussLegendre: 852 case GaussChebyshev: 853 case GaussChebyshev2nd: 854 retVal = (*gaussianQuadrature_)(integrand1(c_inf, f)); 855 break; 856 case Simpson: 857 case Trapezoid: 858 case GaussLobatto: 859 case GaussKronrod: 860 if (maxBound != 0 && maxBound() != Null<Real>()) 861 retVal = (*integrator_)(f, 0.0, maxBound()); 862 else 863 retVal = (*integrator_)(integrand2(c_inf, f), 0.0, 1.0); 864 break; 865 case DiscreteTrapezoid: 866 case DiscreteSimpson: 867 if (maxBound != 0 && maxBound() != Null<Real>()) 868 retVal = (*integrator_)(f, 0.0, maxBound()); 869 else 870 retVal = (*integrator_)(integrand3(c_inf, f), 0.0, 1.0); 871 break; 872 default: 873 QL_FAIL("unknwon integration algorithm"); 874 } 875 876 return retVal; 877 } 878 calculate(Real c_inf,const ext::function<Real (Real)> & f,Real maxBound) const879 Real AnalyticHestonEngine::Integration::calculate( 880 Real c_inf, 881 const ext::function<Real(Real)>& f, 882 Real maxBound) const { 883 884 return AnalyticHestonEngine::Integration::calculate( 885 c_inf, f, 886 ext::bind(&constant<Real, Real>::operator(), 887 constant<Real, Real>(maxBound), 1.0)); 888 } 889 andersenPiterbargIntegrationLimit(Real c_inf,Real epsilon,Real v0,Real t)890 Real AnalyticHestonEngine::Integration::andersenPiterbargIntegrationLimit( 891 Real c_inf, Real epsilon, Real v0, Real t) { 892 893 const Real uMaxGuess = -std::log(epsilon)/c_inf; 894 const Real uMaxStep = 0.1*uMaxGuess; 895 896 const Real uMax = Brent().solve(u_Max(c_inf, epsilon), 897 QL_EPSILON*uMaxGuess, uMaxGuess, uMaxStep); 898 899 const Real uHatMax = Brent().solve(uHat_Max(0.5*v0*t, epsilon), 900 QL_EPSILON*std::sqrt(uMaxGuess), 901 std::sqrt(uMaxGuess), 0.1*std::sqrt(uMaxGuess)); 902 903 return std::max(uMax, uHatMax); 904 } 905 } 906