1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2014 Fabien Le Floc'h 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 /*! \file analytichestonengine.hpp 21 \brief analytic Heston expansion engine 22 */ 23 24 #include <ql/math/functional.hpp> 25 #include <ql/pricingengines/blackformula.hpp> 26 27 #include <ql/instruments/payoffs.hpp> 28 #include <ql/pricingengines/vanilla/hestonexpansionengine.hpp> 29 30 #if defined(QL_PATCH_MSVC) 31 #pragma warning(disable: 4180) 32 #endif 33 34 using std::exp; 35 using std::pow; 36 using std::log; 37 using std::sqrt; 38 39 namespace QuantLib { 40 HestonExpansionEngine(const ext::shared_ptr<HestonModel> & model,HestonExpansionFormula formula)41 HestonExpansionEngine::HestonExpansionEngine( 42 const ext::shared_ptr<HestonModel>& model, 43 HestonExpansionFormula formula) 44 : GenericModelEngine<HestonModel, 45 VanillaOption::arguments, 46 VanillaOption::results>(model), 47 formula_(formula) { 48 } 49 calculate() const50 void HestonExpansionEngine::calculate() const 51 { 52 // this is a european option pricer 53 QL_REQUIRE(arguments_.exercise->type() == Exercise::European, 54 "not an European option"); 55 56 // plain vanilla 57 ext::shared_ptr<PlainVanillaPayoff> payoff = 58 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff); 59 QL_REQUIRE(payoff, "non plain vanilla payoff given"); 60 61 const ext::shared_ptr<HestonProcess>& process = model_->process(); 62 63 const Real riskFreeDiscount = process->riskFreeRate()->discount( 64 arguments_.exercise->lastDate()); 65 const Real dividendDiscount = process->dividendYield()->discount( 66 arguments_.exercise->lastDate()); 67 68 const Real spotPrice = process->s0()->value(); 69 QL_REQUIRE(spotPrice > 0.0, "negative or null underlying given"); 70 71 const Real strikePrice = payoff->strike(); 72 const Real term = process->time(arguments_.exercise->lastDate()); 73 74 //possible optimization: 75 // if term=lastTerm & model=lastModel & formula=lastApprox, reuse approx. 76 const Real forward = spotPrice*dividendDiscount/riskFreeDiscount; 77 Real vol; 78 switch(formula_) { 79 case LPP2: { 80 LPP2HestonExpansion expansion(model_->kappa(), model_->theta(), 81 model_->sigma(), model_->v0(), 82 model_->rho(), term); 83 vol = expansion.impliedVolatility(strikePrice, forward); 84 break; 85 } 86 case LPP3: { 87 LPP3HestonExpansion expansion(model_->kappa(), model_->theta(), 88 model_->sigma(), model_->v0(), 89 model_->rho(), term); 90 vol = expansion.impliedVolatility(strikePrice, forward); 91 break; 92 } 93 case Forde: { 94 FordeHestonExpansion expansion(model_->kappa(), model_->theta(), 95 model_->sigma(), model_->v0(), 96 model_->rho(), term); 97 vol = expansion.impliedVolatility(strikePrice, forward); 98 break; 99 } 100 default: 101 QL_FAIL("unknown expansion formula"); 102 } 103 const Real price = blackFormula(payoff, forward, vol*sqrt(term), 104 riskFreeDiscount, 0); 105 results_.value = price; 106 } 107 LPP2HestonExpansion(Real kappa,Real theta,Real sigma,Real v0,Real rho,Real term)108 LPP2HestonExpansion::LPP2HestonExpansion(Real kappa, Real theta, Real sigma, 109 Real v0, Real rho, Real term) { 110 ekt = exp(kappa*term); 111 e2kt = ekt*ekt; 112 e3kt = e2kt*ekt; 113 e4kt = e2kt*e2kt; 114 coeffs[0] = z0(term, kappa, theta, sigma, v0, rho); 115 coeffs[1] = z1(term, kappa, theta, sigma, v0, rho); 116 coeffs[2] = z2(term, kappa, theta, sigma, v0, rho); 117 } 118 fastpow(Real x,int y)119 static Real fastpow(Real x, int y) { 120 return pow(x,y); 121 } 122 impliedVolatility(const Real strike,const Real forward) const123 Real LPP2HestonExpansion::impliedVolatility(const Real strike, 124 const Real forward) const { 125 Real x = log(strike/forward); 126 Real vol = coeffs[0]+x*(coeffs[1]+(x*coeffs[2])); 127 return std::max(1e-8,vol); 128 } 129 z0(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const130 Real LPP2HestonExpansion::z0(Real t, Real kappa, Real theta, 131 Real delta, Real y, Real rho) const { 132 return (4*pow(delta,2)*kappa*(-theta - 4*ekt*(theta + kappa*t*(theta - y)) + 133 e2kt*((5 - 2*kappa*t)*theta - 2*y) + 2*y)* 134 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) + 135 128*ekt*pow(kappa,3)* 136 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,2) + 137 32*delta*ekt*pow(kappa,2)*rho* 138 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 139 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 140 (-1 + ekt - kappa*t)*y) + 141 pow(delta,2)*ekt*pow(rho,2)* 142 (-theta + kappa*t*theta + (theta - y)/ekt + y)* 143 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 144 (-1 + ekt - kappa*t)*y,2) + 145 (48*pow(delta,2)*e2kt*pow(kappa,2)*pow(rho,2)* 146 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 147 (-1 + ekt - kappa*t)*y,2))/ 148 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) - 149 pow(delta,2)*pow(rho,2)*((1 + ekt*(-1 + kappa*t))*theta + 150 (-1 + ekt)*y)*pow((2 + kappa*t + ekt*(-2 + kappa*t))* 151 theta + (-1 + ekt - kappa*t)*y,2) + 152 2*pow(delta,2)*kappa*((1 + ekt*(-1 + kappa*t))*theta + 153 (-1 + ekt)*y)*(theta - 2*y + 154 e2kt*(-5*theta + 2*kappa*t*theta + 2*y + 155 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 156 4*ekt*(theta + kappa*t*theta - kappa*t*y + 157 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y))) - 158 (8*pow(delta,2)*pow(kappa,2)*((1 + ekt*(-1 + kappa*t))*theta + 159 (-1 + ekt)*y)*(theta - 2*y + 160 e2kt*(-5*theta + 2*kappa*t*theta + 2*y + 161 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 162 4*ekt*(theta + kappa*t*theta - kappa*t*y + 163 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)))) 164 /(-theta + kappa*t*theta + (theta - y)/ekt + y))/ 165 (128.*e3kt*pow(kappa,5)*pow(t,2)* 166 pow((-theta + kappa*t*theta + (theta - y)/ekt + y)/(kappa*t),1.5)); 167 } 168 z1(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const169 Real LPP2HestonExpansion::z1(Real t, Real kappa, Real theta, 170 Real delta, Real y, Real rho) const { 171 return (delta*rho*(-(delta*fastpow(-1 + ekt,2)*rho*(4*theta - y)*y) + 172 2*ekt*fastpow(kappa,3)*fastpow(t,2)*theta* 173 ((2 + 2*ekt + delta*rho*t)*theta - (2 + delta*rho*t)*y) - 174 2*(-1 + ekt)*kappa*(2*theta - y)* 175 ((-1 + ekt)*(-2 + delta*rho*t)*theta + 176 (-2 + 2*ekt + delta*rho*t)*y) + 177 fastpow(kappa,2)*t*((-1 + ekt)* 178 (-4 + delta*rho*t + ekt*(-12 + delta*rho*t))*fastpow(theta,2) + 179 2*(-4 + 4*e2kt + delta*rho*t + 3*delta*ekt*rho*t)*theta* 180 y - (-4 + delta*rho*t + 2*ekt*(2 + delta*rho*t))*fastpow(y,2))))/ 181 (8.*fastpow(kappa,2)*t*sqrt((-theta + kappa*t*theta + (theta - y)/ekt + y)/ 182 (kappa*t))*fastpow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y, 183 2)); 184 } 185 z2(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const186 Real LPP2HestonExpansion::z2(Real t, Real kappa, Real theta, 187 Real delta, Real y, Real rho) const { 188 return (fastpow(delta,2)*sqrt((-theta + kappa*t*theta + (theta - y)/ekt + y)/(kappa*t))* 189 (-12*fastpow(rho,2)*fastpow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 190 (-1 + ekt - kappa*t)*y,2) + 191 (-theta + kappa*t*theta + (theta - y)/ekt + y)* 192 (theta - 2*y + e2kt* 193 (-5*theta + 2*kappa*t*theta + 2*y + 8*fastpow(rho,2)*((-3 + kappa*t)*theta + y)) + 194 4*ekt*(theta + kappa*t*theta - kappa*t*y + 195 fastpow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)))) 196 )/(16.*e2kt*fastpow(-theta + kappa*t*theta + (theta - y)/ekt + y, 197 4)); 198 } 199 FordeHestonExpansion(Real kappa,Real theta,Real sigma,Real v0,Real rho,Real term)200 FordeHestonExpansion::FordeHestonExpansion(Real kappa, Real theta, 201 Real sigma, Real v0, 202 Real rho, Real term) { 203 Real v0Sqrt = sqrt(v0); 204 Real rhoBarSquare = 1 - rho * rho; 205 Real sigma00 = v0Sqrt; 206 Real sigma01 = v0Sqrt * (rho * sigma / (4 * v0)); //term in x 207 Real sigma02 = v0Sqrt * ((1 - 5 * rho * rho / 2) / 24 * sigma * sigma/ (v0 * v0)); //term in x*x 208 Real a00 = -sigma * sigma / 12 * (1 - rho * rho / 4) + v0 * rho * sigma / 4 + kappa / 2 * (theta - v0); 209 Real a01 = rho * sigma / (24 * v0) * (sigma * sigma * rhoBarSquare - 2 * kappa * (theta + v0) + v0 * rho * sigma); //term in x 210 Real a02 = (176 * sigma * sigma - 480 * kappa * theta - 712 * rho * rho * sigma * sigma + 521 * rho * rho * rho * rho * sigma * sigma + 40 * sigma * rho * rho * rho * v0 + 1040 * kappa * theta * rho * rho - 80 * v0 * kappa * rho * rho) * sigma * sigma / (v0 * v0 * 7680); 211 coeffs[0] = sigma00*sigma00+a00*term; 212 coeffs[1] = sigma00*sigma01*2+a01*term; 213 coeffs[2] = sigma00*sigma02*2+sigma01*sigma01+a02*term; 214 coeffs[3] = sigma01*sigma02*2; 215 coeffs[4] = sigma02*sigma02; 216 } 217 impliedVolatility(const Real strike,const Real forward) const218 Real FordeHestonExpansion::impliedVolatility(const Real strike, 219 const Real forward) const { 220 Real x = log(strike/forward); 221 Real var = coeffs[0]+x*(coeffs[1]+(x*(coeffs[2]+x*(coeffs[3]+(x*coeffs[4]))))); 222 var = std::max(1e-8,var); 223 return sqrt(var); 224 } 225 z0(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const226 Real LPP3HestonExpansion::z0(Real t, Real kappa, Real theta, 227 Real delta, Real y, Real rho) const { 228 return (96*pow(delta,2)*ekt*pow(kappa,3)* 229 (-theta - 4*ekt*(theta + kappa*t*(theta - y)) + 230 e2kt*((5 - 2*kappa*t)*theta - 2*y) + 2*y)* 231 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) + 232 3072*e2kt*pow(kappa,5)* 233 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,2) + 234 96*pow(delta,3)*ekt*pow(kappa,2)*rho* 235 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 236 (-2*theta - kappa*t*theta - 2*ekt*(2 + kappa*t)* 237 (2*theta + kappa*t*(theta - y)) + e2kt*((10 - 3*kappa*t)*theta - 3*y) + 238 3*y + 2*kappa*t*y) + 768*delta*e2kt*pow(kappa,4)*rho* 239 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 240 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 241 (-1 + ekt - kappa*t)*y) + 242 6*pow(delta,3)*kappa*rho*(-theta - 4*ekt*(theta + kappa*t*(theta - y)) + 243 e2kt*((5 - 2*kappa*t)*theta - 2*y) + 2*y)* 244 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 245 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 246 (-1 + ekt - kappa*t)*y) + 247 24*pow(delta,2)*e2kt*pow(kappa,2)*pow(rho,2)* 248 (-theta + kappa*t*theta + (theta - y)/ekt + y)* 249 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 250 (-1 + ekt - kappa*t)*y,2) + 251 (1152*pow(delta,2)*e3kt*pow(kappa,4)*pow(rho,2)* 252 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 253 (-1 + ekt - kappa*t)*y,2))/ 254 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) - 255 24*pow(delta,2)*ekt*pow(kappa,2)*pow(rho,2)* 256 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 257 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 258 (-1 + ekt - kappa*t)*y,2) + 259 80*pow(delta,3)*ekt*kappa*pow(rho,3)* 260 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 261 (-1 + ekt - kappa*t)*y,3) + 262 pow(delta,3)*ekt*pow(rho,3)* 263 (-theta + kappa*t*theta + (theta - y)/ekt + y)* 264 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 265 (-1 + ekt - kappa*t)*y,3) - 266 (1440*pow(delta,3)*e3kt*pow(kappa,3)*pow(rho,3)* 267 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 268 (-1 + ekt - kappa*t)*y,3))/ 269 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,2) - 270 (528*pow(delta,3)*e2kt*pow(kappa,2)*pow(rho,3)* 271 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 272 (-1 + ekt - kappa*t)*y,3))/ 273 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) - 274 3*pow(delta,3)*pow(rho,3)*((1 + ekt*(-1 + kappa*t))*theta + 275 (-1 + ekt)*y)*pow((2 + kappa*t + ekt*(-2 + kappa*t))* 276 theta + (-1 + ekt - kappa*t)*y,3) + 277 384*pow(delta,3)*e2kt*pow(kappa,3)*rho* 278 ((2 + kappa*t + 2*ekt*pow(2 + kappa*t,2) + 279 e2kt*(-10 + 3*kappa*t))*theta + 280 (-3 + 3*e2kt - 2*kappa*t - 2*ekt*kappa*t*(2 + kappa*t))*y) - 281 (576*pow(delta,3)*e2kt*pow(kappa,3)*rho* 282 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 283 (-1 + ekt - kappa*t)*y)* 284 ((1 + e2kt*(-5 + 2*kappa*t + 4*pow(rho,2)*(-3 + kappa*t)) + 285 2*ekt*(2 + 2*kappa*t + 286 pow(rho,2)*(6 + 4*kappa*t + pow(kappa,2)*pow(t,2))))*theta + 287 2*(-1 + e2kt*(1 + 2*pow(rho,2)) - 288 ekt*(2*kappa*t + 289 pow(rho,2)*(2 + 2*kappa*t + pow(kappa,2)*pow(t,2))))*y))/ 290 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) + 291 pow(delta,3)*rho*((1 + ekt*(-1 + kappa*t))*theta + 292 (-1 + ekt)*y)*((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 293 (-1 + ekt - kappa*t)*y)* 294 (theta*(12*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2) + 295 8*pow(-1 + ekt,2)*pow(rho,2)*theta - 296 (-1 + ekt)*kappa* 297 (3 + 8*pow(rho,2)*t*theta + ekt*(15 + 8*pow(rho,2)*(9 + t*theta))) 298 + 2*pow(kappa,2)*t*(pow(rho,2)*t*theta + 299 2*ekt*(3 + pow(rho,2)*(12 + t*theta)) + 300 e2kt*(3 + pow(rho,2)*(12 + t*theta)))) - 301 2*(6*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2) + 302 4*pow(-1 + ekt,2)*pow(rho,2)*theta + 303 2*pow(kappa,2)*t*(pow(rho,2)*t*theta + 304 ekt*(3 + pow(rho,2)*(6 + t*theta))) - 305 (-1 + ekt)*kappa* 306 (3 + 6*pow(rho,2)*t*theta + ekt*(3 + 2*pow(rho,2)*(6 + t*theta))))* 307 y + 2*pow(rho,2)*pow(1 - ekt + kappa*t,2)*pow(y,2)) - 308 (40*pow(delta,3)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 309 (-1 + ekt)*y)*((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 310 (-1 + ekt - kappa*t)*y)* 311 (theta*(12*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2) + 312 8*pow(-1 + ekt,2)*pow(rho,2)*theta - 313 (-1 + ekt)*kappa* 314 (3 + 8*pow(rho,2)*t*theta + 315 ekt*(15 + 8*pow(rho,2)*(9 + t*theta))) + 316 2*pow(kappa,2)*t*(pow(rho,2)*t*theta + 317 2*ekt*(3 + pow(rho,2)*(12 + t*theta)) + 318 e2kt*(3 + pow(rho,2)*(12 + t*theta)))) - 319 2*(6*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2) + 320 4*pow(-1 + ekt,2)*pow(rho,2)*theta + 321 2*pow(kappa,2)*t*(pow(rho,2)*t*theta + 322 ekt*(3 + pow(rho,2)*(6 + t*theta))) - 323 (-1 + ekt)*kappa* 324 (3 + 6*pow(rho,2)*t*theta + ekt*(3 + 2*pow(rho,2)*(6 + t*theta))) 325 )*y + 2*pow(rho,2)*pow(1 - ekt + kappa*t,2)*pow(y,2)))/ 326 (-theta + kappa*t*theta + (theta - y)/ekt + y) - 327 12*pow(delta,3)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 328 (-1 + ekt)*y)*(2*theta + kappa*t*theta - y - kappa*t*y + 329 ekt*((-2 + kappa*t)*theta + y))* 330 (theta - 2*y + e2kt* 331 (-5*theta + 2*kappa*t*theta + 2*y + 4*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 332 2*ekt*(2*(theta + kappa*t*(theta - y)) + 333 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y))) + 334 (288*pow(delta,3)*pow(kappa,2)*rho* 335 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 336 (2*theta + kappa*t*theta - y - kappa*t*y + ekt*((-2 + kappa*t)*theta + y))* 337 (theta - 2*y + e2kt* 338 (-5*theta + 2*kappa*t*theta + 2*y + 4*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 339 2*ekt*(2*(theta + kappa*t*(theta - y)) + 340 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)))) 341 /(-theta + kappa*t*theta + (theta - y)/ekt + y) + 342 48*pow(delta,2)*ekt*pow(kappa,3)* 343 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 344 (theta - 2*y + e2kt* 345 (-5*theta + 2*kappa*t*theta + 2*y + 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 346 4*ekt*(theta + kappa*t*theta - kappa*t*y + 347 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y))) - 348 (192*pow(delta,2)*ekt*pow(kappa,4)* 349 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 350 (theta - 2*y + e2kt* 351 (-5*theta + 2*kappa*t*theta + 2*y + 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 352 4*ekt*(theta + kappa*t*theta - kappa*t*y + 353 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)))) 354 /(-theta + kappa*t*theta + (theta - y)/ekt + y) + 355 3*pow(delta,3)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 356 (-1 + ekt)*y)*((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 357 (-1 + ekt - kappa*t)*y)* 358 (theta - 2*y + e2kt* 359 (-5*theta + 2*kappa*t*theta + 2*y + 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 360 4*ekt*(theta + kappa*t*theta - kappa*t*y + 361 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y))) - 362 (12*pow(delta,3)*pow(kappa,2)*rho* 363 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 364 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 365 (-1 + ekt - kappa*t)*y)* 366 (theta - 2*y + e2kt* 367 (-5*theta + 2*kappa*t*theta + 2*y + 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 368 4*ekt*(theta + kappa*t*theta - kappa*t*y + 369 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)))) 370 /(-theta + kappa*t*theta + (theta - y)/ekt + y) + 371 4*pow(delta,3)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 372 (-1 + ekt)*y)*(3*(theta - 2*y)*((2 + kappa*t)*theta - (1 + kappa*t)*y) + 373 3*ekt*(6*pow(theta,2) + theta*y - 2*pow(y,2) + 374 kappa*(13*t*pow(theta,2) + theta*(8 - 18*t*y) + 4*y*(-3 + t*y)) + 375 4*pow(kappa,2)*t*(theta + t*pow(theta,2) - 2*t*theta*y + y*(-2 + t*y))) + 376 3*e3kt*(10*pow(theta,2) + 377 2*pow(kappa,2)*t*theta*(6 + 8*pow(rho,2) + t*theta) - 9*theta*y + 2*pow(y,2) + 378 kappa*(-9*t*pow(theta,2) + 4*(3 + 4*pow(rho,2))*y + 379 theta*(-40 - 64*pow(rho,2) + 4*t*y))) + 380 e2kt*(-54*pow(theta,2) + 381 8*pow(kappa,4)*pow(rho,2)*pow(t,3)*(theta - y) + 39*theta*y - 6*pow(y,2) + 382 24*pow(kappa,3)*pow(t,2)*(theta + 2*pow(rho,2)*theta - (1 + pow(rho,2))*y) + 383 6*pow(kappa,2)*t*(3*t*pow(theta,2) - 8*(1 + pow(rho,2))*y + 384 theta*(16 + 24*pow(rho,2) - 3*t*y)) - 385 3*kappa*(5*t*pow(theta,2) + 2*y*(8*pow(rho,2) + 3*t*y) - 386 theta*(32 + 64*pow(rho,2) + 17*t*y)))) - 387 (48*pow(delta,3)*pow(kappa,2)*rho* 388 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)* 389 (3*(theta - 2*y)*((2 + kappa*t)*theta - (1 + kappa*t)*y) + 390 3*ekt*(6*pow(theta,2) + theta*y - 2*pow(y,2) + 391 kappa*(13*t*pow(theta,2) + theta*(8 - 18*t*y) + 4*y*(-3 + t*y)) + 392 4*pow(kappa,2)*t*(theta + t*pow(theta,2) - 2*t*theta*y + y*(-2 + t*y))) + 393 3*e3kt*(10*pow(theta,2) + 394 2*pow(kappa,2)*t*theta*(6 + 8*pow(rho,2) + t*theta) - 9*theta*y + 395 2*pow(y,2) + kappa*(-9*t*pow(theta,2) + 4*(3 + 4*pow(rho,2))*y + 396 theta*(-40 - 64*pow(rho,2) + 4*t*y))) + 397 e2kt*(-54*pow(theta,2) + 398 8*pow(kappa,4)*pow(rho,2)*pow(t,3)*(theta - y) + 39*theta*y - 6*pow(y,2) + 399 24*pow(kappa,3)*pow(t,2)* 400 (theta + 2*pow(rho,2)*theta - (1 + pow(rho,2))*y) + 401 6*pow(kappa,2)*t*(3*t*pow(theta,2) - 8*(1 + pow(rho,2))*y + 402 theta*(16 + 24*pow(rho,2) - 3*t*y)) - 403 3*kappa*(5*t*pow(theta,2) + 2*y*(8*pow(rho,2) + 3*t*y) - 404 theta*(32 + 64*pow(rho,2) + 17*t*y)))))/ 405 (-theta + kappa*t*theta + (theta - y)/ekt + y) + 406 (240*pow(delta,3)*e2kt*pow(kappa,2)*rho* 407 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 408 (-1 + ekt - kappa*t)*y)* 409 (12*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2)*(theta - y) + 410 2*pow(-1 + ekt,2)*pow(rho,2)*pow(-2*theta + y,2) - 411 (-1 + ekt)*kappa* 412 (8*(1 + ekt)*pow(rho,2)*t*pow(theta,2) + 413 2*y*(-3 - 3*ekt*(1 + 4*pow(rho,2)) + 2*pow(rho,2)*t*y) + 414 theta*(3 - 12*pow(rho,2)*t*y + ekt*(15 + pow(rho,2)*(72 - 4*t*y))) 415 ) + 2*pow(kappa,2)*t*(e2kt*theta* 416 (3 + pow(rho,2)*(12 + t*theta)) + pow(rho,2)*t*pow(theta - y,2) + 417 2*ekt*(pow(rho,2)*t*pow(theta,2) - 3*(y + 2*pow(rho,2)*y) + 418 theta*(3 + pow(rho,2)*(12 - t*y))))))/ 419 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y))/ 420 (3072.*e4kt*pow(kappa,7)*pow(t,2)* 421 pow((-theta + kappa*t*theta + (theta - y)/ekt + y)/(kappa*t),1.5)); 422 } 423 z1(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const424 Real LPP3HestonExpansion::z1(Real t, Real kappa, Real theta, 425 Real delta, Real y, Real rho) const { 426 return (delta*(768*e2kt*pow(kappa,4)*rho* 427 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 428 (-1 + ekt - kappa*t)*y) - 429 (576*delta*e2kt*pow(kappa,3)*pow(rho,2)* 430 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 431 (-1 + ekt - kappa*t)*y,2))/ 432 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) - 433 10*pow(delta,2)*pow(rho,3)*pow((2 + kappa*t + ekt*(-2 + kappa*t))* 434 theta + (-1 + ekt - kappa*t)*y,3) + 435 (6*pow(delta,2)*kappa*pow(rho,3)* 436 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 437 (-1 + ekt - kappa*t)*y,3))/ 438 (-theta + kappa*t*theta + (theta - y)/ekt + y) - 439 (3360*pow(delta,2)*e3kt*pow(kappa,3)*pow(rho,3)* 440 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 441 (-1 + ekt - kappa*t)*y,3))/ 442 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,3) - 443 (288*pow(delta,2)*e2kt*pow(kappa,2)*pow(rho,3)* 444 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 445 (-1 + ekt - kappa*t)*y,3))/ 446 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,2) + 447 (234*pow(delta,2)*ekt*kappa*pow(rho,3)* 448 pow((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 449 (-1 + ekt - kappa*t)*y,3))/ 450 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) - 451 96*delta*ekt*pow(kappa,3)* 452 ((1 + 4*ekt*(1 + kappa*t) + e2kt*(-5 + 2*kappa*t))*theta + 453 2*(-1 + e2kt - 2*ekt*kappa*t)*y) - 454 12*pow(delta,2)*kappa*rho*((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 455 (-1 + ekt - kappa*t)*y)* 456 ((1 + 4*ekt*(1 + kappa*t) + e2kt*(-5 + 2*kappa*t))*theta + 457 2*(-1 + e2kt - 2*ekt*kappa*t)*y) - 458 192*pow(delta,2)*ekt*pow(kappa,2)*rho* 459 ((2 + kappa*t + 2*ekt*pow(2 + kappa*t,2) + 460 e2kt*(-10 + 3*kappa*t))*theta + 461 (-3 + 3*e2kt - 2*kappa*t - 2*ekt*kappa*t*(2 + kappa*t))*y) 462 - (12*pow(delta,2)*ekt*pow(kappa,2)*rho* 463 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 464 (-1 + ekt - kappa*t)*y)* 465 ((1 + e2kt*(-5 + 2*kappa*t + 8*pow(rho,2)*(-3 + kappa*t)) + 466 4*ekt*(1 + kappa*t + 467 pow(rho,2)*(6 + 4*kappa*t + pow(kappa,2)*pow(t,2))))*theta + 468 2*(-1 + e2kt*(1 + 4*pow(rho,2)) - 469 2*ekt*(kappa*t + 470 pow(rho,2)*(2 + 2*kappa*t + pow(kappa,2)*pow(t,2))))*y))/ 471 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) + 472 (576*pow(delta,2)*ekt*pow(kappa,2)*rho* 473 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 474 (-1 + ekt - kappa*t)*y)* 475 ((1 + e2kt*(-5 + 2*kappa*t + 4*pow(rho,2)*(-3 + kappa*t)) + 476 2*ekt*(2 + 2*kappa*t + 477 pow(rho,2)*(6 + 4*kappa*t + pow(kappa,2)*pow(t,2))))*theta + 478 2*(-1 + e2kt*(1 + 2*pow(rho,2)) - 479 ekt*(2*kappa*t + 480 pow(rho,2)*(2 + 2*kappa*t + pow(kappa,2)*pow(t,2))))*y))/ 481 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) + 482 (5*pow(delta,2)*rho*((1 + ekt*(-1 + kappa*t))*theta + 483 (-1 + ekt)*y)* 484 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 485 (-1 + ekt - kappa*t)*y)* 486 (theta*(12*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2) + 487 8*pow(-1 + ekt,2)*pow(rho,2)*theta - 488 (-1 + ekt)*kappa* 489 (3 + 8*pow(rho,2)*t*theta + 490 ekt*(15 + 8*pow(rho,2)*(9 + t*theta))) + 491 2*pow(kappa,2)*t*(pow(rho,2)*t*theta + 492 2*ekt*(3 + pow(rho,2)*(12 + t*theta)) + 493 e2kt*(3 + pow(rho,2)*(12 + t*theta)))) - 494 2*(6*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2) + 495 4*pow(-1 + ekt,2)*pow(rho,2)*theta + 496 2*pow(kappa,2)*t*(pow(rho,2)*t*theta + 497 ekt*(3 + pow(rho,2)*(6 + t*theta))) - 498 (-1 + ekt)*kappa* 499 (3 + 6*pow(rho,2)*t*theta + 500 ekt*(3 + 2*pow(rho,2)*(6 + t*theta))))*y + 501 2*pow(rho,2)*pow(1 - ekt + kappa*t,2)*pow(y,2)))/ 502 (ekt*(-theta + kappa*t*theta + (theta - y)/ekt + y)) - 503 (48*pow(delta,2)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 504 (-1 + ekt)*y)* 505 (2*theta + kappa*t*theta - y - kappa*t*y + 506 ekt*((-2 + kappa*t)*theta + y))* 507 (theta - 2*y + e2kt* 508 (-5*theta + 2*kappa*t*theta + 2*y + 4*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 509 2*ekt*(2*(theta + kappa*t*(theta - y)) + 510 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)) 511 ))/(ekt*(-theta + kappa*t*theta + (theta - y)/ekt + y)) + 512 (96*delta*pow(kappa,3)*((1 + ekt*(-1 + kappa*t))*theta + 513 (-1 + ekt)*y)* 514 (theta - 2*y + e2kt* 515 (-5*theta + 2*kappa*t*theta + 2*y + 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 516 4*ekt*(theta + kappa*t*theta - kappa*t*y + 517 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)) 518 ))/(-theta + kappa*t*theta + (theta - y)/ekt + y) + 519 (9*pow(delta,2)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 520 (-1 + ekt)*y)* 521 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 522 (-1 + ekt - kappa*t)*y)* 523 (theta - 2*y + e2kt* 524 (-5*theta + 2*kappa*t*theta + 2*y + 8*pow(rho,2)*((-3 + kappa*t)*theta + y)) + 525 4*ekt*(theta + kappa*t*theta - kappa*t*y + 526 pow(rho,2)*((6 + kappa*t*(4 + kappa*t))*theta - (2 + kappa*t*(2 + kappa*t))*y)) 527 ))/(ekt*(-theta + kappa*t*theta + (theta - y)/ekt + y)) - 528 (48*pow(delta,2)*ekt*pow(kappa,2)*rho* 529 (3*(theta - 2*y)*((2 + kappa*t)*theta - (1 + kappa*t)*y) + 530 3*ekt*(6*pow(theta,2) + theta*y - 2*pow(y,2) + 531 kappa*(13*t*pow(theta,2) + theta*(8 - 18*t*y) + 4*y*(-3 + t*y)) + 532 4*pow(kappa,2)*t*(theta + t*pow(theta,2) - 2*t*theta*y + y*(-2 + t*y))) + 533 3*e3kt*(10*pow(theta,2) + 534 2*pow(kappa,2)*t*theta*(6 + 8*pow(rho,2) + t*theta) - 9*theta*y + 535 2*pow(y,2) + kappa*(-9*t*pow(theta,2) + 4*(3 + 4*pow(rho,2))*y + 536 theta*(-40 - 64*pow(rho,2) + 4*t*y))) + 537 e2kt*(-54*pow(theta,2) + 538 8*pow(kappa,4)*pow(rho,2)*pow(t,3)*(theta - y) + 39*theta*y - 539 6*pow(y,2) + 24*pow(kappa,3)*pow(t,2)* 540 (theta + 2*pow(rho,2)*theta - (1 + pow(rho,2))*y) + 541 6*pow(kappa,2)*t*(3*t*pow(theta,2) - 8*(1 + pow(rho,2))*y + 542 theta*(16 + 24*pow(rho,2) - 3*t*y)) - 543 3*kappa*(5*t*pow(theta,2) + 2*y*(8*pow(rho,2) + 3*t*y) - 544 theta*(32 + 64*pow(rho,2) + 17*t*y)))))/ 545 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y) + 546 (12*pow(delta,2)*kappa*rho*((1 + ekt*(-1 + kappa*t))*theta + 547 (-1 + ekt)*y)* 548 (3*(theta - 2*y)*((2 + kappa*t)*theta - (1 + kappa*t)*y) + 549 3*ekt*(6*pow(theta,2) + theta*y - 2*pow(y,2) + 550 kappa*(13*t*pow(theta,2) + theta*(8 - 18*t*y) + 4*y*(-3 + t*y)) + 551 4*pow(kappa,2)*t*(theta + t*pow(theta,2) - 2*t*theta*y + y*(-2 + t*y))) + 552 3*e3kt*(10*pow(theta,2) + 553 2*pow(kappa,2)*t*theta*(6 + 8*pow(rho,2) + t*theta) - 9*theta*y + 554 2*pow(y,2) + kappa*(-9*t*pow(theta,2) + 4*(3 + 4*pow(rho,2))*y + 555 theta*(-40 - 64*pow(rho,2) + 4*t*y))) + 556 e2kt*(-54*pow(theta,2) + 557 8*pow(kappa,4)*pow(rho,2)*pow(t,3)*(theta - y) + 39*theta*y - 558 6*pow(y,2) + 24*pow(kappa,3)*pow(t,2)* 559 (theta + 2*pow(rho,2)*theta - (1 + pow(rho,2))*y) + 560 6*pow(kappa,2)*t*(3*t*pow(theta,2) - 8*(1 + pow(rho,2))*y + 561 theta*(16 + 24*pow(rho,2) - 3*t*y)) - 562 3*kappa*(5*t*pow(theta,2) + 2*y*(8*pow(rho,2) + 3*t*y) - 563 theta*(32 + 64*pow(rho,2) + 17*t*y)))))/ 564 (ekt*(-theta + kappa*t*theta + (theta - y)/ekt + y)) + 565 (240*pow(delta,2)*e2kt*pow(kappa,2)*rho* 566 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 567 (-1 + ekt - kappa*t)*y)* 568 (12*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2)*(theta - y) + 569 2*pow(-1 + ekt,2)*pow(rho,2)*pow(-2*theta + y,2) - 570 (-1 + ekt)*kappa* 571 (8*(1 + ekt)*pow(rho,2)*t*pow(theta,2) + 572 2*y*(-3 - 3*ekt*(1 + 4*pow(rho,2)) + 2*pow(rho,2)*t*y) + 573 theta*(3 - 12*pow(rho,2)*t*y + 574 ekt*(15 + pow(rho,2)*(72 - 4*t*y)))) + 575 2*pow(kappa,2)*t*(e2kt*theta*(3 + pow(rho,2)*(12 + t*theta)) + 576 pow(rho,2)*t*pow(theta - y,2) + 577 2*ekt*(pow(rho,2)*t*pow(theta,2) - 3*(y + 2*pow(rho,2)*y) + 578 theta*(3 + pow(rho,2)*(12 - t*y))))))/ 579 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,2) - 580 (120*pow(delta,2)*ekt*kappa*rho* 581 ((2 + kappa*t + ekt*(-2 + kappa*t))*theta + 582 (-1 + ekt - kappa*t)*y)* 583 (12*ekt*pow(kappa,3)*pow(rho,2)*pow(t,2)*(theta - y) + 584 2*pow(-1 + ekt,2)*pow(rho,2)*pow(-2*theta + y,2) - 585 (-1 + ekt)*kappa* 586 (8*(1 + ekt)*pow(rho,2)*t*pow(theta,2) + 587 2*y*(-3 - 3*ekt*(1 + 4*pow(rho,2)) + 2*pow(rho,2)*t*y) + 588 theta*(3 - 12*pow(rho,2)*t*y + 589 ekt*(15 + pow(rho,2)*(72 - 4*t*y)))) + 590 2*pow(kappa,2)*t*(e2kt*theta*(3 + pow(rho,2)*(12 + t*theta)) + 591 pow(rho,2)*t*pow(theta - y,2) + 592 2*ekt*(pow(rho,2)*t*pow(theta,2) - 3*(y + 2*pow(rho,2)*y) + 593 theta*(3 + pow(rho,2)*(12 - t*y))))))/ 594 ((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y)))/ 595 (1536.*e3kt*pow(kappa,6)*pow(t,2)* 596 pow((-theta + kappa*t*theta + (theta - y)/ekt + y)/(kappa*t),1.5)); 597 } 598 z2(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const599 Real LPP3HestonExpansion::z2(Real t, Real kappa, Real theta, 600 Real delta, Real y, Real rho) const{ 601 return (pow(delta,2)*(8*e3kt*pow(kappa,5)*pow(rho,2)*pow(t,4)*(2 + delta*rho*t)* 602 pow(theta,2)*(theta - y) - delta*pow(-1 + ekt,3)*rho* 603 (2*(-1 + ekt*(-5 + 24*pow(rho,2)))*pow(theta,3) + 604 (7 + ekt*(3 + 56*pow(rho,2)))*pow(theta,2)*y - 605 3*(1 + ekt*(-3 + 8*pow(rho,2)))*theta*pow(y,2) + 606 2*(-1 + ekt*(-1 + 2*pow(rho,2)))*pow(y,3)) - 607 pow(-1 + ekt,2)*kappa* 608 ((-4 + delta*rho*t - 8*ekt* 609 (2 - 12*pow(rho,2) - 4*delta*rho*t + 25*delta*pow(rho,3)*t) + 610 e2kt*(20 - 96*pow(rho,2) + 3*delta*rho*t + 56*delta*pow(rho,3)*t) 611 )*pow(theta,3) - 2*(-8 + 2*delta*rho*t + 612 e2kt*(24 - 80*pow(rho,2) - 9*delta*rho*t + 613 24*delta*pow(rho,3)*t) - 614 4*ekt*(4 - 20*pow(rho,2) - 10*delta*rho*t + 39*delta*pow(rho,3)*t) 615 )*pow(theta,2)*y + (5*(-4 + delta*rho*t) + 616 ekt*(-16 + 80*pow(rho,2) + 57*delta*rho*t - 617 140*delta*pow(rho,3)*t) + 618 2*e2kt*(18 - 40*pow(rho,2) - 3*delta*rho*t + 619 6*delta*pow(rho,3)*t))*theta*pow(y,2) + 620 2*(4 + e2kt*(-4 + 8*pow(rho,2)) - delta*rho*t + 621 ekt*rho*(-8*rho - 7*delta*t + 14*delta*pow(rho,2)*t))*pow(y,3)) + 622 ekt*(-1 + ekt)*pow(kappa,2)*t* 623 ((-24 + 128*pow(rho,2) + 9*delta*rho*t - 144*delta*pow(rho,3)*t - 624 4*ekt*(6 - 8*pow(rho,2) - 9*delta*rho*t + 6*delta*pow(rho,3)*t) + 625 e2kt*(48 - 160*pow(rho,2) - 9*delta*rho*t + 626 24*delta*pow(rho,3)*t))*pow(theta,3) - 627 (-72 + 320*pow(rho,2) + 27*delta*rho*t - 360*delta*pow(rho,3)*t - 628 ekt*rho*(160*rho - 81*delta*t + 348*delta*pow(rho,2)*t) + 629 2*e2kt*(36 - 80*pow(rho,2) - 3*delta*rho*t + 630 6*delta*pow(rho,3)*t))*pow(theta,2)*y - 631 2*(32 - 128*pow(rho,2) + 12*e2kt*(-1 + 2*pow(rho,2)) - 632 15*delta*rho*t + 144*delta*pow(rho,3)*t + 633 2*ekt*(-10 + 52*pow(rho,2) - 13*delta*rho*t + 634 58*delta*pow(rho,3)*t))*theta*pow(y,2) + 635 4*(4 - 16*pow(rho,2) - 3*delta*rho*t + 18*delta*pow(rho,3)*t + 636 ekt*(-4 + 16*pow(rho,2) - 2*delta*rho*t + 11*delta*pow(rho,3)*t))* 637 pow(y,3)) - 4*e2kt*pow(kappa,4)*pow(t,3)*theta* 638 (2*e2kt*(-1 + 2*pow(rho,2))*pow(theta,2) + 639 pow(rho,2)*(4 + 13*delta*rho*t)*pow(theta - y,2) + 640 ekt*((-4 + 16*pow(rho,2) - 2*delta*rho*t + 9*delta*pow(rho,3)*t)* 641 pow(theta,2) + (4 - 32*pow(rho,2) + 2*delta*rho*t - 19*delta*pow(rho,3)*t)* 642 theta*y + 4*pow(rho,2)*(2 + delta*rho*t)*pow(y,2))) - 643 2*ekt*pow(kappa,3)*pow(t,2)* 644 (-4*pow(rho,2)*(-4 + 3*delta*rho*t)*pow(theta - y,3) + 645 e3kt*pow(theta,2)* 646 ((18 - 40*pow(rho,2) - delta*rho*t + 2*delta*pow(rho,3)*t)*theta + 647 12*(-1 + 2*pow(rho,2))*y) + 648 2*ekt*((-9 + 36*pow(rho,2) + 19*delta*pow(rho,3)*t)*pow(theta,3) + 649 2*(9 - 30*pow(rho,2) + 7*delta*pow(rho,3)*t)*pow(theta,2)*y + 650 (-8 + 20*pow(rho,2) + delta*rho*t - 46*delta*pow(rho,3)*t)*theta*pow(y,2) + 651 pow(rho,2)*(4 + 13*delta*rho*t)*pow(y,3)) + 652 e2kt*(8*theta*y*(-3*theta + 2*y) + 653 delta*rho*t*theta*(7*pow(theta,2) - 23*theta*y + 8*pow(y,2)) - 654 8*pow(rho,2)*(6*pow(theta,3) - 18*pow(theta,2)*y + 11*theta*pow(y,2) - 655 pow(y,3)) + 4*delta*pow(rho,3)*t* 656 (-13*pow(theta,3) + 31*pow(theta,2)*y - 14*theta*pow(y,2) + pow(y,3))))))/ 657 (64.*pow(kappa,2)*t*sqrt((-theta + kappa*t*theta + (theta - y)/ekt + y)/ 658 (kappa*t))*pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y, 659 4)); 660 } 661 z3(Real t,Real kappa,Real theta,Real delta,Real y,Real rho) const662 Real LPP3HestonExpansion::z3(Real t, Real kappa, Real theta, 663 Real delta, Real y, Real rho) const { 664 return (pow(delta,3)*ekt*rho*((-15*(2 + kappa*t) + 665 3*e4kt*(50 - 79*kappa*t + 35*pow(kappa,2)*pow(t,2) - 666 6*pow(kappa,3)*pow(t,3) + 667 8*pow(rho,2)*(-18 + 15*kappa*t - 6*pow(kappa,2)*pow(t,2) + 668 pow(kappa,3)*pow(t,3))) + 669 ekt*(-3*(20 + 86*kappa*t + 29*pow(kappa,2)*pow(t,2)) + 670 pow(rho,2)*(432 + 936*kappa*t + 552*pow(kappa,2)*pow(t,2) + 671 92*pow(kappa,3)*pow(t,3))) + 672 e2kt*(360 + 324*kappa*t - 261*pow(kappa,2)*pow(t,2) - 673 48*pow(kappa,3)*pow(t,3) - 674 4*pow(rho,2)*(324 + 378*kappa*t - 12*pow(kappa,2)*pow(t,2) - 675 2*pow(kappa,3)*pow(t,3) + 23*pow(kappa,4)*pow(t,4))) + 676 e3kt*(3*(-140 + 62*kappa*t + 81*pow(kappa,2)*pow(t,2) - 677 38*pow(kappa,3)*pow(t,3) + 8*pow(kappa,4)*pow(t,4)) + 678 4*pow(rho,2)*(324 + 54*kappa*t - 114*pow(kappa,2)*pow(t,2) + 679 77*pow(kappa,3)*pow(t,3) - 19*pow(kappa,4)*pow(t,4) + 680 2*pow(kappa,5)*pow(t,5))))*pow(theta,3) + 681 (15*(7 + 4*kappa*t) + 3*e4kt* 682 (-79 + 70*kappa*t - 18*pow(kappa,2)*pow(t,2) + 683 24*pow(rho,2)*(5 - 4*kappa*t + pow(kappa,2)*pow(t,2))) - 684 3*ekt*(26 - 200*kappa*t - 87*pow(kappa,2)*pow(t,2) + 685 4*pow(rho,2)*(30 + 142*kappa*t + 115*pow(kappa,2)*pow(t,2) + 686 23*pow(kappa,3)*pow(t,3))) + 687 2*e2kt*(3*(-66 - 195*kappa*t + 63*pow(kappa,2)*pow(t,2) + 688 16*pow(kappa,3)*pow(t,3)) + 689 4*pow(rho,2)*(135 + 390*kappa*t - 9*pow(kappa,2)*pow(t,2) - 690 48*pow(kappa,3)*pow(t,3) + 23*pow(kappa,4)*pow(t,4))) + 691 e3kt*(606 + 300*kappa*t - 585*pow(kappa,2)*pow(t,2) + 692 210*pow(kappa,3)*pow(t,3) - 24*pow(kappa,4)*pow(t,4) - 693 4*pow(rho,2)*(270 + 282*kappa*t - 345*pow(kappa,2)*pow(t,2) + 694 153*pow(kappa,3)*pow(t,3) - 29*pow(kappa,4)*pow(t,4) + 695 2*pow(kappa,5)*pow(t,5))))*pow(theta,2)*y + 696 (-93 - 75*kappa*t + 3*e4kt* 697 (35 - 18*kappa*t + 24*pow(rho,2)*(-2 + kappa*t)) + 698 3*ekt*(58 - 123*kappa*t - 86*pow(kappa,2)*pow(t,2) + 699 4*pow(rho,2)*(12 + 80*kappa*t + 92*pow(kappa,2)*pow(t,2) + 700 23*pow(kappa,3)*pow(t,3))) + 701 e3kt*(-3*(74 + 137*kappa*t - 100*pow(kappa,2)*pow(t,2) + 702 16*pow(kappa,3)*pow(t,3)) - 703 16*pow(rho,2)*(-27 - 51*kappa*t + 45*pow(kappa,2)*pow(t,2) - 704 12*pow(kappa,3)*pow(t,3) + pow(kappa,4)*pow(t,4))) + 705 e2kt*(36 + 909*kappa*t - 42*pow(kappa,2)*pow(t,2) - 706 60*pow(kappa,3)*pow(t,3) - 707 4*pow(rho,2)*(108 + 462*kappa*t + 96*pow(kappa,2)*pow(t,2) - 708 117*pow(kappa,3)*pow(t,3) + 23*pow(kappa,4)*pow(t,4))))*theta*pow(y,2) 709 + 2*(9 + 3*e4kt*(-3 + 4*pow(rho,2)) + 15*kappa*t + 710 e2kt*(-3*kappa*t*(33 + 10*kappa*t) + 711 pow(rho,2)*(36 + 192*kappa*t + 96*pow(kappa,2)*pow(t,2) - 712 46*pow(kappa,3)*pow(t,3))) + 713 e3kt*(18 + 57*kappa*t - 12*pow(kappa,2)*pow(t,2) - 714 2*pow(rho,2)*(18 + 48*kappa*t - 21*pow(kappa,2)*pow(t,2) + 715 2*pow(kappa,3)*pow(t,3))) + 716 ekt*(3*(-6 + 9*kappa*t + 14*pow(kappa,2)*pow(t,2)) - 717 2*pow(rho,2)*(6 + 48*kappa*t + 69*pow(kappa,2)*pow(t,2) + 718 23*pow(kappa,3)*pow(t,3))))*pow(y,3)))/ 719 (96.*kappa*t*sqrt((-theta + kappa*t*theta + (theta - y)/ekt + y)/(kappa*t))* 720 pow((1 + ekt*(-1 + kappa*t))*theta + (-1 + ekt)*y,5)); 721 } 722 LPP3HestonExpansion(Real kappa,Real theta,Real sigma,Real v0,Real rho,Real term)723 LPP3HestonExpansion::LPP3HestonExpansion(Real kappa, Real theta, Real sigma, 724 Real v0, Real rho, Real term) { 725 ekt = exp(kappa*term); 726 e2kt = ekt*ekt; 727 e3kt = e2kt*ekt; 728 e4kt = e2kt*e2kt; 729 coeffs[0] = z0(term, kappa, theta, sigma, v0, rho); 730 coeffs[1] = z1(term, kappa, theta, sigma, v0, rho); 731 coeffs[2] = z2(term, kappa, theta, sigma, v0, rho); 732 coeffs[3] = z3(term, kappa, theta, sigma, v0, rho); 733 } 734 impliedVolatility(const Real strike,const Real forward) const735 Real LPP3HestonExpansion::impliedVolatility(const Real strike, 736 const Real forward) const { 737 Real x = log(strike/forward); 738 Real vol = coeffs[0]+x*(coeffs[1]+x*(coeffs[2]+x*(coeffs[3]))); 739 return std::max(1e-8,vol); 740 } 741 } 742