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