QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
coshestonengine.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2017 Klaus Spanderen
5 Copyright (C) 2022 Ignacio Anguita
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
23
24namespace QuantLib {
25
27 const ext::shared_ptr<HestonModel>& model, Real L, Size N)
30 VanillaOption::results>(model),
31 L_(L), N_(N),
32 kappa_(model_->kappa()),
33 theta_(model_->theta()),
34 sigma_(model_->sigma()),
35 rho_ (model_->rho()) ,
36 v0_ (model_->v0()) { }
37
38
40 kappa_ = model_->kappa();
41 theta_ = model_->theta();
42 sigma_ = model_->sigma();
43 rho_ = model_->rho();
44 v0_ = model_->v0();
45
49 }
50
51
53
54 // this is a european option pricer
55 QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
56 "not an European option");
57
58 // plain vanilla
59 const ext::shared_ptr<PlainVanillaPayoff> payoff =
60 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
61 QL_REQUIRE(payoff, "non plain vanilla payoff given");
62
63 const ext::shared_ptr<HestonProcess> process = model_->process();
64
65 const Date maturityDate = arguments_.exercise->lastDate();
66 const Time maturity = process->time(maturityDate);
67
68 const Real cum1 = c1(maturity);
69 const Real w = std::sqrt(std::fabs(c2(maturity))
70 // the 4th order doesn't necessarily improve the precision
71 // + std::sqrt(std::fabs(c4(maturity)))
72 );
73
74 const Real k = payoff->strike();
75 const Real spot = process->s0()->value();
76 QL_REQUIRE(spot > 0.0, "negative or null underlying given");
77
78 const DiscountFactor df
79 = process->riskFreeRate()->discount(maturityDate);
80 const DiscountFactor qf
81 = process->dividendYield()->discount(maturityDate);
82 const Real fwd = spot*qf/df;
83 const Real x = std::log(fwd/k);
84
85 const Real a = x + cum1 - L_*w;
86 const Real b = x + cum1 + L_*w;
87
88 // Check if it exceeds the truncation bound
89
90 if (x >= b/2 || x <= a/2) {
91 //returns lower/upper bounds
92 if (payoff->optionType() == Option::Put)
93 results_.value = std::max(-spot*qf+k*df,0.0);
94 else if (payoff->optionType() == Option::Call)
95 results_.value = std::max(spot*qf-k*df,0.0);
96 else
97 QL_FAIL("unknown payoff type");
98 return;
99 }
100
101
102 const Real d = 1.0/(b-a);
103
104 const Real expA = std::exp(a);
105 Real s = chF(0, maturity).real()*(expA-1-a)*d;
106
107 for (Size n=1; n < N_; ++n) {
108 const Real r = n*M_PI*d;
109 const Real U_n = 2.0*d*( 1.0/(1.0 + r*r)
110 *(expA + r*std::sin(r*a) - std::cos(r*a)) - 1.0/r*std::sin(r*a));
111
112 s += U_n*(chF(r, maturity)
113 *std::exp(std::complex<Real>(0, r*(x-a)))).real();
114 }
115
116 if (payoff->optionType() == Option::Put)
117 results_.value = k*df*s;
118 else if (payoff->optionType() == Option::Call) {
119 const DiscountFactor qf
120 = process->dividendYield()->discount(maturityDate);
121 results_.value = spot*qf - k*df*(1-s);
122 }
123 else
124 QL_FAIL("unknown payoff type");
125 }
126
128 return std::log( model_->process()->dividendYield()->discount(t)
129 / model_->process()->riskFreeRate()->discount(t));
130 }
131
132 std::complex<Real> COSHestonEngine::chF(Real u, Real t)
133 const {
134 const Real sigma2 = sigma_*sigma_;
135
136 const std::complex<Real> D = std::sqrt(
137 squared(std::complex<Real>(kappa_, -rho_*sigma_*u))
138 + std::complex<Real>(u*u, u)*sigma2);
139
140 const std::complex<Real> g(kappa_, -rho_*sigma_*u);
141
142 const std::complex<Real> G = (g-D)/(g+D);
143
144 return std::exp(
145 v0_/(sigma2)*(1.0-std::exp(-D*t))/(1.0-G*std::exp(-D*t))
146 *(g-D) + kappa_*theta_/sigma2*((g-D)*t
147 -2.0*std::log((1.0-G*std::exp(-D*t))/(1.0-G)))
148 );
149 }
150
151 /*
152 Mathematica program to calculate the n-th cumulant
153
154 d[z_] := Sqrt[(kappa -i*rho*sigma*z)^2 + (z*z+i*z)*sigma^2]
155
156 g[z_] := (kappa -i*rho*sigma*z - d[z])/(kappa -i*rho*sigma*z + d[z])
157
158 phi[z_] := Exp[ v0/(sigma^2)*(1-Exp[-d[z]*t])/(1-g[z]*Exp[-d[z]*t])
159 *(kappa -i*rho*sigma*z - d[z]) + kappa*theta/sigma^2
160 *((kappa -i*rho*sigma*z-d[z])*t
161 -2*Log[(1-g[z]*Exp[-d[z]*t])/(1-g[z]) ]) ]
162
163 e[z_] := Log[phi[-i*z]]
164
165 // for C++ formatting see
166 // http://mathematica.stackexchange.com/questions/114175/cform-not-getting-exp-pow-log-functions
167 cpp = RawBoxes[Replace[ToBoxes@#,
168 InterpretationBox[a_, b_, c___] :> With[{aa =
169 StringReplace[ a, {"Sqrt" -> "std::sqrt", "Power(E," -> "std::exp(",
170 "Power" -> "std::pow"}]}, aa], {0, Infinity}]] &;
171
172 c[n_] := cpp@CForm[FullSimplify[Derivative[n][e][0],
173 kappa > 0 && theta > 0 && v0 > 0 && sigma > 0 &&
174 rho [Element] {-1, 1} && i^2 == -1]]
175 */
176
178 return (-theta_ + std::exp(kappa_*t)
179 *( theta_ - kappa_*t*theta_ -
180 v0_) + v0_)/(2*std::exp(kappa_*t)*kappa_);
181 }
182
184 const Real sigma2 = sigma_*sigma_;
185 const Real kappa2 = kappa_*kappa_;
186 const Real kappa3 = kappa2*kappa_;
187
188 return (sigma2*(theta_ - 2*v0_) +
189 std::exp(2*kappa_*t)*(8*kappa3*t*theta_ -
190 8*kappa2*(theta_ + rho_*sigma_*t*theta_ - v0_) +
191 sigma2*(-5*theta_ + 2*v0_) + 2*kappa_*sigma_*(8*rho_*theta_ +
192 sigma_*t*theta_ - 4*rho_*v0_)) +
193 4*std::exp(kappa_*t)*(sigma2*theta_ -
194 2*kappa2*(-1 + rho_*sigma_*t)*(theta_ - v0_) +
195 kappa_*sigma_*(sigma_*t*(theta_ - v0_) + 2*rho_*(-2*theta_ +
196 v0_))))/(8.*std::exp(2*kappa_*t)*kappa3);
197 }
198
200 const Real sigma2 = sigma_*sigma_;
201 const Real sigma3 = sigma2*sigma_;
202 const Real kappa2 = kappa_*kappa_;
203 const Real kappa3 = kappa2*kappa_;
204 const Real kappa4 = kappa3*kappa_;
205 const Real rho2 = rho_*rho_;
206
207 return
208 -(sigma_*(sigma3*(theta_ - 3*v0_) +
209 std::exp(3*kappa_*t)*(2*(-11*sigma3 -
210 24*kappa4*rho_*t + 3*kappa_*sigma2*(20*rho_ +
211 sigma_*t) - 6*kappa2*sigma_*(5 + 3*rho_*(4*rho_ + sigma_*t)) +
212 12*kappa3*(sigma_*t + 2*rho_*(2 + rho_*sigma_*t)))*theta_ -
213 6*(2*kappa_*rho_ - sigma_)*(4*kappa2 - 4*kappa_*rho_*sigma_ +
214 sigma2)*v0_) + 6*std::exp(kappa_*t)*sigma_*(-2*kappa2*(-1 +
215 rho_*sigma_*t)*(theta_ - 2*v0_) + sigma2*(theta_ - v0_) +
217 3*std::exp(2*kappa_*t)*(2*kappa_*sigma2*(-16*rho_*theta_ +
218 sigma_*t*(3*theta_ - v0_)) + 8*kappa4*rho_*t*(-2 +
219 rho_*sigma_*t)*(theta_ - v0_) + sigma3*(5*theta_ + v0_) +
220 8*kappa3*(-(rho_*(4 + sigma2*t*t)*theta_) + 2*sigma_*t*(theta_ - v0_) +
221 2*rho2*sigma_*t*(2*theta_ - v0_) + rho_*(2 +
222 sigma2*t*t)*v0_) + 2*kappa2*sigma_*((8
223 + 24*rho2 - 16*rho_*sigma_*t + sigma2*t*t)*theta_ - (8*rho2 -
224 8*rho_*sigma_*t + sigma2*t*t)*v0_))))/(16.*std::exp(3*kappa_*t)*
225 kappa_*kappa4);
226 }
227
229 const Real sigma2 = sigma_*sigma_;
230 const Real sigma3 = sigma2*sigma_;
231 const Real sigma4 = sigma2*sigma2;
232 const Real kappa2 = kappa_*kappa_;
233 const Real kappa3 = kappa2*kappa_;
234 const Real kappa4 = kappa2*kappa2;
235 const Real kappa5 = kappa2*kappa3;
236 const Real kappa6 = kappa3*kappa3;
237 const Real kappa7 = kappa4*kappa3;
238 const Real rho2 = rho_*rho_;
239 const Real rho3 = rho2*rho_;
240 const Time t2 = t*t;
241 const Time t3 = t2*t;
242
243 return
244 (sigma2*(3*sigma4*(theta_ - 4*v0_) +
245 3*std::exp(4*kappa_*t)*((-93*sigma4 +
246 64*kappa5*(t + 4*rho2*t) +
247 4*kappa_*sigma3*(176*rho_ + 5*sigma_*t) -
248 32*kappa2*sigma2*(11 + 50*rho2 +
249 5*rho_*sigma_*t) + 32*kappa3*sigma_*(3*sigma_*t + 4*rho_*(10 +
250 8*rho2 + 3*rho_*sigma_*t)) - 32*kappa4*(5 +
251 4*rho_*(6*rho_ + (3 + 2*rho2)*sigma_*t)))*theta_ +
252 4*(4*kappa2 - 4*kappa_*rho_*sigma_ +
253 sigma2)*(4*kappa2*(1 + 4*rho2) -
254 20*kappa_*rho_*sigma_ + 5*sigma2)*v0_) +
255 24*std::exp(kappa_*t)*sigma2*(-2*kappa2*(-1 +
256 rho_*sigma_*t)*(theta_ - 3*v0_) + sigma2*(theta_ - 2*v0_) +
258 3*sigma_*t*v0_)) + 12*std::exp(2*kappa_*t)*(sigma4*(7*theta_ -
259 4*v0_) + 8*kappa4*(1 + 2*rho_*sigma_*t*(-2 +
260 rho_*sigma_*t))*(theta_ - 2*v0_) +
261 2*kappa_*sigma3*(-24*rho_*theta_ + 5*sigma_*t*theta_ +
262 20*rho_*v0_ - 6*sigma_*t*v0_) + 4*kappa2*sigma2*((6
263 + 20*rho2 - 14*rho_*sigma_*t +
264 sigma2*t2)*theta_ - 2*(3 + 12*rho2 -
265 10*rho_*sigma_*t + sigma2*t2)*v0_) +
266 8*kappa3*sigma_*((3*sigma_*t + 2*rho_*(-4 + sigma_*t*(4*rho_ -
267 sigma_*t)))*theta_ + 2*(-3*sigma_*t + 2*rho_*(3 + sigma_*t*(-3*rho_ +
268 sigma_*t)))*v0_)) -
269 8*std::exp(3*kappa_*t)*(16*kappa6*rho2*t2*(-3 + rho_*sigma_*t)*(theta_ - v0_) - 3*sigma4*(7*theta_ +
270 2*v0_) + 2*kappa3*sigma_*((192*(rho_ + rho3) -
271 6*(9 + 40*rho2)*sigma_*t + 42*rho_*sigma2*t2 -
272 sigma3*t3)*theta_ + (-48*rho3 + 18*(1
273 + 4*rho2)*sigma_*t - 24*rho_*sigma2*t2
274 + sigma3*t3)*v0_) + 12*kappa4*((-4 -
275 24*rho2 + 8*rho_*(4 + 3*rho2)*sigma_*t - (3 +
276 14*rho2)*sigma2*t2 + rho_*sigma3*t3)*theta_ + (8*rho2 -
277 8*rho_*(2 + rho2)*sigma_*t + (3 +
278 8*rho2)*sigma2*t2 - rho_*sigma3*t3)*v0_) -
279 6*kappa2*sigma2*((15 + 80*rho2 -
280 35*rho_*sigma_*t + 2*sigma2*t2)*theta_ + (3 +
281 sigma_*t*(7*rho_ - sigma_*t))*v0_) + 24*kappa5*t*((-2 +
282 rho_*(4*sigma_*t + rho_*(-8 + sigma_*t*(4*rho_ - sigma_*t))))*theta_ + (2 +
283 rho_*(-4*sigma_*t + rho_*(4 + sigma_*t*(-2*rho_ + sigma_*t))))*v0_) +
284 3*kappa_*sigma3*(sigma_*t*(-9*theta_ + v0_) + 10*rho_*(6*theta_
285 + v0_)))))/(64.*std::exp(4*kappa_*t)*kappa7);
286 }
287
289 return c1(t);
290 }
292 return c2(t);
293 }
295 return c3(t)/std::pow(c2(t), 1.5);
296 }
298 return c4(t)/squared(c2(t));
299 }
300
301}
302
const Real kappa_
void calculate() const override
std::complex< Real > chF(Real u, Real t) const
Real kurtosis(Time t) const
COSHestonEngine(const ext::shared_ptr< HestonModel > &model, Real L=16, Size N=200)
Concrete date class.
Definition: date.hpp:125
Base class for some pricing engine on a particular model.
Heston model for the stochastic volatility of an asset.
Definition: hestonmodel.hpp:42
basic option arguments
Definition: option.hpp:57
Vanilla option (no discrete dividends, no barriers) on a single asset.
Heston engine based on Fourier-Cosine series expansions.
const DefaultType & t
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
ext::function< Real(Real)> b
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real DiscountFactor
discount factor between dates
Definition: types.hpp:66
std::size_t Size
size of a container
Definition: types.hpp:58
Real kappa
Real theta
Real v0
Real rho
Real sigma
ext::shared_ptr< QuantLib::Payoff > payoff
functionals and combinators not included in the STL
#define M_PI
Definition: any.hpp:35
T squared(T x)
Definition: functional.hpp:37
ext::shared_ptr< YieldTermStructure > r