QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
integralhestonvarianceoptionengine.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) 2008 Lorella Fatone
5 Copyright (C) 2008 Francesca Mariani
6 Copyright (C) 2008 Maria Cristina Recchioni
7 Copyright (C) 2008 Francesco Zirilli
8 Copyright (C) 2008 StatPro Italia srl
9
10 This file is part of QuantLib, a free-software/open-source library
11 for financial quantitative analysts and developers - http://quantlib.org/
12
13 QuantLib is free software: you can redistribute it and/or modify it
14 under the terms of the QuantLib license. You should have received a
15 copy of the license along with this program; if not, please email
16 <quantlib-dev@lists.sf.net>. The license is also available online at
17 <http://quantlib.org/license.shtml>.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 FOR A PARTICULAR PURPOSE. See the license for more details.
22*/
23
24#include <ql/errors.hpp>
26#include <ql/functional.hpp>
27#include <complex>
28#include <utility>
29#include <memory>
30
31namespace QuantLib {
32
33 namespace {
34
35 /*
36 *****************************************************************
37 **
38 ** Parameters defining the initial condition of the Heston model
39 ** and the European call option
40 **
41 *****************************************************************
42 */
43 /*
44 *****************************************************************
45 ** Assign: v0, eprice, tau, rtax
46 ******************************************************************
47 ******************************************************************
48 ** v0: initial variance
49 ** eprice: realized variance strike price
50 ** tau: time to maturity
51 * rtax: risk free interest rate
52 ****************************************************************
53 */
54
55 typedef std::complex<Real> Complex;
56
57 Real IvopOneDim(Real eps, Real chi, Real theta, Real /*rho*/,
58 Real v0, Real eprice, Time tau, Real rtax)
59 {
60 Real ss=0.0;
61 std::unique_ptr<double[]> xiv(new double[2048*2048+1]);
62 double nris=0.0;
63 int j=0,mm=0;
64 double pi=0,pi2=0;
65 double dstep=0;
66 Real option=0, impart=0;
67
68 std::unique_ptr<Complex[]> ff(new Complex[2048*2048]);
69 Complex xi;
70 Complex ui,beta,zita,gamma,csum,vero;
71 Complex contrib, caux, caux1,caux2,caux3;
72
73 ui=Complex(0.0,1.0);
74
75 /*
76 **********************************************************
77 ** i0: initial integrated variance i0=0
78 **********************************************************
79 */
80 Real i0=0.0;
81 //s=2.0*chi*theta/(eps*eps)-1.0;
82
83 //s=s+1;
84
85 /*
86 *************************************************
87 ** Start integration procedure
88 *************************************************
89 */
90
91 pi= 3.14159265358979324;
92 pi2=2.0*pi;
93 Real s=2.0*chi*theta/(eps*eps)-1.0;
94 /*
95 ****************************************
96 ** Note that s must be greater than zero
97 ****************************************
98 */
99
100 if(s<=0)
101 {
102 QL_FAIL("this parameter must be greater than zero-> " << s);
103 }
104
105 ss=s+1;
106
107 /*
108 *************************************************
109 ** Start integration procedure
110 *************************************************
111
112 **************************************************************
113 ** The oscillatory integral that approximates the price of
114 ** the realized variance option is computed using the method
115 ** proposed by Bailey, Swarztrauber in the paper published in
116 ** Siam Journal on Scientific Computing Vol 15(5) 1994
117 ** p. 1105-1110
118 **************************************************************
119
120 **************************************************************
121 ** dstep: real number, generally a power of two, that must be
122 ** assigned to determine the grid of
123 ** integration. Hint: dstep=256 or 512 (dstep<=2048)
124 **************************************************************
125 */
126 dstep=256.0;
127 nris=std::sqrt(pi2)/dstep;
128 mm=(int)(pi2/(nris*nris));
129
130 /*
131 ******************************************
132 ** Definition of the integration grid **
133 ******************************************
134 */
135 for (j=0;j<=mm-1;j++)
136 {
137 xiv[j+1]=(j-mm/2.0)*nris;
138 }
139
140 for (j=0;j<=mm-1;j++)
141 {
142 xi=xiv[j+1];
143 caux=chi*chi;
144 caux1=2.0*eps*eps;
145 caux1=caux1*xi;
146 caux1=caux1*ui;
147 caux2=caux1+caux;
148
149 zita=0.5*std::sqrt(caux2);
150
151 caux1=std::exp(-2.0*tau*zita);
152
153 beta=0.5*chi+zita;
154 beta=beta+caux1*(zita-0.5*chi);
155 gamma=1.0-caux1;
156
157 caux=-ss*tau;
158 caux2=caux*(zita-0.5*chi);
159 caux=ss*std::log(2.0*(zita/beta));
160 caux3=-v0*ui*xi*(gamma/beta);
161 caux=caux+caux3;
162 caux=caux+caux2;
163
164 ff[j+1]=std::exp(caux);
165 if(std::sqrt(std::imag(xi)*std::imag(xi)+std::real(xi)*std::real(xi))>1.e-06)
166 {
167 contrib=-eprice/(ui*xi);
168 caux=ui*xi;
169 caux=caux*eprice;
170 caux=std::exp(caux);
171 caux=caux-1.0;
172 caux2=ui*xi*ui*xi;
173 contrib=contrib+caux/caux2;
174 }
175 else
176 {
177 contrib=eprice*eprice*0.5;
178 }
179 ff[j+1]=ff[j+1]*contrib;
180 }
181 csum=0.0;
182 for (j=0;j<=mm-1;j++)
183 {
184 caux=std::pow(-1.0,j);
185 caux2=-2.0*pi*(double)mm*(double)j*0.5/(double)mm;
186 caux3=ui*caux2;
187 csum=csum+ff[j+1]*caux*std::exp(caux3);
188 }
189 csum=csum*std::sqrt(std::pow(-1.0,mm))*nris/pi2;
190 vero=i0-eprice+theta*tau+(1.0-std::exp(-chi*tau))*(v0-theta)/chi;
191 csum=csum+vero;
192 option=std::exp(-rtax*tau)*std::real(csum);
193 impart=std::imag(csum);
194 QL_ENSURE(impart <= 1e-12,
195 "imaginary part option (must be zero) = " << impart);
196 return option;
197 }
198
199
200
201 Real IvopTwoDim(Real eps, Real chi, Real theta, Real /*rho*/,
202 Real v0, Time tau, Real rtax,
203 const ext::function<Real(Real)>& payoff) {
204
205 Real ss=0.0;
206 std::unique_ptr<double[]> xiv(new double[2048*2048+1]);
207 std::unique_ptr<double[]> ivet(new double[2048 * 2048 + 1]);
208 double nris=0.0;
209 int j=0,mm=0,k=0;
210 double pi=0,pi2=0;
211
212 double dstep=0;
213 Real ip=0;
214 Real payoffval=0;
215 Real option=0/*, impart=0*/;
216
217 Real sumr=0;//,sumi=0;
218 Complex dxi,z;
219
220 std::unique_ptr<Complex[]> ff(new Complex[2048*2048]);
221 Complex xi;
222 Complex ui,beta,zita,gamma,csum;
223 Complex caux,caux1,caux2,caux3;
224
225 ui=Complex(0.0,1.0);
226
227 /*
228 **********************************************************
229 ** i0: initial integrated variance i0=0
230 **********************************************************
231 */
232 Real i0=0.0;
233
234 /*
235 *************************************************
236 ** Start integration procedure
237 *************************************************
238 */
239
240 pi= 3.14159265358979324;
241 pi2=2.0*pi;
242
243 Real s=2.0*chi*theta/(eps*eps)-1.0;
244 /*
245 ****************************************
246 ** Note that s must be greater than zero
247 ****************************************
248 */
249
250 if(s<=0)
251 {
252 QL_FAIL("this parameter must be greater than zero-> " << s);
253 }
254
255 ss=s+1;
256
257 /*
258 *************************************************
259 ** Start integration procedure
260 *************************************************
261
262 **************************************************************
263 ** The oscillatory integral that approximates the price of
264 ** the realized variance option is computed using the method
265 ** proposed by Bailey, Swarztrauber in the paper published in
266 ** Siam Journal on Scientific Computing Vol 15(5) 1994
267 ** p. 1105-1110
268 **************************************************************
269
270 **************************************************************
271 ** dstep: real number, generally a power of two that must be
272 ** assigned to determine the grid of
273 ** integration. Hint: dstep=256 or 512 (dstep<=2048)
274 **************************************************************
275 */
276 dstep=64.0;
277 nris=std::sqrt(pi2)/dstep;
278 mm=(int)(pi2/(nris*nris));
279
280 /*
281 ******************************************
282 ** Definition of the integration grid **
283 ******************************************
284 */
285
286 for (j=0;j<=mm-1;j++)
287 {
288 xiv[j+1]=(j-mm/2.0)*nris;
289 ivet[j+1]=(j-mm/2.0)*pi2/((double)mm*nris);
290 }
291
292 for (j=0;j<=mm-1;j++)
293 {
294 xi=xiv[j+1];
295
296 caux=chi*chi;
297 caux1=2.0*eps*eps;
298 caux1=caux1*xi;
299 caux1=caux1*ui;
300 caux2=caux1+caux;
301
302 zita=0.5*std::sqrt(caux2);
303 caux1=std::exp(-2.0*tau*zita);
304
305 beta=0.5*chi+zita;
306 beta=beta+caux1*(zita-0.5*chi);
307
308 gamma=1.0-caux1;
309
310 caux=-ss*tau;
311 caux2=caux*(zita-0.5*chi);
312 caux=ss*std::log(2.0*(zita/beta));
313 caux3=-v0*ui*xi*(gamma/beta);
314 caux=caux+caux3;
315 caux=caux+caux2;
316 ff[j+1]=std::exp(caux);
317 }
318
319 sumr=0.0;
320 //sumi=0.0;
321 for (k=0;k<=mm-1;k++)
322 {
323 ip=i0-ivet[k+1];
324 payoffval=payoff(ip);
325
326 dxi=2.0*pi*(double)k/(double)mm*ui;
327 csum=0.0;
328 for (j=0;j<=mm-1;j++)
329 {
330 z=-(double)j*dxi;
331 caux=std::pow(-1.0,j);
332 csum=csum+ff[j+1]*caux*std::exp(z);
333 }
334 csum=csum*std::pow(-1.0,k)*nris/pi2;
335
336 sumr=sumr+payoffval*std::real(csum);
337 //sumi=sumi+payoffval*std::imag(csum);
338 }
339 sumr=sumr*nris;
340 //sumi=sumi*nris;
341
342 option=std::exp(-rtax*tau)*sumr;
343 //impart=sumi;
344 //QL_ENSURE(impart <= 1e-3,
345 // "imaginary part option (must be close to zero) = " << impart);
346 return option;
347 }
348
349 struct payoff_adapter {
350 ext::shared_ptr<QuantLib::Payoff> payoff;
351 explicit payoff_adapter(ext::shared_ptr<QuantLib::Payoff> payoff)
352 : payoff(std::move(payoff)) {}
353 Real operator()(Real S) const {
354 return (*payoff)(S);
355 }
356 };
357
358 }
359
361 ext::shared_ptr<HestonProcess> process)
362 : process_(std::move(process)) {
364 }
365
367
368 QL_REQUIRE(process_->dividendYield().empty(),
369 "this engine does not manage dividend yields");
370
371 Handle<YieldTermStructure> riskFreeRate = process_->riskFreeRate();
372
373 Real epsilon = process_->sigma();
374 Real chi = process_->kappa();
375 Real theta = process_->theta();
376 Real rho = process_->rho();
377 Real v0 = process_->v0();
378
379 Time tau = riskFreeRate->dayCounter().yearFraction(
380 Settings::instance().evaluationDate(),
382 Rate r = riskFreeRate->zeroRate(arguments_.maturityDate,
383 riskFreeRate->dayCounter(),
384 Continuous);
385
386 ext::shared_ptr<PlainVanillaPayoff> plainPayoff =
387 ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
388 if ((plainPayoff != nullptr) && plainPayoff->optionType() == Option::Call) {
389 // a specialization for Call options is available
390 Real strike = plainPayoff->strike();
391 results_.value = IvopOneDim(epsilon, chi, theta, rho,
392 v0, strike, tau, r)
394 } else {
395 results_.value = IvopTwoDim(epsilon, chi, theta, rho, v0, tau, r,
396 payoff_adapter(arguments_.payoff))
398 }
399 }
400
401}
402
Shared handle to an observable.
Definition: handle.hpp:41
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
static Settings & instance()
access to the unique instance
Definition: singleton.hpp:104
ext::shared_ptr< Payoff > payoff
Classes and functions for error handling.
#define QL_ENSURE(condition, message)
throw an error if the given post-condition is not verified
Definition: errors.hpp:130
#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
Maps function, bind and cref to either the boost or std implementation.
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
Real Rate
interest rates
Definition: types.hpp:70
Real theta
Real v0
Real rho
ext::shared_ptr< QuantLib::Payoff > payoff
integral Heston-model variance-option engine
Definition: any.hpp:35
STL namespace.
ext::shared_ptr< YieldTermStructure > r
Real beta
Definition: sabr.cpp:200