1	/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3	/*
4	Copyright (C) 2005, 2007, 2009, 2010, 2012, 2014, 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
21	#include "hestonmodel.hpp"
22	#include "utilities.hpp"
23	#include <ql/experimental/exoticoptions/analyticpdfhestonengine.hpp>
24	#include <ql/instruments/dividendbarrieroption.hpp>
25	#include <ql/instruments/dividendvanillaoption.hpp>
26	#include <ql/math/integrals/gausslobattointegral.hpp>
27	#include <ql/math/optimization/differentialevolution.hpp>
28	#include <ql/math/optimization/levenbergmarquardt.hpp>
29	#include <ql/math/randomnumbers/rngtraits.hpp>
30	#include <ql/math/functional.hpp>
31	#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
32	#include <ql/methods/montecarlo/pathgenerator.hpp>
33	#include <ql/models/equity/hestonmodel.hpp>
34	#include <ql/models/equity/hestonmodelhelper.hpp>
35	#include <ql/models/equity/piecewisetimedependenthestonmodel.hpp>
36	#include <ql/pricingengines/barrier/fdblackscholesbarrierengine.hpp>
37	#include <ql/pricingengines/barrier/fdhestonbarrierengine.hpp>
38	#include <ql/pricingengines/blackformula.hpp>
39	#include <ql/pricingengines/vanilla/analyticdividendeuropeanengine.hpp>
40	#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
41	#include <ql/pricingengines/vanilla/analyticptdhestonengine.hpp>
42	#include <ql/pricingengines/vanilla/coshestonengine.hpp>
43	#include <ql/pricingengines/vanilla/exponentialfittinghestonengine.hpp>
44	#include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>
45	#include <ql/pricingengines/vanilla/fdhestonvanillaengine.hpp>
46	#include <ql/pricingengines/vanilla/hestonexpansionengine.hpp>
47	#include <ql/pricingengines/vanilla/mceuropeanhestonengine.hpp>
48	#include <ql/processes/hestonprocess.hpp>
49	#include <ql/quotes/simplequote.hpp>
50	#include <ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp>
51	#include <ql/termstructures/yield/flatforward.hpp>
52	#include <ql/termstructures/yield/zerocurve.hpp>
53	#include <ql/time/calendars/nullcalendar.hpp>
54	#include <ql/time/calendars/target.hpp>
55	#include <ql/time/daycounters/actual360.hpp>
56	#include <ql/time/daycounters/actual365fixed.hpp>
57	#include <ql/time/daycounters/actualactual.hpp>
58	#include <ql/time/period.hpp>
59	#include <cmath>
60	#include <utility>
61
62	using namespace QuantLib;
63	using namespace boost::unit_test_framework;
64
65	namespace {
66
67	struct CalibrationMarketData {
68	Handle<Quote> s0;
69	Handle<YieldTermStructure> riskFreeTS, dividendYield;
70	std::vector<ext::shared_ptr<CalibrationHelper> > options;
71	};
72
73	CalibrationMarketData getDAXCalibrationMarketData() {
74	/* this example is taken from A. Sepp
75	Pricing European-Style Options under Jump Diffusion Processes
76	with Stochstic Volatility: Applications of Fourier Transform
77	http://math.ut.ee/~spartak/papers/stochjumpvols.pdf
78	*/
79
80	Date settlementDate(Settings::instance().evaluationDate());
81
82	DayCounter dayCounter = Actual365Fixed();
83	Calendar calendar = TARGET();
84
85	Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 };
86	Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 };
87
88	std::vector<Date> dates;
89	std::vector<Rate> rates;
90	dates.push_back(x: settlementDate);
91	rates.push_back(x: 0.0357);
92	Size i;
93	for (i = 0; i < 8; ++i) {
94	dates.push_back(x: settlementDate + t[i]);
95	rates.push_back(x: r[i]);
96	}
97	Handle<YieldTermStructure> riskFreeTS(
98	ext::make_shared<ZeroCurve>(args&: dates, args&: rates, args&: dayCounter));
99
100	Handle<YieldTermStructure> dividendYield(
101	flatRate(today: settlementDate, forward: 0.0, dc: dayCounter));
102
103	Volatility v[] =
104	{ 0.6625,0.4875,0.4204,0.3667,0.3431,0.3267,0.3121,0.3121,
105	0.6007,0.4543,0.3967,0.3511,0.3279,0.3154,0.2984,0.2921,
106	0.5084,0.4221,0.3718,0.3327,0.3155,0.3027,0.2919,0.2889,
107	0.4541,0.3869,0.3492,0.3149,0.2963,0.2926,0.2819,0.2800,
108	0.4060,0.3607,0.3330,0.2999,0.2887,0.2811,0.2751,0.2775,
109	0.3726,0.3396,0.3108,0.2781,0.2788,0.2722,0.2661,0.2686,
110	0.3550,0.3277,0.3012,0.2781,0.2781,0.2661,0.2661,0.2681,
111	0.3428,0.3209,0.2958,0.2740,0.2688,0.2627,0.2580,0.2620,
112	0.3302,0.3062,0.2799,0.2631,0.2573,0.2533,0.2504,0.2544,
113	0.3343,0.2959,0.2705,0.2540,0.2504,0.2464,0.2448,0.2462,
114	0.3460,0.2845,0.2624,0.2463,0.2425,0.2385,0.2373,0.2422,
115	0.3857,0.2860,0.2578,0.2399,0.2357,0.2327,0.2312,0.2351,
116	0.3976,0.2860,0.2607,0.2356,0.2297,0.2268,0.2241,0.2320 };
117
118	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 4468.17));
119	Real strike[] = { 3400,3600,3800,4000,4200,4400,
120	4500,4600,4800,5000,5200,5400,5600 };
121
122	std::vector<ext::shared_ptr<CalibrationHelper> > options;
123
124	for (Size s = 0; s < 13; ++s) {
125	for (Size m = 0; m < 8; ++m) {
126	Handle<Quote> vol(ext::make_shared<SimpleQuote>(args&: v[s*8+m]));
127
128	Period maturity((int)((t[m]+3)/7.), Weeks); // round to weeks
129	options.push_back(x: ext::make_shared<HestonModelHelper>(args&: maturity, args&: calendar,
130	args&: s0, args&: strike[s], args&: vol,
131	args&: riskFreeTS, args&: dividendYield,
132	args: BlackCalibrationHelper::ImpliedVolError));
133	}
134	}
135
136	CalibrationMarketData marketData = { .s0: s0, .riskFreeTS: riskFreeTS, .dividendYield: dividendYield, .options: options };
137
138	return marketData;
139	}
140
141	}
142
143
144	void HestonModelTest::testBlackCalibration() {
145	BOOST_TEST_MESSAGE(
146	"Testing Heston model calibration using a flat volatility surface...");
147
148	/* calibrate a Heston model to a constant volatility surface without
149	smile. expected result is a vanishing volatility of the volatility.
150	In addition theta and v0 should be equal to the constant variance */
151
152	Date today = Date::todaysDate();
153	Settings::instance().evaluationDate() = today;
154
155	DayCounter dayCounter = Actual360();
156	Calendar calendar = NullCalendar();
157
158	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.04, dc: dayCounter));
159	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.50, dc: dayCounter));
160
161	std::vector<Period> optionMaturities = {1 * Months, 2 * Months, 3 * Months, 6 * Months,
162	9 * Months, 1 * Years, 2 * Years};
163
164	std::vector<ext::shared_ptr<CalibrationHelper> > options;
165	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.0));
166	Handle<Quote> vol(ext::make_shared<SimpleQuote>(args: 0.1));
167	Volatility volatility = vol->value();
168
169	for (auto& optionMaturitie : optionMaturities) {
170	for (Real moneyness = -1.0; moneyness < 2.0; moneyness += 1.0) {
171	const Time tau = dayCounter.yearFraction(
172	d1: riskFreeTS->referenceDate(),
173	d2: calendar.advance(date: riskFreeTS->referenceDate(), period: optionMaturitie));
174	const Real fwdPrice =
175	s0->value() * dividendTS->discount(t: tau) / riskFreeTS->discount(t: tau);
176	const Real strikePrice = fwdPrice * std::exp(x: -moneyness * volatility * std::sqrt(x: tau));
177
178	options.push_back(x: ext::make_shared<HestonModelHelper>(
179	args&: optionMaturitie, args&: calendar, args&: s0, args: strikePrice, args&: vol, args&: riskFreeTS, args&: dividendTS));
180	}
181	}
182
183	for (Real sigma = 0.1; sigma < 0.7; sigma += 0.2) {
184	const Real v0=0.01;
185	const Real kappa=0.2;
186	const Real theta=0.02;
187	const Real rho=-0.75;
188
189	ext::shared_ptr<HestonProcess> process(
190	ext::make_shared<HestonProcess>(args&: riskFreeTS, args&: dividendTS,
191	args&: s0, args: v0, args: kappa, args: theta, args&: sigma, args: rho));
192
193	ext::shared_ptr<HestonModel> model(ext::make_shared<HestonModel>(args&: process));
194	ext::shared_ptr<PricingEngine> engine(
195	ext::make_shared<AnalyticHestonEngine>(args&: model, args: 96));
196
197	for (auto& option : options)
198	ext::dynamic_pointer_cast<BlackCalibrationHelper>(r: option)->setPricingEngine(engine);
199
200	LevenbergMarquardt om(1e-8, 1e-8, 1e-8);
201	model->calibrate(options, method&: om, endCriteria: EndCriteria(400, 40, 1.0e-8,
202	1.0e-8, 1.0e-8));
203
204	Real tolerance = 3.0e-3;
205
206	if (model->sigma() > tolerance) {
207	BOOST_ERROR("Failed to reproduce expected sigma"
208	<< "\n calculated: " << model->sigma()
209	<< "\n expected: " << 0.0
210	<< "\n tolerance: " << tolerance);
211	}
212
213	if (std::fabs(x: model->kappa()
214	*(model->theta()-volatility*volatility)) > tolerance) {
215	BOOST_ERROR("Failed to reproduce expected theta"
216	<< "\n calculated: " << model->theta()
217	<< "\n expected: " << volatility*volatility);
218	}
219
220	if (std::fabs(x: model->v0()-volatility*volatility) > tolerance) {
221	BOOST_ERROR("Failed to reproduce expected v0"
222	<< "\n calculated: " << model->v0()
223	<< "\n expected: " << volatility*volatility);
224	}
225	}
226	}
227
228
229	void HestonModelTest::testDAXCalibration() {
230
231	BOOST_TEST_MESSAGE(
232	"Testing Heston model calibration using DAX volatility data...");
233
234	Date settlementDate(5, July, 2002);
235	Settings::instance().evaluationDate() = settlementDate;
236
237	CalibrationMarketData marketData = getDAXCalibrationMarketData();
238
239	const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;
240	const Handle<YieldTermStructure> dividendTS = marketData.dividendYield;
241	const Handle<Quote> s0 = marketData.s0;
242
243	const std::vector<ext::shared_ptr<CalibrationHelper> >& options = marketData.options;
244
245	const Real v0=0.1;
246	const Real kappa=1.0;
247	const Real theta=0.1;
248	const Real sigma=0.5;
249	const Real rho=-0.5;
250
251	const ext::shared_ptr<HestonProcess> process(
252	ext::make_shared<HestonProcess>(
253	args: riskFreeTS, args: dividendTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
254
255	const ext::shared_ptr<HestonModel> model(
256	ext::make_shared<HestonModel>(args: process));
257
258	const ext::shared_ptr<PricingEngine> engines[] = {
259	ext::make_shared<AnalyticHestonEngine>(args: model, args: 64),
260	ext::make_shared<COSHestonEngine>(args: model, args: 12, args: 75),
261	ext::make_shared<ExponentialFittingHestonEngine>(args: model)
262	};
263
264	const Array params = model->params();
265	for (const auto& engine : engines) {
266	model->setParams(params);
267	for (const auto& option : options)
268	ext::dynamic_pointer_cast<BlackCalibrationHelper>(r: option)->setPricingEngine(engine);
269
270	LevenbergMarquardt om(1e-8, 1e-8, 1e-8);
271	model->calibrate(options, method&: om,
272	endCriteria: EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8));
273
274	Real sse = 0;
275	for (Size i = 0; i < 13*8; ++i) {
276	const Real diff = options[i]->calibrationError()*100.0;
277	sse += diff*diff;
278	}
279	Real expected = 177.2; //see article by A. Sepp.
280	if (std::fabs(x: sse - expected) > 1.0) {
281	BOOST_FAIL("Failed to reproduce calibration error"
282	<< "\n calculated: " << sse
283	<< "\n expected: " << expected);
284	}
285	}
286	}
287
288	void HestonModelTest::testAnalyticVsBlack() {
289	BOOST_TEST_MESSAGE("Testing analytic Heston engine against Black formula...");
290
291	Date settlementDate = Date::todaysDate();
292	Settings::instance().evaluationDate() = settlementDate;
293	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
294	Date exerciseDate = settlementDate + 6*Months;
295
296	ext::shared_ptr<StrikedTypePayoff> payoff(
297	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 30));
298	ext::shared_ptr<Exercise> exercise(
299	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
300
301	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.1, dc: dayCounter));
302	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.04, dc: dayCounter));
303
304	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 32.0));
305
306	const Real v0=0.05;
307	const Real kappa=5.0;
308	const Real theta=0.05;
309	const Real sigma=1.0e-4;
310	const Real rho=0.0;
311
312	ext::shared_ptr<HestonProcess> process(ext::make_shared<HestonProcess>(
313	args&: riskFreeTS, args&: dividendTS, args&: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
314
315	VanillaOption option(payoff, exercise);
316	ext::shared_ptr<PricingEngine> engine(
317	ext::make_shared<AnalyticHestonEngine>(
318	args: ext::make_shared<HestonModel>(args&: process), args: 144));
319
320	option.setPricingEngine(engine);
321	Real calculated = option.NPV();
322
323	Real yearFraction = dayCounter.yearFraction(d1: settlementDate, d2: exerciseDate);
324	Real forwardPrice = 32*std::exp(x: (0.1-0.04)*yearFraction);
325	Real expected = blackFormula(optionType: payoff->optionType(), strike: payoff->strike(),
326	forward: forwardPrice, stdDev: std::sqrt(x: 0.05*yearFraction)) *
327	std::exp(x: -0.1*yearFraction);
328	Real error = std::fabs(x: calculated - expected);
329	Real tolerance = 2.0e-7;
330	if (error > tolerance) {
331	BOOST_FAIL("failed to reproduce Black price with AnalyticHestonEngine"
332	<< "\n calculated: " << calculated
333	<< "\n expected: " << expected
334	<< "\n error: " << std::scientific << error);
335	}
336
337	engine =
338	ext::make_shared<FdHestonVanillaEngine>(
339	args: ext::make_shared<HestonModel>(args&: process),
340	args: 200,args: 200,args: 100);
341	option.setPricingEngine(engine);
342
343	calculated = option.NPV();
344	error = std::fabs(x: calculated - expected);
345	tolerance = 1.0e-3;
346	if (error > tolerance) {
347	BOOST_FAIL("failed to reproduce Black price with FdHestonVanillaEngine"
348	<< "\n calculated: " << calculated
349	<< "\n expected: " << expected
350	<< "\n error: " << std::scientific << error);
351	}
352
353	}
354
355
356	void HestonModelTest::testAnalyticVsCached() {
357	BOOST_TEST_MESSAGE("Testing analytic Heston engine against cached values...");
358
359	Date settlementDate(27, December, 2004);
360	Settings::instance().evaluationDate() = settlementDate;
361	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
362	Date exerciseDate(28, March, 2005);
363
364	ext::shared_ptr<StrikedTypePayoff> payoff(
365	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 1.05));
366	ext::shared_ptr<Exercise> exercise(
367	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
368
369	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.0225, dc: dayCounter));
370	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.02, dc: dayCounter));
371
372	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.0));
373	const Real v0 = 0.1;
374	const Real kappa = 3.16;
375	const Real theta = 0.09;
376	const Real sigma = 0.4;
377	const Real rho = -0.2;
378
379	ext::shared_ptr<HestonProcess> process(ext::make_shared<HestonProcess>(
380	args&: riskFreeTS, args&: dividendTS, args&: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
381
382	VanillaOption option(payoff, exercise);
383
384	ext::shared_ptr<AnalyticHestonEngine> engine(
385	ext::make_shared<AnalyticHestonEngine>(
386	args: ext::make_shared<HestonModel>(args&: process), args: 64));
387
388	option.setPricingEngine(engine);
389
390	Real expected1 = 0.0404774515;
391	Real calculated1 = option.NPV();
392	Real tolerance = 1.0e-8;
393
394	if (std::fabs(x: calculated1 - expected1) > tolerance) {
395	BOOST_ERROR("Failed to reproduce cached analytic price"
396	<< "\n calculated: " << calculated1
397	<< "\n expected: " << expected1);
398	}
399
400
401	// reference values from www.wilmott.com, technical forum
402	// search for "Heston or VG price check"
403
404	Real K[] = {0.9,1.0,1.1};
405	Real expected2[] = { 0.1330371,0.0641016, 0.0270645 };
406	Real calculated2[6];
407
408	Size i;
409	for (i = 0; i < 6; ++i) {
410	Date exerciseDate(8+i/3, September, 2005);
411
412	ext::shared_ptr<StrikedTypePayoff> payoff(
413	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args&: K[i%3]));
414	ext::shared_ptr<Exercise> exercise(
415	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
416
417	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.05, dc: dayCounter));
418	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.02, dc: dayCounter));
419
420	Real s = riskFreeTS->discount(t: 0.7)/dividendTS->discount(t: 0.7);
421	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args&: s));
422
423	ext::shared_ptr<HestonProcess> process(
424	ext::make_shared<HestonProcess>(
425	args&: riskFreeTS, args&: dividendTS, args&: s0, args: 0.09, args: 1.2, args: 0.08, args: 1.8, args: -0.45));
426
427	VanillaOption option(payoff, exercise);
428
429	ext::shared_ptr<PricingEngine> engine(
430	ext::make_shared<AnalyticHestonEngine>(
431	args: ext::make_shared<HestonModel>(args&: process)));
432
433	option.setPricingEngine(engine);
434	calculated2[i] = option.NPV();
435	}
436
437	// we are after the value for T=0.7
438	Time t1 = dayCounter.yearFraction(d1: settlementDate, d2: Date(8, September,2005));
439	Time t2 = dayCounter.yearFraction(d1: settlementDate, d2: Date(9, September,2005));
440
441	for (i = 0; i < 3; ++i) {
442	const Real interpolated =
443	calculated2[i]+(calculated2[i+3]-calculated2[i])/(t2-t1)*(0.7-t1);
444
445	if (std::fabs(x: interpolated - expected2[i]) > 100*tolerance) {
446	BOOST_ERROR("Failed to reproduce cached analytic prices:"
447	<< "\n calculated: " << interpolated
448	<< "\n expected: " << expected2[i] );
449	}
450	}
451	}
452
453
454	void HestonModelTest::testMcVsCached() {
455	BOOST_TEST_MESSAGE(
456	"Testing Monte Carlo Heston engine against cached values...");
457
458	Date settlementDate(27, December, 2004);
459	Settings::instance().evaluationDate() = settlementDate;
460
461	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
462	Date exerciseDate(28, March, 2005);
463
464	ext::shared_ptr<StrikedTypePayoff> payoff(
465	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 1.05));
466	ext::shared_ptr<Exercise> exercise(
467	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
468
469	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.7, dc: dayCounter));
470	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.4, dc: dayCounter));
471
472	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.05));
473
474	ext::shared_ptr<HestonProcess> process(
475	ext::make_shared<HestonProcess>(
476	args&: riskFreeTS, args&: dividendTS, args&: s0, args: 0.3, args: 1.16, args: 0.2, args: 0.8, args: 0.8,
477	args: HestonProcess::QuadraticExponentialMartingale));
478
479	VanillaOption option(payoff, exercise);
480
481	ext::shared_ptr<PricingEngine> engine;
482	engine = MakeMCEuropeanHestonEngine<PseudoRandom>(process)
483	.withStepsPerYear(steps: 11)
484	.withAntitheticVariate()
485	.withSamples(samples: 50000)
486	.withSeed(seed: 1234);
487
488	option.setPricingEngine(engine);
489
490	Real expected = 0.0632851308977151;
491	Real calculated = option.NPV();
492	Real errorEstimate = option.errorEstimate();
493	Real tolerance = 7.5e-4;
494
495	if (std::fabs(x: calculated - expected) > 2.34*errorEstimate) {
496	BOOST_ERROR("Failed to reproduce cached price"
497	<< "\n calculated: " << calculated
498	<< "\n expected: " << expected
499	<< " +/- " << errorEstimate);
500	}
501
502	if (errorEstimate > tolerance) {
503	BOOST_ERROR("failed to reproduce error estimate"
504	<< "\n calculated: " << errorEstimate
505	<< "\n expected: " << tolerance);
506	}
507	}
508
509	void HestonModelTest::testFdBarrierVsCached() {
510	BOOST_TEST_MESSAGE("Testing FD barrier Heston engine against cached values...");
511
512	DayCounter dc = Actual360();
513	Date today = Date::todaysDate();
514
515	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
516	Handle<YieldTermStructure> rTS(flatRate(today, forward: 0.08, dc));
517	Handle<YieldTermStructure> qTS(flatRate(today, forward: 0.04, dc));
518
519	Date exDate = today + 180;
520	ext::shared_ptr<Exercise> exercise(
521	ext::make_shared<EuropeanExercise>(args&: exDate));
522
523	ext::shared_ptr<StrikedTypePayoff> payoff(
524	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 90.0));
525
526	ext::shared_ptr<HestonProcess> process(
527	ext::make_shared<HestonProcess>(
528	args&: rTS, args&: qTS, args&: s0, args: 0.25*0.25, args: 1.0, args: 0.25*0.25, args: 0.001, args: 0.0));
529
530	ext::shared_ptr<PricingEngine> engine;
531	engine = ext::make_shared<FdHestonBarrierEngine>(
532	args: ext::make_shared<HestonModel>(args&: process),
533	args: 200,args: 400,args: 100);
534
535	BarrierOption option(Barrier::DownOut, 95.0, 3.0, payoff, exercise);
536	option.setPricingEngine(engine);
537
538	Real calculated = option.NPV();
539	Real expected = 9.0246;
540	Real error = std::fabs(x: calculated-expected);
541	if (error > 1.0e-3) {
542	BOOST_FAIL("failed to reproduce cached price with FD Barrier engine"
543	<< "\n calculated: " << calculated
544	<< "\n expected: " << expected
545	<< "\n error: " << std::scientific << error);
546	}
547
548	option = BarrierOption(Barrier::DownIn, 95.0, 3.0, payoff, exercise);
549	option.setPricingEngine(engine);
550
551	calculated = option.NPV();
552	expected = 7.7627;
553	error = std::fabs(x: calculated-expected);
554	if (error > 1.0e-3) {
555	BOOST_FAIL("failed to reproduce cached price with FD Barrier engine"
556	<< "\n calculated: " << calculated
557	<< "\n expected: " << expected
558	<< "\n error: " << std::scientific << error);
559	}
560	}
561
562	void HestonModelTest::testFdVanillaVsCached() {
563	BOOST_TEST_MESSAGE("Testing FD vanilla Heston engine against cached values...");
564
565	Date settlementDate(27, December, 2004);
566	Settings::instance().evaluationDate() = settlementDate;
567
568	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
569	Date exerciseDate(28, March, 2005);
570
571	ext::shared_ptr<StrikedTypePayoff> payoff(
572	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 1.05));
573	ext::shared_ptr<Exercise> exercise(
574	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
575
576	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.7, dc: dayCounter));
577	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.4, dc: dayCounter));
578
579	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.05));
580
581	VanillaOption option(payoff, exercise);
582
583	ext::shared_ptr<HestonProcess> process(
584	ext::make_shared<HestonProcess>(
585	args&: riskFreeTS, args&: dividendTS, args&: s0, args: 0.3, args: 1.16, args: 0.2, args: 0.8, args: 0.8));
586
587	option.setPricingEngine(
588	MakeFdHestonVanillaEngine(ext::make_shared<HestonModel>(args&: process))
589	.withTGrid(tGrid: 100)
590	.withXGrid(xGrid: 200)
591	.withVGrid(vGrid: 100)
592	);
593
594	Real expected = 0.06325;
595	Real calculated = option.NPV();
596	Real error = std::fabs(x: calculated - expected);
597	Real tolerance = 1.0e-4;
598
599	if (error > tolerance) {
600	BOOST_FAIL("failed to reproduce cached price with FD engine"
601	<< "\n calculated: " << calculated
602	<< "\n expected: " << expected
603	<< "\n error: " << std::scientific << error);
604	}
605	}
606
607	void HestonModelTest::testFdVanillaWithDividendsVsCached() {
608	BOOST_TEST_MESSAGE("Testing FD vanilla Heston engine for discrete dividends...");
609
610	Date settlementDate(27, December, 2004);
611	Settings::instance().evaluationDate() = settlementDate;
612
613	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
614
615	auto payoff = ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 95.0);
616
617	auto s0 = Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100.0));
618	auto riskFreeTS = Handle<YieldTermStructure>(flatRate(forward: 0.05, dc: dayCounter));
619	auto dividendTS = Handle<YieldTermStructure>(flatRate(forward: 0.0, dc: dayCounter));
620
621	auto exerciseDate = Date(28, March, 2006);
622	auto exercise = ext::make_shared<EuropeanExercise>(args&: exerciseDate);
623
624	std::vector<Date> dividendDates;
625	std::vector<Real> dividends;
626	for (Date d = settlementDate + 3*Months;
627	d < exercise->lastDate();
628	d += 6*Months) {
629	dividendDates.push_back(x: d);
630	dividends.push_back(x: 1.0);
631	}
632
633	QL_DEPRECATED_DISABLE_WARNING
634	DividendVanillaOption divOption(payoff, exercise,
635	dividendDates, dividends);
636	QL_DEPRECATED_ENABLE_WARNING
637	auto process = ext::make_shared<HestonProcess>(
638	args&: riskFreeTS, args&: dividendTS, args&: s0, args: 0.04, args: 1.0, args: 0.04, args: 0.001, args: 0.0);
639	divOption.setPricingEngine(
640	MakeFdHestonVanillaEngine(ext::make_shared<HestonModel>(args&: process))
641	.withTGrid(tGrid: 200)
642	.withXGrid(xGrid: 400)
643	.withVGrid(vGrid: 100)
644	);
645
646	Real calculated = divOption.NPV();
647	// Value calculated with an independent FD framework, validated with
648	// an independent MC framework
649	Real expected = 12.946;
650	Real error = std::fabs(x: calculated - expected);
651	Real tolerance = 5.0e-3;
652
653	if (error > tolerance) {
654	BOOST_FAIL("failed to reproduce discrete dividend price with FD engine"
655	<< "\n calculated: " << calculated
656	<< "\n expected: " << expected
657	<< "\n error: " << std::scientific << error);
658	}
659
660
661	VanillaOption option(payoff, exercise);
662	option.setPricingEngine(
663	MakeFdHestonVanillaEngine(ext::make_shared<HestonModel>(args&: process))
664	.withCashDividends(dividendDates, dividendAmounts: dividends)
665	.withTGrid(tGrid: 200)
666	.withXGrid(xGrid: 400)
667	.withVGrid(vGrid: 100)
668	);
669
670	calculated = option.NPV();
671	error = std::fabs(x: calculated - expected);
672
673	if (error > tolerance) {
674	BOOST_FAIL("failed to reproduce discrete dividend price with FD engine"
675	<< "\n calculated: " << calculated
676	<< "\n expected: " << expected
677	<< "\n error: " << std::scientific << error);
678	}
679	}
680
681	void HestonModelTest::testFdAmerican() {
682	BOOST_TEST_MESSAGE("Testing FD vanilla Heston engine for american exercise...");
683
684	Date settlementDate(27, December, 2004);
685	Settings::instance().evaluationDate() = settlementDate;
686
687	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
688
689	auto s0 = Handle<Quote>(ext::make_shared<SimpleQuote>(args: 100.0));
690	auto riskFreeTS = Handle<YieldTermStructure>(flatRate(forward: 0.05, dc: dayCounter));
691	auto dividendTS = Handle<YieldTermStructure>(flatRate(forward: 0.03, dc: dayCounter));
692	auto process = ext::make_shared<HestonProcess>(
693	args&: riskFreeTS, args&: dividendTS, args&: s0, args: 0.04, args: 1.0, args: 0.04, args: 0.001, args: 0.0);
694	auto payoff = ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: 95.0);
695	auto exerciseDate = Date(28, March, 2006);
696	auto exercise = ext::make_shared<AmericanExercise>(
697	args&: settlementDate, args&: exerciseDate);
698	auto option = VanillaOption(payoff, exercise);
699	option.setPricingEngine(
700	MakeFdHestonVanillaEngine(ext::make_shared<HestonModel>(args&: process))
701	.withTGrid(tGrid: 200)
702	.withXGrid(xGrid: 400)
703	.withVGrid(vGrid: 100)
704	);
705	Real calculated = option.NPV();
706
707	Handle<BlackVolTermStructure> volTS(flatVol(today: settlementDate, volatility: 0.2,
708	dc: dayCounter));
709	ext::shared_ptr<BlackScholesMertonProcess> ref_process(
710	ext::make_shared<BlackScholesMertonProcess>(args&: s0, args&: dividendTS, args&: riskFreeTS, args&: volTS));
711	ext::shared_ptr<PricingEngine> ref_engine(
712	ext::make_shared<FdBlackScholesVanillaEngine>(args&: ref_process, args: 200, args: 400));
713	option.setPricingEngine(ref_engine);
714	Real expected = option.NPV();
715
716	Real error = std::fabs(x: calculated - expected);
717	Real tolerance = 1.0e-3;
718
719	if (error > tolerance) {
720	BOOST_FAIL("failed to reproduce american option price with FD engine"
721	<< "\n calculated: " << calculated
722	<< "\n expected: " << expected
723	<< "\n error: " << std::scientific << error);
724	}
725	}
726
727	namespace {
728	struct HestonProcessDiscretizationDesc {
729	HestonProcess::Discretization discretization;
730	Size nSteps;
731	std::string name;
732	};
733	}
734
735	void HestonModelTest::testKahlJaeckelCase() {
736	BOOST_TEST_MESSAGE(
737	"Testing MC and FD Heston engines for the Kahl-Jaeckel example...");
738
739	/* Example taken from Wilmott mag (Sept. 2005).
740	"Not-so-complex logarithms in the Heston model",
741	Example was also discussed within the Wilmott thread
742	"QuantLib code is very high quatlity"
743	*/
744
745	Date settlementDate(30, March, 2007);
746	Settings::instance().evaluationDate() = settlementDate;
747
748	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
749	Date exerciseDate(30, March, 2017);
750
751	const ext::shared_ptr<StrikedTypePayoff> payoff(
752	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 200));
753	const ext::shared_ptr<Exercise> exercise(
754	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
755
756	VanillaOption option(payoff, exercise);
757
758
759	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.0, dc: dayCounter));
760	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.0, dc: dayCounter));
761
762	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100));
763
764	const Real v0 = 0.16;
765	const Real theta = v0;
766	const Real kappa = 1.0;
767	const Real sigma = 2.0;
768	const Real rho =-0.8;
769
770
771	const HestonProcessDiscretizationDesc descriptions[] = {
772	{ .discretization: HestonProcess::NonCentralChiSquareVariance, .nSteps: 10,
773	.name: "NonCentralChiSquareVariance" },
774	{ .discretization: HestonProcess::QuadraticExponentialMartingale, .nSteps: 100,
775	.name: "QuadraticExponentialMartingale" },
776	};
777
778	const Real tolerance = 0.2;
779	const Real expected = 4.95212;
780
781	for (const auto& description : descriptions) {
782	const ext::shared_ptr<HestonProcess> process(ext::make_shared<HestonProcess>(
783	args&: riskFreeTS, args&: dividendTS, args&: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho, args: description.discretization));
784
785	const ext::shared_ptr<PricingEngine> engine =
786	MakeMCEuropeanHestonEngine<PseudoRandom>(process)
787	.withSteps(steps: description.nSteps)
788	.withAntitheticVariate()
789	.withAbsoluteTolerance(tolerance)
790	.withSeed(seed: 1234);
791	option.setPricingEngine(engine);
792
793	const Real calculated = option.NPV();
794	const Real errorEstimate = option.errorEstimate();
795
796	if (std::fabs(x: calculated - expected) > 2.34*errorEstimate) {
797	BOOST_ERROR("Failed to reproduce cached price with MC engine"
798	<< "\n discretization: " << description.name
799	<< "\n expected: " << expected
800	<< "\n calculated: " << calculated << " +/- " << errorEstimate);
801	}
802
803	if (errorEstimate > tolerance) {
804	BOOST_ERROR("failed to reproduce error estimate with MC engine"
805	<< "\n discretization: " << description.name << "\n calculated : "
806	<< errorEstimate << "\n expected : " << tolerance);
807	}
808	}
809
810	option.setPricingEngine(
811	MakeMCEuropeanHestonEngine<LowDiscrepancy>(
812	ext::make_shared<HestonProcess>(
813	args&: riskFreeTS, args&: dividendTS, args&: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho,
814	args: HestonProcess::BroadieKayaExactSchemeLaguerre))
815	.withSteps(steps: 1)
816	.withSamples(samples: 1023));
817
818	Real calculated = option.NPV();
819	if (std::fabs(x: calculated - expected) > 0.5*tolerance) {
820	BOOST_ERROR("Failed to reproduce cached price with MC engine"
821	<< "\n discretization: BroadieKayaExactSchemeLobatto"
822	<< "\n calculated: " << calculated
823	<< "\n expected: " << expected
824	<< "\n tolerance: " << tolerance);
825	}
826
827
828	const ext::shared_ptr<HestonModel> hestonModel(
829	ext::make_shared<HestonModel>(
830	args: ext::make_shared<HestonProcess>(
831	args&: riskFreeTS, args&: dividendTS, args&: s0, args: v0,
832	args: kappa, args: theta, args: sigma, args: rho)));
833
834	option.setPricingEngine(
835	ext::make_shared<FdHestonVanillaEngine>(args: hestonModel, args: 200, args: 401, args: 101));
836
837	calculated = option.NPV();
838	Real error = std::fabs(x: calculated - expected);
839	if (error > 5.0e-2) {
840	BOOST_FAIL("failed to reproduce cached price with FD engine"
841	<< "\n calculated: " << calculated
842	<< "\n expected: " << expected
843	<< "\n error: " << std::scientific << error);
844	}
845
846	option.setPricingEngine(
847	ext::make_shared<AnalyticHestonEngine>(args: hestonModel, args: 1e-6, args: 1000));
848
849	calculated = option.NPV();
850	error = std::fabs(x: calculated - expected);
851
852	if (error > 0.00002) {
853	BOOST_FAIL("failed to reproduce cached price with "
854	"GaussLobatto engine"
855	<< "\n calculated: " << calculated
856	<< "\n expected: " << expected
857	<< "\n error: " << std::scientific << error);
858	}
859
860	option.setPricingEngine(
861	ext::make_shared<COSHestonEngine>(args: hestonModel, args: 16, args: 400));
862	calculated = option.NPV();
863	error = std::fabs(x: calculated - expected);
864
865	if (error > 0.00002) {
866	BOOST_FAIL("failed to reproduce cached price with "
867	"Cosine engine"
868	<< "\n calculated: " << calculated
869	<< "\n expected: " << expected
870	<< "\n error: " << std::scientific << error);
871	}
872
873	option.setPricingEngine(
874	ext::make_shared<ExponentialFittingHestonEngine>(args: hestonModel));
875	calculated = option.NPV();
876	error = std::fabs(x: calculated - expected);
877
878	if (error > 0.00002) {
879	BOOST_FAIL("failed to reproduce cached price with "
880	"exponential fitting Heston engine"
881	<< "\n calculated: " << calculated
882	<< "\n expected: " << expected
883	<< "\n error: " << std::scientific << error);
884	}
885	}
886
887	namespace {
888	struct HestonParameter {
889	Real v0, kappa, theta, sigma, rho; };
890	}
891
892	void HestonModelTest::testDifferentIntegrals() {
893	BOOST_TEST_MESSAGE(
894	"Testing different numerical Heston integration algorithms...");
895
896	const Date settlementDate(27, December, 2004);
897	Settings::instance().evaluationDate() = settlementDate;
898
899	const DayCounter dayCounter = ActualActual(ActualActual::ISDA);
900
901	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.05, dc: dayCounter));
902	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.03, dc: dayCounter));
903
904	const Real strikes[] = { 0.5, 0.7, 1.0, 1.25, 1.5, 2.0 };
905	const Integer maturities[] = { 1, 2, 3, 12, 60, 120, 360};
906	const Option::Type types[] ={ Option::Put, Option::Call };
907
908	const HestonParameter equityfx = { .v0: 0.07, .kappa: 2.0, .theta: 0.04, .sigma: 0.55, .rho: -0.8 };
909	const HestonParameter highCorr = { .v0: 0.07, .kappa: 1.0, .theta: 0.04, .sigma: 0.55, .rho: 0.995 };
910	const HestonParameter lowVolOfVol = { .v0: 0.07, .kappa: 1.0, .theta: 0.04, .sigma: 0.025, .rho: -0.75 };
911	const HestonParameter highVolOfVol = { .v0: 0.07, .kappa: 1.0, .theta: 0.04, .sigma: 5.0, .rho: -0.75 };
912	const HestonParameter kappaEqSigRho = { .v0: 0.07, .kappa: 0.4, .theta: 0.04, .sigma: 0.5, .rho: 0.8 };
913
914	std::vector<HestonParameter> params = {
915	equityfx,
916	highCorr,
917	lowVolOfVol,
918	highVolOfVol,
919	kappaEqSigRho
920	};
921	const Real tol[] = { 1e-3, 1e-3, 0.2, 0.01, 1e-3 };
922
923	for (std::vector<HestonParameter>::const_iterator iter = params.begin();
924	iter != params.end(); ++iter) {
925
926	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.0));
927	ext::shared_ptr<HestonProcess> process(
928	ext::make_shared<HestonProcess>(
929	args&: riskFreeTS, args&: dividendTS,
930	args&: s0, args: iter->v0, args: iter->kappa,
931	args: iter->theta, args: iter->sigma, args: iter->rho));
932
933	ext::shared_ptr<HestonModel> model(
934	ext::make_shared<HestonModel>(args&: process));
935
936	ext::shared_ptr<AnalyticHestonEngine> lobattoEngine(
937	ext::make_shared<AnalyticHestonEngine>(args&: model, args: 1e-10,
938	args: 1000000));
939	ext::shared_ptr<AnalyticHestonEngine> laguerreEngine(
940	ext::make_shared<AnalyticHestonEngine>(args&: model, args: 128));
941	ext::shared_ptr<AnalyticHestonEngine> legendreEngine(
942	ext::make_shared<AnalyticHestonEngine>(
943	args&: model, args: AnalyticHestonEngine::Gatheral,
944	args: AnalyticHestonEngine::Integration::gaussLegendre(integrationOrder: 512)));
945	ext::shared_ptr<AnalyticHestonEngine> chebyshevEngine(
946	ext::make_shared<AnalyticHestonEngine>(
947	args&: model, args: AnalyticHestonEngine::Gatheral,
948	args: AnalyticHestonEngine::Integration::gaussChebyshev(integrationOrder: 512)));
949	ext::shared_ptr<AnalyticHestonEngine> chebyshev2ndEngine(
950	ext::make_shared<AnalyticHestonEngine>(
951	args&: model, args: AnalyticHestonEngine::Gatheral,
952	args: AnalyticHestonEngine::Integration::gaussChebyshev2nd(integrationOrder: 512)));
953
954	Real maxLegendreDiff = 0.0;
955	Real maxChebyshevDiff = 0.0;
956	Real maxChebyshev2ndDiff= 0.0;
957	Real maxLaguerreDiff = 0.0;
958
959	for (int maturitie : maturities) {
960	ext::shared_ptr<Exercise> exercise(
961	ext::make_shared<EuropeanExercise>(args: settlementDate + Period(maturitie, Months)));
962
963	for (Real strike : strikes) {
964	for (auto type : types) {
965
966	ext::shared_ptr<StrikedTypePayoff> payoff(
967	ext::make_shared<PlainVanillaPayoff>(args&: type, args&: strike));
968
969	VanillaOption option(payoff, exercise);
970
971	option.setPricingEngine(lobattoEngine);
972	const Real lobattoNPV = option.NPV();
973
974	option.setPricingEngine(laguerreEngine);
975	const Real laguerre = option.NPV();
976
977	option.setPricingEngine(legendreEngine);
978	const Real legendre = option.NPV();
979
980	option.setPricingEngine(chebyshevEngine);
981	const Real chebyshev = option.NPV();
982
983	option.setPricingEngine(chebyshev2ndEngine);
984	const Real chebyshev2nd = option.NPV();
985
986	maxLaguerreDiff
987	= std::max(a: maxLaguerreDiff,
988	b: std::fabs(x: lobattoNPV-laguerre));
989	maxLegendreDiff
990	= std::max(a: maxLegendreDiff,
991	b: std::fabs(x: lobattoNPV-legendre));
992	maxChebyshevDiff
993	= std::max(a: maxChebyshevDiff,
994	b: std::fabs(x: lobattoNPV-chebyshev));
995	maxChebyshev2ndDiff
996	= std::max(a: maxChebyshev2ndDiff,
997	b: std::fabs(x: lobattoNPV-chebyshev2nd));
998	}
999	}
1000	}
1001	const Real maxDiff = std::max(a: std::max(
1002	a: std::max(a: maxLaguerreDiff,b: maxLegendreDiff),
1003	b: maxChebyshevDiff), b: maxChebyshev2ndDiff);
1004
1005	const Real tr = tol[iter - params.begin()];
1006	if (maxDiff > tr) {
1007	BOOST_ERROR("Failed to reproduce Heston pricing values "
1008	"within given tolerance"
1009	<< "\n maxDifference: " << maxDiff
1010	<< "\n tolerance: " << tr);
1011	}
1012	}
1013	}
1014
1015	void HestonModelTest::testMultipleStrikesEngine() {
1016	BOOST_TEST_MESSAGE("Testing multiple-strikes FD Heston engine...");
1017
1018	Date settlementDate(27, December, 2004);
1019	Settings::instance().evaluationDate() = settlementDate;
1020
1021	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
1022	Date exerciseDate(28, March, 2006);
1023
1024	ext::shared_ptr<Exercise> exercise(
1025	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
1026
1027	Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.06, dc: dayCounter));
1028	Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.02, dc: dayCounter));
1029
1030	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.05));
1031
1032	ext::shared_ptr<HestonProcess> process(
1033	ext::make_shared<HestonProcess>(
1034	args&: riskFreeTS, args&: dividendTS, args&: s0, args: 0.16, args: 2.5, args: 0.09, args: 0.8, args: -0.8));
1035	ext::shared_ptr<HestonModel> model(
1036	ext::make_shared<HestonModel>(args&: process));
1037
1038	std::vector<Real> strikes = {1.0, 0.5, 0.75, 1.5, 2.0};
1039
1040	ext::shared_ptr<FdHestonVanillaEngine> singleStrikeEngine(
1041	ext::make_shared<FdHestonVanillaEngine>(args&: model, args: 20, args: 400, args: 50));
1042	ext::shared_ptr<FdHestonVanillaEngine> multiStrikeEngine(
1043	ext::make_shared<FdHestonVanillaEngine>(args&: model, args: 20, args: 400, args: 50));
1044	multiStrikeEngine->enableMultipleStrikesCaching(strikes);
1045
1046	Real relTol = 5e-3;
1047	for (Real& strike : strikes) {
1048	ext::shared_ptr<StrikedTypePayoff> payoff(
1049	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args&: strike));
1050
1051	VanillaOption aOption(payoff, exercise);
1052	aOption.setPricingEngine(multiStrikeEngine);
1053
1054	Real npvCalculated = aOption.NPV();
1055	Real deltaCalculated = aOption.delta();
1056	Real gammaCalculated = aOption.gamma();
1057	Real thetaCalculated = aOption.theta();
1058
1059	aOption.setPricingEngine(singleStrikeEngine);
1060	Real npvExpected = aOption.NPV();
1061	Real deltaExpected = aOption.delta();
1062	Real gammaExpected = aOption.gamma();
1063	Real thetaExpected = aOption.theta();
1064
1065	if (std::fabs(x: npvCalculated-npvExpected)/npvExpected > relTol) {
1066	BOOST_FAIL("failed to reproduce price with FD multi strike engine"
1067	<< "\n calculated: " << npvCalculated
1068	<< "\n expected: " << npvExpected
1069	<< "\n error: " << std::scientific << relTol);
1070	}
1071	if (std::fabs(x: deltaCalculated-deltaExpected)/deltaExpected > relTol) {
1072	BOOST_FAIL("failed to reproduce delta with FD multi strike engine"
1073	<< "\n calculated: " << deltaCalculated
1074	<< "\n expected: " << deltaExpected
1075	<< "\n error: " << std::scientific << relTol);
1076	}
1077	if (std::fabs(x: gammaCalculated-gammaExpected)/gammaExpected > relTol) {
1078	BOOST_FAIL("failed to reproduce gamma with FD multi strike engine"
1079	<< "\n calculated: " << gammaCalculated
1080	<< "\n expected: " << gammaExpected
1081	<< "\n error: " << std::scientific << relTol);
1082	}
1083	if (std::fabs(x: thetaCalculated-thetaExpected)/thetaExpected > relTol) {
1084	BOOST_FAIL("failed to reproduce theta with FD multi strike engine"
1085	<< "\n calculated: " << thetaCalculated
1086	<< "\n expected: " << thetaExpected
1087	<< "\n error: " << std::scientific << relTol);
1088	}
1089	}
1090	}
1091
1092
1093
1094	void HestonModelTest::testAnalyticPiecewiseTimeDependent() {
1095	BOOST_TEST_MESSAGE("Testing analytic piecewise time dependent Heston prices...");
1096
1097	Date settlementDate(27, December, 2004);
1098	Settings::instance().evaluationDate() = settlementDate;
1099	DayCounter dayCounter = ActualActual(ActualActual::ISDA);
1100	Date exerciseDate(28, March, 2005);
1101
1102	ext::shared_ptr<StrikedTypePayoff> payoff(
1103	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 1.0));
1104	ext::shared_ptr<Exercise> exercise(
1105	ext::make_shared<EuropeanExercise>(args&: exerciseDate));
1106
1107	std::vector<Date> dates = {settlementDate, {1, January, 2007}};
1108	std::vector<Rate> irates = {0.0, 0.2};
1109	Handle<YieldTermStructure> riskFreeTS(
1110	ext::make_shared<ZeroCurve>(args&: dates, args&: irates, args&: dayCounter));
1111
1112	std::vector<Rate> qrates = {0.0, 0.3};
1113	Handle<YieldTermStructure> dividendTS(
1114	ext::make_shared<ZeroCurve>(args&: dates, args&: qrates, args&: dayCounter));
1115
1116	const Real v0 = 0.1;
1117	Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.0));
1118
1119	ConstantParameter theta(0.09, PositiveConstraint());
1120	ConstantParameter kappa(3.16, PositiveConstraint());
1121	ConstantParameter sigma(4.40, PositiveConstraint());
1122	ConstantParameter rho (-0.8, BoundaryConstraint(-1.0, 1.0));
1123
1124	ext::shared_ptr<PiecewiseTimeDependentHestonModel> model =
1125	ext::make_shared<PiecewiseTimeDependentHestonModel>(
1126	args&: riskFreeTS, args&: dividendTS,
1127	args&: s0, args: v0, args&: theta, args&: kappa,
1128	args&: sigma, args&: rho, args: TimeGrid(20.0, 2));
1129
1130	VanillaOption option(payoff, exercise);
1131
1132	ext::shared_ptr<HestonProcess> hestonProcess(
1133	ext::make_shared<HestonProcess>(
1134	args&: riskFreeTS, args&: dividendTS, args&: s0, args: v0,
1135	args: kappa(0.0), args: theta(0.0), args: sigma(0.0), args: rho(0.0)));
1136	ext::shared_ptr<HestonModel> hestonModel =
1137	ext::make_shared<HestonModel>(args&: hestonProcess);
1138	option.setPricingEngine(
1139	ext::make_shared<AnalyticHestonEngine>(args&: hestonModel));
1140
1141	const Real expected = option.NPV();
1142
1143	option.setPricingEngine(ext::shared_ptr<PricingEngine>(
1144	new AnalyticPTDHestonEngine(model)));
1145
1146	const Real calculatedGatheral = option.NPV();
1147	if (std::fabs(x: calculatedGatheral-expected) > 1e-12) {
1148	BOOST_ERROR("failed to reproduce Heston prices with Gatheral ChF"
1149	<< "\n calculated: " << calculatedGatheral
1150	<< "\n expected: " << expected);
1151	}
1152
1153	option.setPricingEngine(ext::shared_ptr<PricingEngine>(
1154	new AnalyticPTDHestonEngine(
1155	model,
1156	AnalyticPTDHestonEngine::AndersenPiterbarg,
1157	AnalyticPTDHestonEngine::Integration::gaussLaguerre(integrationOrder: 164))));
1158	const Real calculatedAndersenPiterbarg = option.NPV();
1159
1160	if (std::fabs(x: calculatedAndersenPiterbarg-expected) > 1e-8) {
1161	BOOST_ERROR("failed to reproduce Heston prices Andersen-Piterbarg"
1162	<< "\n calculated: " << calculatedAndersenPiterbarg
1163	<< "\n expected: " << expected);
1164	}
1165	}
1166
1167	void HestonModelTest::testDAXCalibrationOfTimeDependentModel() {
1168	BOOST_TEST_MESSAGE(
1169	"Testing time-dependent Heston model calibration...");
1170
1171	Date settlementDate(5, July, 2002);
1172	Settings::instance().evaluationDate() = settlementDate;
1173
1174	CalibrationMarketData marketData = getDAXCalibrationMarketData();
1175
1176	const Handle<YieldTermStructure> riskFreeTS = marketData.riskFreeTS;
1177	const Handle<YieldTermStructure> dividendTS = marketData.dividendYield;
1178	const Handle<Quote> s0 = marketData.s0;
1179
1180	const std::vector<ext::shared_ptr<CalibrationHelper> >& options = marketData.options;
1181
1182	std::vector<Time> modelTimes = {0.25, 10.0};
1183	const TimeGrid modelGrid(modelTimes.begin(), modelTimes.end());
1184
1185	const Real v0=0.1;
1186	ConstantParameter sigma( 0.5, PositiveConstraint());
1187	ConstantParameter theta( 0.1, PositiveConstraint());
1188	ConstantParameter rho( -0.5, BoundaryConstraint(-1.0, 1.0));
1189
1190	std::vector<Time> pTimes(1, 0.25);
1191	PiecewiseConstantParameter kappa(pTimes, PositiveConstraint());
1192
1193	for (Size i=0; i < pTimes.size()+1; ++i) {
1194	kappa.setParam(i, x: 10.0);
1195	}
1196
1197	ext::shared_ptr<PiecewiseTimeDependentHestonModel> model =
1198	ext::make_shared<PiecewiseTimeDependentHestonModel>(
1199	args: riskFreeTS, args: dividendTS,
1200	args: s0, args: v0, args&: theta, args&: kappa,
1201	args&: sigma, args&: rho, args: modelGrid);
1202
1203	const ext::shared_ptr<PricingEngine> engines[] = {
1204	ext::make_shared<AnalyticPTDHestonEngine>(args&: model),
1205	ext::make_shared<AnalyticPTDHestonEngine>(
1206	args&: model,
1207	args: AnalyticPTDHestonEngine::AndersenPiterbarg,
1208	args: AnalyticPTDHestonEngine::Integration::gaussLaguerre(integrationOrder: 64)),
1209	ext::make_shared<AnalyticPTDHestonEngine>(
1210	args&: model,
1211	args: AnalyticPTDHestonEngine::AndersenPiterbarg,
1212	args: AnalyticPTDHestonEngine::Integration::discreteTrapezoid(evaluation: 72))
1213	};
1214
1215	for (const auto& engine : engines) {
1216	for (const auto& option : options)
1217	ext::dynamic_pointer_cast<BlackCalibrationHelper>(r: option)->setPricingEngine(engine);
1218
1219	LevenbergMarquardt om(1e-8, 1e-8, 1e-8);
1220	model->calibrate(options, method&: om,
1221	endCriteria: EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8));
1222
1223	Real sse = 0;
1224	for (Size i = 0; i < 13*8; ++i) {
1225	const Real diff = options[i]->calibrationError()*100.0;
1226	sse += diff*diff;
1227	}
1228
1229	Real expected = 74.4;
1230	if (std::fabs(x: sse - expected) > 1.0) {
1231	BOOST_ERROR("Failed to reproduce calibration error"
1232	<< "\n calculated: " << sse
1233	<< "\n expected: " << expected);
1234	}
1235	}
1236	}
1237
1238	void HestonModelTest::testAlanLewisReferencePrices() {
1239	BOOST_TEST_MESSAGE("Testing Alan Lewis reference prices...");
1240
1241	/*
1242	* testing Alan Lewis reference prices posted in
1243	* http://wilmott.com/messageview.cfm?catid=34&threadid=90957
1244	*/
1245
1246	const Date settlementDate(5, July, 2002);
1247	Settings::instance().evaluationDate() = settlementDate;
1248
1249	const Date maturityDate(5, July, 2003);
1250	const ext::shared_ptr<Exercise> exercise(
1251	ext::make_shared<EuropeanExercise>(args: maturityDate));
1252
1253	const DayCounter dayCounter = Actual365Fixed();
1254	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.01, dc: dayCounter));
1255	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.02, dc: dayCounter));
1256
1257	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
1258
1259	const Real v0 = 0.04;
1260	const Real rho = -0.5;
1261	const Real sigma = 1.0;
1262	const Real kappa = 4.0;
1263	const Real theta = 0.25;
1264
1265	const ext::shared_ptr<HestonProcess> process(
1266	ext::make_shared<HestonProcess>(
1267	args: riskFreeTS, args: dividendTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1268	const ext::shared_ptr<HestonModel> model(
1269	ext::make_shared<HestonModel>(args: process));
1270
1271	const ext::shared_ptr<PricingEngine> laguerreEngine(
1272	ext::make_shared<AnalyticHestonEngine>(args: model, args: 128U));
1273
1274	const ext::shared_ptr<PricingEngine> gaussLobattoEngine(
1275	ext::make_shared<AnalyticHestonEngine>(args: model, QL_EPSILON, args: 100000U));
1276
1277	const ext::shared_ptr<PricingEngine> cosEngine(
1278	ext::make_shared<COSHestonEngine>(args: model, args: 20, args: 400));
1279
1280	const ext::shared_ptr<PricingEngine> exponentialFittingEngine(
1281	ext::make_shared<ExponentialFittingHestonEngine>(args: model));
1282
1283	const ext::shared_ptr<PricingEngine> andersenPiterbargEngine(
1284	new AnalyticHestonEngine(
1285	model,
1286	AnalyticHestonEngine::AndersenPiterbarg,
1287	AnalyticHestonEngine::Integration::discreteTrapezoid(evaluation: 92),
1288	QL_EPSILON));
1289
1290	const Real strikes[] = { 80, 90, 100, 110, 120 };
1291	const Option::Type types[] = { Option::Put, Option::Call };
1292	const ext::shared_ptr<PricingEngine> engines[]
1293	= { laguerreEngine, gaussLobattoEngine,
1294	cosEngine, andersenPiterbargEngine, exponentialFittingEngine };
1295
1296	const Real expectedResults[][2] = {
1297	{ 7.958878113256768285213263077598987193482161301733,
1298	26.774758743998854221382195325726949201687074848341 },
1299	{ 12.017966707346304987709573290236471654992071308187,
1300	20.933349000596710388139445766564068085476194042256 },
1301	{ 17.055270961270109413522653999411000974895436309183,
1302	16.070154917028834278213466703938231827658768230714 },
1303	{ 23.017825898442800538908781834822560777763225722188,
1304	12.132211516709844867860534767549426052805766831181 },
1305	{ 29.811026202682471843340682293165857439167301370697,
1306	9.024913483457835636553375454092357136489051667150 }
1307	};
1308
1309	const Real tol = 1e-12; // 3e-15 works on linux/ia32,
1310	// but keep some buffer for other platforms
1311
1312	for (Size i=0; i < LENGTH(strikes); ++i) {
1313	const Real strike = strikes[i];
1314
1315	for (Size j=0; j < LENGTH(types); ++j) {
1316	const Option::Type type = types[j];
1317
1318	for (Size k=0; k < LENGTH(engines); ++k) {
1319	const ext::shared_ptr<PricingEngine> engine = engines[k];
1320
1321	const ext::shared_ptr<StrikedTypePayoff> payoff(
1322	ext::make_shared<PlainVanillaPayoff>(args: type, args: strike));
1323
1324	VanillaOption option(payoff, exercise);
1325	option.setPricingEngine(engine);
1326
1327	const Real expected = expectedResults[i][j];
1328	const Real calculated = option.NPV();
1329	const Real relError = std::fabs(x: calculated-expected)/expected;
1330
1331	if (relError > tol || std::isnan(x: calculated)) {
1332	BOOST_ERROR(
1333	"failed to reproduce Alan Lewis Reference prices "
1334	<< "\n strike : " << strike
1335	<< "\n option type: " << type
1336	<< "\n engine type: " << k
1337	<< "\n rel. error : " << relError);
1338	}
1339	}
1340	}
1341	}
1342	}
1343
1344	void HestonModelTest::testAnalyticPDFHestonEngine() {
1345	BOOST_TEST_MESSAGE("Testing analytic PDF Heston engine...");
1346
1347	const Date settlementDate(5, January, 2014);
1348	Settings::instance().evaluationDate() = settlementDate;
1349
1350	const DayCounter dayCounter = Actual365Fixed();
1351	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.07, dc: dayCounter));
1352	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.185, dc: dayCounter));
1353
1354	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
1355
1356	const Real v0 = 0.1;
1357	const Real rho = -0.5;
1358	const Real sigma = 1.0;
1359	const Real kappa = 4.0;
1360	const Real theta = 0.05;
1361
1362	const ext::shared_ptr<HestonModel> model(
1363	ext::make_shared<HestonModel>(
1364	args: ext::make_shared<HestonProcess>(args: riskFreeTS, args: dividendTS,
1365	args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
1366
1367	const Real tol = 1e-6;
1368	const ext::shared_ptr<AnalyticPDFHestonEngine> pdfEngine(
1369	ext::make_shared<AnalyticPDFHestonEngine>(args: model, args: tol));
1370
1371	const ext::shared_ptr<PricingEngine> analyticEngine(
1372	ext::make_shared<AnalyticHestonEngine>(args: model, args: 178));
1373
1374	const Date maturityDate(5, July, 2014);
1375	const Time maturity = dayCounter.yearFraction(d1: settlementDate, d2: maturityDate);
1376	const ext::shared_ptr<Exercise> exercise(
1377	ext::make_shared<EuropeanExercise>(args: maturityDate));
1378
1379	// 1. check a plain vanilla call option
1380	for (Real strike=40; strike < 190; strike+=20) {
1381	const ext::shared_ptr<StrikedTypePayoff> vanillaPayoff(
1382	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args&: strike));
1383
1384	VanillaOption planVanillaOption(vanillaPayoff, exercise);
1385
1386	planVanillaOption.setPricingEngine(pdfEngine);
1387	const Real calculated = planVanillaOption.NPV();
1388
1389	planVanillaOption.setPricingEngine(analyticEngine);
1390	const Real expected = planVanillaOption.NPV();
1391
1392	if (std::fabs(x: calculated-expected) > 3*tol) {
1393	BOOST_FAIL(
1394	"failed to reproduce plain vanilla european prices with"
1395	" the analytic probability density engine"
1396	<< "\n strike : " << strike
1397	<< "\n expected : " << expected
1398	<< "\n calculated : " << calculated
1399	<< "\n diff : " << std::fabs(calculated-expected)
1400	<< "\n tol ; " << tol);
1401	}
1402	}
1403
1404	// 2. digital call option (approx. with a call spread)
1405	for (Real strike=40; strike < 190; strike+=10) {
1406	VanillaOption digitalOption(
1407	ext::make_shared<CashOrNothingPayoff>(args: Option::Call, args&: strike, args: 1.0),
1408	exercise);
1409	digitalOption.setPricingEngine(pdfEngine);
1410	const Real calculated = digitalOption.NPV();
1411
1412	const Real eps = 0.01;
1413	VanillaOption longCall(
1414	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike-eps),
1415	exercise);
1416	longCall.setPricingEngine(analyticEngine);
1417
1418	VanillaOption shortCall(
1419	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike+eps),
1420	exercise);
1421	shortCall.setPricingEngine(analyticEngine);
1422
1423	const Real expected = (longCall.NPV() - shortCall.NPV())/(2*eps);
1424	if (std::fabs(x: calculated-expected) > tol) {
1425	BOOST_FAIL(
1426	"failed to reproduce european digital prices with"
1427	" the analytic probability density engine"
1428	<< "\n strike : " << strike
1429	<< "\n expected : " << expected
1430	<< "\n calculated : " << calculated
1431	<< "\n diff : " << std::fabs(calculated-expected)
1432	<< "\n tol : " << tol);
1433	}
1434
1435	const DiscountFactor d = riskFreeTS->discount(d: maturityDate);
1436	const Real expectedCDF = 1.0 - expected/d;
1437	const Real calculatedCDF = pdfEngine->cdf(X: strike, t: maturity);
1438
1439	if (std::fabs(x: expectedCDF - calculatedCDF) > tol) {
1440	BOOST_FAIL(
1441	"failed to reproduce cumulative distribution function"
1442	<< "\n strike : " << strike
1443	<< "\n expected CDF : " << expectedCDF
1444	<< "\n calculated CDF: " << calculatedCDF
1445	<< "\n diff : "
1446	<< std::fabs(calculatedCDF-expectedCDF)
1447	<< "\n tol : " << tol);
1448
1449	}
1450	}
1451	}
1452
1453	void HestonModelTest::testExpansionOnAlanLewisReference() {
1454	BOOST_TEST_MESSAGE("Testing expansion on Alan Lewis reference prices...");
1455
1456	const Date settlementDate(5, July, 2002);
1457	Settings::instance().evaluationDate() = settlementDate;
1458
1459	const Date maturityDate(5, July, 2003);
1460	const ext::shared_ptr<Exercise> exercise =
1461	ext::make_shared<EuropeanExercise>(args: maturityDate);
1462
1463	const DayCounter dayCounter = Actual365Fixed();
1464	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.01, dc: dayCounter));
1465	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.02, dc: dayCounter));
1466
1467	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
1468
1469	const Real v0 = 0.04;
1470	const Real rho = -0.5;
1471	const Real sigma = 1.0;
1472	const Real kappa = 4.0;
1473	const Real theta = 0.25;
1474
1475	const ext::shared_ptr<HestonProcess> process =
1476	ext::make_shared<HestonProcess>(args: riskFreeTS, args: dividendTS, args: s0, args: v0,
1477	args: kappa, args: theta, args: sigma, args: rho);
1478	const ext::shared_ptr<HestonModel> model =
1479	ext::make_shared<HestonModel>(args: process);
1480
1481	const ext::shared_ptr<PricingEngine> lpp2Engine =
1482	ext::make_shared<HestonExpansionEngine>(args: model,
1483	args: HestonExpansionEngine::LPP2);
1484	//don't test Forde as it does not behave well on this example
1485	const ext::shared_ptr<PricingEngine> lpp3Engine =
1486	ext::make_shared<HestonExpansionEngine>(args: model,
1487	args: HestonExpansionEngine::LPP3);
1488
1489	const Real strikes[] = { 80, 90, 100, 110, 120 };
1490	const Option::Type types[] = { Option::Put, Option::Call };
1491	const ext::shared_ptr<PricingEngine> engines[]
1492	= { lpp2Engine, lpp3Engine };
1493
1494	const Real expectedResults[][2] = {
1495	{ 7.958878113256768285213263077598987193482161301733,
1496	26.774758743998854221382195325726949201687074848341 },
1497	{ 12.017966707346304987709573290236471654992071308187,
1498	20.933349000596710388139445766564068085476194042256 },
1499	{ 17.055270961270109413522653999411000974895436309183,
1500	16.070154917028834278213466703938231827658768230714 },
1501	{ 23.017825898442800538908781834822560777763225722188,
1502	12.132211516709844867860534767549426052805766831181 },
1503	{ 29.811026202682471843340682293165857439167301370697,
1504	9.024913483457835636553375454092357136489051667150 }
1505	};
1506
1507	const Real tol[2] = {1.003e-2, 3.645e-3};
1508
1509	for (Size i=0; i < LENGTH(strikes); ++i) {
1510	const Real strike = strikes[i];
1511
1512	for (Size j=0; j < LENGTH(types); ++j) {
1513	const Option::Type type = types[j];
1514
1515	for (Size k=0; k < LENGTH(engines); ++k) {
1516	const ext::shared_ptr<PricingEngine> engine = engines[k];
1517
1518	const ext::shared_ptr<StrikedTypePayoff> payoff =
1519	ext::make_shared<PlainVanillaPayoff>(args: type, args: strike);
1520
1521	VanillaOption option(payoff, exercise);
1522	option.setPricingEngine(engine);
1523
1524	const Real expected = expectedResults[i][j];
1525	const Real calculated = option.NPV();
1526	const Real relError = std::fabs(x: calculated-expected)/expected;
1527
1528	if (relError > tol[k]) {
1529	BOOST_ERROR(
1530	"failed to reproduce Alan Lewis Reference prices "
1531	<< "\n strike : " << strike
1532	<< "\n option type: " << type
1533	<< "\n engine type: " << k
1534	<< "\n rel. error : " << relError);
1535	}
1536	}
1537	}
1538	}
1539	}
1540
1541	void HestonModelTest::testExpansionOnFordeReference() {
1542	BOOST_TEST_MESSAGE("Testing expansion on Forde reference prices...");
1543
1544	const Real forward = 100.0;
1545	const Real v0 = 0.04;
1546	const Real rho = -0.4;
1547	const Real sigma = 0.2;
1548	const Real kappa = 1.15;
1549	const Real theta = 0.04;
1550
1551	const Real terms[] = {0.1, 1.0, 5.0, 10.0};
1552
1553	const Real strikes[] = { 60, 80, 90, 100, 110, 120, 140 };
1554
1555	const Real referenceVols[][7] = {
1556	{0.27284673574924445, 0.22360758200372477, 0.21023988547031242, 0.1990674789471587, 0.19118230678920461, 0.18721342919371017, 0.1899869903378507},
1557	{0.25200775151345, 0.2127275920953156, 0.20286528150874591, 0.19479398358151515, 0.18872591728967686, 0.18470857955411824, 0.18204457060905446},
1558	{0.21637821506229973, 0.20077227130455172, 0.19721753043236154, 0.1942233023784151, 0.191693211401571, 0.18955229722896752, 0.18491727548069495},
1559	{0.20672925973965342, 0.198583062164427, 0.19668274423922746, 0.1950420231354201, 0.193610364344706, 0.1923502827886502, 0.18934360917857015}
1560	};
1561
1562	const Real tol[][4] = {
1563	{0.06, 0.03, 0.03, 0.02},
1564	{0.15, 0.08, 0.04, 0.02},
1565	{0.06, 0.08, 1.0, 1.0} //forde breaks down for long maturities
1566	};
1567	const Real tolAtm[][4] = {
1568	{4e-6, 7e-4, 2e-3, 9e-4},
1569	{7e-6, 4e-4, 9e-4, 4e-4},
1570	{4e-4, 3e-2, 0.28, 1.0}
1571	};
1572	for (Size j=0; j < LENGTH(terms); ++j) {
1573	const Real term = terms[j];
1574	const ext::shared_ptr<HestonExpansion> lpp2 =
1575	ext::make_shared<LPP2HestonExpansion>(args: kappa, args: theta, args: sigma,
1576	args: v0, args: rho, args: term);
1577	const ext::shared_ptr<HestonExpansion> lpp3 =
1578	ext::make_shared<LPP3HestonExpansion>(args: kappa, args: theta, args: sigma,
1579	args: v0, args: rho, args: term);
1580	const ext::shared_ptr<HestonExpansion> forde =
1581	ext::make_shared<FordeHestonExpansion>(args: kappa, args: theta, args: sigma,
1582	args: v0, args: rho, args: term);
1583	const ext::shared_ptr<HestonExpansion> expansions[] = { lpp2, lpp3, forde };
1584	for (Size i=0; i < LENGTH(strikes); ++i) {
1585	const Real strike = strikes[i];
1586	for (Size k=0; k < LENGTH(expansions); ++k) {
1587	const ext::shared_ptr<HestonExpansion> expansion = expansions[k];
1588
1589	const Real expected = referenceVols[j][i];
1590	const Real calculated = expansion->impliedVolatility(strike, forward);
1591	const Real relError = std::fabs(x: calculated-expected)/expected;
1592	const Real refTol = strike == forward ? tolAtm[k][j] : tol[k][j];
1593	if (relError > refTol) {
1594	BOOST_ERROR(
1595	"failed to reproduce Forde reference vols "
1596	<< "\n strike : " << strike
1597	<< "\n expansion type: " << k
1598	<< "\n rel. error : " << relError);
1599	}
1600	}
1601	}
1602	}
1603	}
1604
1605	namespace {
1606	void reportOnIntegrationMethodTest(VanillaOption& option,
1607	const ext::shared_ptr<HestonModel>& model,
1608	const AnalyticHestonEngine::Integration& integration,
1609	AnalyticHestonEngine::ComplexLogFormula formula,
1610	bool isAdaptive,
1611	Real expected,
1612	Real tol,
1613	Size valuations,
1614	const std::string& method) {
1615
1616	if (integration.isAdaptiveIntegration() != isAdaptive)
1617	BOOST_ERROR(method << " is not an adaptive integration routine");
1618
1619	const ext::shared_ptr<AnalyticHestonEngine> engine =
1620	ext::make_shared<AnalyticHestonEngine>(
1621	args: model, args&: formula, args: integration, args: 1e-9);
1622
1623	option.setPricingEngine(engine);
1624	const Real calculated = option.NPV();
1625
1626	const Real error = std::fabs(x: calculated - expected);
1627
1628	if (std::isnan(x: error) || error > tol) {
1629	BOOST_ERROR("failed to reproduce simple Heston Pricing with "
1630	<< "\n integration method: " << method
1631	<< std::setprecision(12)
1632	<< "\n expected : " << expected
1633	<< "\n calculated : " << calculated
1634	<< "\n error : " << error);
1635	}
1636
1637	if ( valuations != Null<Size>()
1638	&& valuations != engine->numberOfEvaluations()) {
1639	BOOST_ERROR("nubmer of function evaluations does not match "
1640	<< "\n integration method : " << method
1641	<< "\n expected function calls : " << valuations
1642	<< "\n number of function calls: "
1643	<< engine->numberOfEvaluations());
1644	}
1645	}
1646	}
1647
1648	void HestonModelTest::testAllIntegrationMethods() {
1649	BOOST_TEST_MESSAGE("Testing semi-analytic Heston pricing with all "
1650	"integration methods...");
1651
1652	const Date settlementDate(7, February, 2017);
1653	Settings::instance().evaluationDate() = settlementDate;
1654
1655	const DayCounter dayCounter = Actual365Fixed();
1656	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.05, dc: dayCounter));
1657	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.075, dc: dayCounter));
1658
1659	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
1660
1661	const Real v0 = 0.1;
1662	const Real rho = -0.75;
1663	const Real sigma = 0.4;
1664	const Real kappa = 4.0;
1665	const Real theta = 0.05;
1666
1667	const ext::shared_ptr<HestonModel> model =
1668	ext::make_shared<HestonModel>(
1669	args: ext::make_shared<HestonProcess>(
1670	args: riskFreeTS, args: dividendTS,
1671	args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1672
1673	const ext::shared_ptr<StrikedTypePayoff> payoff =
1674	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: s0->value());
1675
1676	const Date maturityDate = settlementDate + Period(1, Years);
1677	const ext::shared_ptr<Exercise> exercise =
1678	ext::make_shared<EuropeanExercise>(args: maturityDate);
1679
1680	VanillaOption option(payoff, exercise);
1681
1682	const Real tol = 1e-8;
1683	const Real expected = 10.147041515497;
1684
1685	// Gauss-Laguerre with Gatheral logarithm integration method
1686	reportOnIntegrationMethodTest(option, model,
1687	integration: AnalyticHestonEngine::Integration::gaussLaguerre(),
1688	formula: AnalyticHestonEngine::Gatheral,
1689	isAdaptive: false, expected, tol, valuations: 256, method: "Gauss-Laguerre with Gatheral logarithm");
1690
1691	// Gauss-Laguerre with branch correction integration method
1692	reportOnIntegrationMethodTest(option, model,
1693	integration: AnalyticHestonEngine::Integration::gaussLaguerre(),
1694	formula: AnalyticHestonEngine::BranchCorrection,
1695	isAdaptive: false, expected, tol, valuations: 256, method: "Gauss-Laguerre with branch correction");
1696
1697	// Gauss-Laguerre with Andersen-Piterbarg integration method
1698	reportOnIntegrationMethodTest(option, model,
1699	integration: AnalyticHestonEngine::Integration::gaussLaguerre(),
1700	formula: AnalyticHestonEngine::AndersenPiterbarg,
1701	isAdaptive: false, expected, tol, valuations: 128,
1702	method: "Gauss-Laguerre with Andersen Piterbarg control variate");
1703
1704	// Gauss-Legendre with Gatheral logarithm integration method
1705	reportOnIntegrationMethodTest(option, model,
1706	integration: AnalyticHestonEngine::Integration::gaussLegendre(),
1707	formula: AnalyticHestonEngine::Gatheral,
1708	isAdaptive: false, expected, tol, valuations: 256, method: "Gauss-Legendre with Gatheral logarithm");
1709
1710	// Gauss-Legendre with branch correction integration method
1711	reportOnIntegrationMethodTest(option, model,
1712	integration: AnalyticHestonEngine::Integration::gaussLegendre(),
1713	formula: AnalyticHestonEngine::BranchCorrection,
1714	isAdaptive: false, expected, tol, valuations: 256, method: "Gauss-Legendre with branch correction");
1715
1716	// Gauss-Legendre with Andersen-Piterbarg integration method
1717	reportOnIntegrationMethodTest(option, model,
1718	integration: AnalyticHestonEngine::Integration::gaussLegendre(integrationOrder: 256),
1719	formula: AnalyticHestonEngine::AndersenPiterbarg,
1720	isAdaptive: false, expected, tol: 1e-4, valuations: 256,
1721	method: "Gauss-Legendre with Andersen Piterbarg control variate");
1722
1723	// Gauss-Chebyshev with Gatheral logarithm integration method
1724	reportOnIntegrationMethodTest(option, model,
1725	integration: AnalyticHestonEngine::Integration::gaussChebyshev(integrationOrder: 512),
1726	formula: AnalyticHestonEngine::Gatheral,
1727	isAdaptive: false, expected, tol: 1e-4, valuations: 1024, method: "Gauss-Chebyshev with Gatheral logarithm");
1728
1729	// Gauss-Chebyshev with branch correction integration method
1730	reportOnIntegrationMethodTest(option, model,
1731	integration: AnalyticHestonEngine::Integration::gaussChebyshev(integrationOrder: 512),
1732	formula: AnalyticHestonEngine::BranchCorrection,
1733	isAdaptive: false, expected, tol: 1e-4, valuations: 1024, method: "Gauss-Chebyshev with branch correction");
1734
1735	// Gauss-Chebyshev with Andersen-Piterbarg integration method
1736	reportOnIntegrationMethodTest(option, model,
1737	integration: AnalyticHestonEngine::Integration::gaussChebyshev(integrationOrder: 512),
1738	formula: AnalyticHestonEngine::AndersenPiterbarg,
1739	isAdaptive: false, expected, tol: 1e-4, valuations: 512,
1740	method: "Gauss-Laguerre with Andersen Piterbarg control variate");
1741
1742	// Gauss-Chebyshev2nd with Gatheral logarithm integration method
1743	reportOnIntegrationMethodTest(option, model,
1744	integration: AnalyticHestonEngine::Integration::gaussChebyshev2nd(integrationOrder: 512),
1745	formula: AnalyticHestonEngine::Gatheral,
1746	isAdaptive: false, expected, tol: 2e-4, valuations: 1024,
1747	method: "Gauss-Chebyshev2nd with Gatheral logarithm");
1748
1749	// Gauss-Chebyshev with branch correction integration method
1750	reportOnIntegrationMethodTest(option, model,
1751	integration: AnalyticHestonEngine::Integration::gaussChebyshev2nd(integrationOrder: 512),
1752	formula: AnalyticHestonEngine::BranchCorrection,
1753	isAdaptive: false, expected, tol: 2e-4, valuations: 1024,
1754	method: "Gauss-Chebyshev2nd with branch correction");
1755
1756	// Gauss-Chebyshev with Andersen-Piterbarg integration method
1757	reportOnIntegrationMethodTest(option, model,
1758	integration: AnalyticHestonEngine::Integration::gaussChebyshev2nd(integrationOrder: 512),
1759	formula: AnalyticHestonEngine::AndersenPiterbarg,
1760	isAdaptive: false, expected, tol: 2e-4, valuations: 512,
1761	method: "Gauss-Chebyshev2nd with Andersen Piterbarg control variate");
1762
1763	// Discrete Simpson rule with Gatheral logarithm
1764	reportOnIntegrationMethodTest(option, model,
1765	integration: AnalyticHestonEngine::Integration::discreteSimpson(evaluation: 512),
1766	formula: AnalyticHestonEngine::Gatheral,
1767	isAdaptive: false, expected, tol, valuations: 1024,
1768	method: "Discrete Simpson rule with Gatheral logarithm");
1769
1770	// Discrete Simpson rule with Andersen-Piterbarg integration method
1771	reportOnIntegrationMethodTest(option, model,
1772	integration: AnalyticHestonEngine::Integration::discreteSimpson(evaluation: 64),
1773	formula: AnalyticHestonEngine::AndersenPiterbarg,
1774	isAdaptive: false, expected, tol, valuations: 64,
1775	method: "Discrete Simpson rule with Andersen Piterbarg control variate");
1776
1777	// Discrete Trapezoid rule with Gatheral logarithm
1778	reportOnIntegrationMethodTest(option, model,
1779	integration: AnalyticHestonEngine::Integration::discreteTrapezoid(evaluation: 512),
1780	formula: AnalyticHestonEngine::Gatheral,
1781	isAdaptive: false, expected, tol: 2e-4, valuations: 1024,
1782	method: "Discrete Trapezoid rule with Gatheral logarithm");
1783
1784	// Discrete Trapezoid rule with Andersen-Piterbarg integration method
1785	reportOnIntegrationMethodTest(option, model,
1786	integration: AnalyticHestonEngine::Integration::discreteTrapezoid(evaluation: 64),
1787	formula: AnalyticHestonEngine::AndersenPiterbarg,
1788	isAdaptive: false, expected, tol, valuations: 64,
1789	method: "Discrete Trapezoid rule with Andersen Piterbarg control variate");
1790
1791	// Gauss-Lobatto with Gatheral logarithm
1792	reportOnIntegrationMethodTest(option, model,
1793	integration: AnalyticHestonEngine::Integration::gaussLobatto(relTolerance: tol, absTolerance: Null<Real>()),
1794	formula: AnalyticHestonEngine::Gatheral,
1795	isAdaptive: true, expected, tol, valuations: Null<Size>(),
1796	method: "Gauss-Lobatto with Gatheral logarithm");
1797
1798	// Gauss-Lobatto with Andersen-Piterbarg integration method
1799	reportOnIntegrationMethodTest(option, model,
1800	integration: AnalyticHestonEngine::Integration::gaussLobatto(relTolerance: tol, absTolerance: Null<Real>()),
1801	formula: AnalyticHestonEngine::AndersenPiterbarg,
1802	isAdaptive: true, expected, tol, valuations: Null<Size>(),
1803	method: "Gauss-Lobatto with Andersen Piterbarg control variate");
1804
1805	// Gauss-Konrod with Gatheral logarithm
1806	reportOnIntegrationMethodTest(option, model,
1807	integration: AnalyticHestonEngine::Integration::gaussKronrod(absTolerance: tol),
1808	formula: AnalyticHestonEngine::Gatheral,
1809	isAdaptive: true, expected, tol, valuations: Null<Size>(),
1810	method: "Gauss-Konrod with Gatheral logarithm");
1811
1812	// Gauss-Konrod with Andersen-Piterbarg integration method
1813	reportOnIntegrationMethodTest(option, model,
1814	integration: AnalyticHestonEngine::Integration::gaussKronrod(absTolerance: tol),
1815	formula: AnalyticHestonEngine::AndersenPiterbarg,
1816	isAdaptive: true, expected, tol, valuations: Null<Size>(),
1817	method: "Gauss-Konrod with Andersen Piterbarg control variate");
1818
1819	// Simpson with Gatheral logarithm
1820	reportOnIntegrationMethodTest(option, model,
1821	integration: AnalyticHestonEngine::Integration::simpson(absTolerance: tol),
1822	formula: AnalyticHestonEngine::Gatheral,
1823	isAdaptive: true, expected, tol: 1e-6, valuations: Null<Size>(),
1824	method: "Simpson with Gatheral logarithm");
1825
1826	// Simpson with Andersen-Piterbarg integration method
1827	reportOnIntegrationMethodTest(option, model,
1828	integration: AnalyticHestonEngine::Integration::simpson(absTolerance: tol),
1829	formula: AnalyticHestonEngine::AndersenPiterbarg,
1830	isAdaptive: true, expected, tol: 1e-6, valuations: Null<Size>(),
1831	method: "Simpson with Andersen Piterbarg control variate");
1832
1833	// Trapezoid with Gatheral logarithm
1834	reportOnIntegrationMethodTest(option, model,
1835	integration: AnalyticHestonEngine::Integration::trapezoid(absTolerance: tol),
1836	formula: AnalyticHestonEngine::Gatheral,
1837	isAdaptive: true, expected, tol: 1e-6, valuations: Null<Size>(),
1838	method: "Trapezoid with Gatheral logarithm");
1839
1840	// Trapezoid with Andersen-Piterbarg integration method
1841	reportOnIntegrationMethodTest(option, model,
1842	integration: AnalyticHestonEngine::Integration::trapezoid(absTolerance: tol),
1843	formula: AnalyticHestonEngine::AndersenPiterbarg,
1844	isAdaptive: true, expected, tol: 1e-6, valuations: Null<Size>(),
1845	method: "Trapezoid with Andersen Piterbarg control variate");
1846	}
1847
1848	namespace {
1849	class LogCharacteristicFunction {
1850	public:
1851	LogCharacteristicFunction(Size n, Time t, ext::shared_ptr<COSHestonEngine> engine)
1852	: t_(t), alpha_(0.0, 1.0), engine_(std::move(engine)) {
1853	for (Size i=1; i < n; ++i, alpha_*=std::complex<Real>(0,1));
1854	}
1855
1856	Real operator()(Real u) const {
1857	return (std::log(z: engine_->chF(u, t: t_))/alpha_).real();
1858	}
1859
1860	private:
1861	const Time t_;
1862	std::complex<Real> alpha_;
1863	const ext::shared_ptr<COSHestonEngine> engine_;
1864	};
1865	}
1866
1867	void HestonModelTest::testCosHestonCumulants() {
1868	BOOST_TEST_MESSAGE("Testing Heston COS cumulants...");
1869
1870	const Date settlementDate(7, February, 2017);
1871	Settings::instance().evaluationDate() = settlementDate;
1872
1873	const DayCounter dayCounter = Actual365Fixed();
1874	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.15, dc: dayCounter));
1875	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.075, dc: dayCounter));
1876
1877	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
1878
1879	const Real v0 = 0.1;
1880	const Real rho = -0.75;
1881	const Real sigma = 0.4;
1882	const Real kappa = 4.0;
1883	const Real theta = 0.25;
1884
1885	const ext::shared_ptr<HestonModel> model =
1886	ext::make_shared<HestonModel>(
1887	args: ext::make_shared<HestonProcess>(
1888	args: riskFreeTS, args: dividendTS,
1889	args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1890
1891	const ext::shared_ptr<COSHestonEngine> cosEngine =
1892	ext::make_shared<COSHestonEngine>(args: model);
1893
1894	const Real tol = 1e-7;
1895	const NumericalDifferentiation::Scheme central(
1896	NumericalDifferentiation::Central);
1897
1898	for (Time t=0.01; t < 41.0; t+=t) {
1899	const Real nc1 = NumericalDifferentiation(
1900	ext::function<Real(Real)>(
1901	LogCharacteristicFunction(1, t, cosEngine)),
1902	1, 1e-5, 5, central)(0.0);
1903
1904	const Real c1 = cosEngine->c1(t);
1905
1906	if (std::fabs(x: nc1 - c1) > tol) {
1907	BOOST_ERROR(" failed to reproduce first cumulant"
1908	<< "\n expected: " << nc1
1909	<< "\n calculated: " << c1
1910	<< "\n difference: " << std::fabs(nc1 - c1));
1911	}
1912
1913	const Real nc2 = NumericalDifferentiation(
1914	ext::function<Real(Real)>(
1915	LogCharacteristicFunction(2, t, cosEngine)),
1916	2, 1e-2, 5, central)(0.0);
1917
1918	const Real c2 = cosEngine->c2(t);
1919
1920	if (std::fabs(x: nc2 - c2) > tol) {
1921	BOOST_ERROR(" failed to reproduce second cumulant"
1922	<< "\n expected: " << nc2
1923	<< "\n calculated: " << c2
1924	<< "\n difference: " << std::fabs(nc2 - c2));
1925	}
1926
1927	const Real nc3 = NumericalDifferentiation(
1928	ext::function<Real(Real)>(
1929	LogCharacteristicFunction(3, t, cosEngine)),
1930	3, 5e-3, 7, central)(0.0);
1931
1932	const Real c3 = cosEngine->c3(t);
1933
1934	if (std::fabs(x: nc3 - c3) > tol) {
1935	BOOST_ERROR(" failed to reproduce third cumulant"
1936	<< "\n expected: " << nc3
1937	<< "\n calculated: " << c3
1938	<< "\n difference: " << std::fabs(nc3 - c3));
1939	}
1940
1941	const Real nc4 = NumericalDifferentiation(
1942	ext::function<Real(Real)>(
1943	LogCharacteristicFunction(4, t, cosEngine)),
1944	4, 5e-2, 9, central)(0.0);
1945
1946	const Real c4 = cosEngine->c4(t);
1947
1948	if (std::fabs(x: nc4 - c4) > 10*tol) {
1949	BOOST_ERROR(" failed to reproduce 4th cumulant"
1950	<< "\n expected: " << nc4
1951	<< "\n calculated: " << c4
1952	<< "\n difference: " << std::fabs(nc4 - c4));
1953	}
1954	}
1955	}
1956
1957	void HestonModelTest::testCosHestonEngine() {
1958	BOOST_TEST_MESSAGE("Testing Heston pricing via COS method...");
1959
1960	const Date settlementDate(7, February, 2017);
1961	Settings::instance().evaluationDate() = settlementDate;
1962
1963	const DayCounter dayCounter = Actual365Fixed();
1964	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.15, dc: dayCounter));
1965	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.07, dc: dayCounter));
1966
1967	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
1968
1969	const Real v0 = 0.1;
1970	const Real rho = -0.75;
1971	const Real sigma = 1.8;
1972	const Real kappa = 4.0;
1973	const Real theta = 0.22;
1974
1975	const ext::shared_ptr<HestonModel> model =
1976	ext::make_shared<HestonModel>(
1977	args: ext::make_shared<HestonProcess>(
1978	args: riskFreeTS, args: dividendTS,
1979	args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1980
1981	const Date maturityDate = settlementDate + Period(1, Years);
1982
1983	const ext::shared_ptr<Exercise> exercise =
1984	ext::make_shared<EuropeanExercise>(args: maturityDate);
1985
1986	const ext::shared_ptr<PricingEngine> cosEngine(
1987	ext::make_shared<COSHestonEngine>(args: model, args: 25, args: 600));
1988
1989	const ext::shared_ptr<StrikedTypePayoff> payoffs[] = {
1990	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0->value()+20),
1991	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0->value()+150),
1992	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: s0->value()-20),
1993	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: s0->value()-90)
1994	};
1995
1996	const Real expected[] = {
1997	9.364410588426075, 0.01036797658132471,
1998	5.319092971836708, 0.01032681906278383 };
1999
2000	const Real tol = 1e-10;
2001
2002	for (Size i=0; i < LENGTH(payoffs); ++i) {
2003	VanillaOption option(payoffs[i], exercise);
2004
2005	option.setPricingEngine(cosEngine);
2006	const Real calculated = option.NPV();
2007
2008	const Real error = std::fabs(x: expected[i] - calculated);
2009
2010	if (error > tol) {
2011	BOOST_ERROR(" failed to reproduce prices with COSHestonEngine"
2012	<< "\n expected: " << expected[i]
2013	<< "\n calculated: " << calculated
2014	<< "\n difference: " << error);
2015	}
2016	}
2017	}
2018
2019	void HestonModelTest::testCosHestonEngineTruncation() {
2020	BOOST_TEST_MESSAGE("Testing Heston pricing via COS method outside truncation bounds...");
2021
2022	const Date todaysDate(22, August, 2022);
2023	const Date maturity(23, August, 2022);
2024	Settings::instance().evaluationDate() = todaysDate;
2025
2026	Option::Type type(Option::Call);
2027	Real underlying = 100;
2028	Real strike = 200;
2029	Rate dividendYield = 0;
2030	Rate riskFreeRate = 0;
2031	DayCounter dayCounter = Actual365Fixed();
2032
2033	ext::shared_ptr<Exercise> europeanExercise(new EuropeanExercise(maturity));
2034	Handle<Quote> underlyingH(ext::shared_ptr<Quote>(new SimpleQuote(underlying)));
2035	Handle<YieldTermStructure> riskFreeTS(
2036	ext::shared_ptr<YieldTermStructure>(
2037	new FlatForward(todaysDate, riskFreeRate, dayCounter)));
2038	Handle<YieldTermStructure> dividendTS(
2039	ext::shared_ptr<YieldTermStructure>(
2040	new FlatForward(todaysDate, dividendYield, dayCounter)));
2041
2042	ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(type, strike));
2043	VanillaOption europeanOption(payoff, europeanExercise);
2044
2045	ext::shared_ptr<HestonProcess> hestonProcess(
2046	new HestonProcess(riskFreeTS, dividendTS, underlyingH, .007, .8, .007, .1, -.2));
2047	ext::shared_ptr<HestonModel> hestonModel(
2048	new HestonModel(hestonProcess));
2049
2050	europeanOption.setPricingEngine(ext::shared_ptr<PricingEngine>(
2051	new COSHestonEngine(hestonModel)));
2052
2053	const Real tol = 1e-7;
2054	const Real error = std::fabs(x: europeanOption.NPV() - 0.0);
2055
2056	if (error > tol) {
2057	BOOST_ERROR(" failed to reproduce prices with COSHestonEngine"
2058	<< "\n expected: " << 0.0
2059	<< "\n calculated: " << europeanOption.NPV()
2060	<< "\n difference: " << error);
2061	}
2062
2063	}
2064
2065	void HestonModelTest::testCharacteristicFct() {
2066	BOOST_TEST_MESSAGE("Testing Heston characteristic function...");
2067
2068	const Date settlementDate(30, March, 2017);
2069	Settings::instance().evaluationDate() = settlementDate;
2070
2071	const DayCounter dayCounter = Actual365Fixed();
2072	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.35, dc: dayCounter));
2073	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.17, dc: dayCounter));
2074
2075	const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
2076
2077	const Real v0 = 0.1;
2078	const Real rho = -0.85;
2079	const Real sigma = 0.8;
2080	const Real kappa = 2.0;
2081	const Real theta = 0.15;
2082
2083	const ext::shared_ptr<HestonModel> model =
2084	ext::make_shared<HestonModel>(
2085	args: ext::make_shared<HestonProcess>(
2086	args: riskFreeTS, args: dividendTS,
2087	args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
2088
2089	const Real u[] = { 1.0, 0.45, 3,4 };
2090	const Real t[] = { 0.01, 23.2, 3.2};
2091
2092	const COSHestonEngine cosEngine(model);
2093	const AnalyticHestonEngine analyticEngine(model);
2094
2095	constexpr double tol = 100*QL_EPSILON;
2096	for (Real i : u) {
2097	for (Real j : t) {
2098	const std::complex<Real> c = cosEngine.chF(u: i, t: j);
2099	const std::complex<Real> a = analyticEngine.chF(z: i, t: j);
2100
2101	const Real error = std::abs(z: a-c);
2102	if (error > tol) {
2103	BOOST_ERROR(" failed to reproduce prices with characteristic Fct"
2104	<< "\n Cos Engine: " << c
2105	<< "\n analytic engine: " << a
2106	<< "\n difference: " << error);
2107	}
2108	}
2109	}
2110	}
2111
2112	void HestonModelTest::testAndersenPiterbargPricing() {
2113	BOOST_TEST_MESSAGE("Testing Andersen-Piterbarg method to "
2114	"price under the Heston model...");
2115
2116	const Date settlementDate(30, March, 2017);
2117	Settings::instance().evaluationDate() = settlementDate;
2118
2119	const DayCounter dayCounter = Actual365Fixed();
2120	const Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.10, dc: dayCounter));
2121	const Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.06, dc: dayCounter));
2122
2123	const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
2124
2125	const Real v0 = 0.1;
2126	const Real rho = 0.80;
2127	const Real sigma = 0.75;
2128	const Real kappa = 1.0;
2129	const Real theta = 0.1;
2130
2131	const ext::shared_ptr<HestonModel> model =
2132	ext::make_shared<HestonModel>(
2133	args: ext::make_shared<HestonProcess>(
2134	args: riskFreeTS, args: dividendTS,
2135	args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
2136
2137	const ext::shared_ptr<AnalyticHestonEngine>
2138	andersenPiterbargLaguerreEngine(
2139	ext::make_shared<AnalyticHestonEngine>(
2140	args: model,
2141	args: AnalyticHestonEngine::AndersenPiterbarg,
2142	args: AnalyticHestonEngine::Integration::gaussLaguerre()));
2143
2144	const ext::shared_ptr<AnalyticHestonEngine>
2145	andersenPiterbargLobattoEngine(
2146	ext::make_shared<AnalyticHestonEngine>(
2147	args: model,
2148	args: AnalyticHestonEngine::AndersenPiterbarg,
2149	args: AnalyticHestonEngine::Integration::gaussLobatto(
2150	relTolerance: Null<Real>(), absTolerance: 1e-9, maxEvaluations: 10000), args: 1e-9));
2151
2152	const ext::shared_ptr<AnalyticHestonEngine>
2153	andersenPiterbargSimpsonEngine(
2154	ext::make_shared<AnalyticHestonEngine>(
2155	args: model,
2156	args: AnalyticHestonEngine::AndersenPiterbarg,
2157	args: AnalyticHestonEngine::Integration::discreteSimpson(evaluation: 256),
2158	args: 1e-8));
2159
2160	const ext::shared_ptr<AnalyticHestonEngine>
2161	andersenPiterbargTrapezoidEngine(
2162	ext::make_shared<AnalyticHestonEngine>(
2163	args: model,
2164	args: AnalyticHestonEngine::AndersenPiterbarg,
2165	args: AnalyticHestonEngine::Integration::discreteTrapezoid(evaluation: 164),
2166	args: 1e-8));
2167
2168	const ext::shared_ptr<AnalyticHestonEngine>
2169	andersenPiterbargTrapezoidEngine2(
2170	ext::make_shared<AnalyticHestonEngine>(
2171	args: model,
2172	args: AnalyticHestonEngine::AndersenPiterbarg,
2173	args: AnalyticHestonEngine::Integration::trapezoid(absTolerance: 1e-8, maxEvaluations: 256),
2174	args: 1e-8));
2175
2176	const ext::shared_ptr<ExponentialFittingHestonEngine>
2177	andersenPiterbargExponentialFittingEngine(
2178	ext::make_shared<ExponentialFittingHestonEngine>(args: model));
2179
2180	const ext::shared_ptr<PricingEngine> engines[] = {
2181	andersenPiterbargLaguerreEngine,
2182	andersenPiterbargLobattoEngine,
2183	andersenPiterbargSimpsonEngine,
2184	andersenPiterbargTrapezoidEngine,
2185	andersenPiterbargTrapezoidEngine2,
2186	andersenPiterbargExponentialFittingEngine
2187	};
2188
2189	const std::string algos[] = {
2190	"Gauss-Laguerre", "Gauss-Lobatto",
2191	"Discrete Simpson", "Discrete Trapezoid", "Trapezoid",
2192	"Exponential Fitting"
2193	};
2194
2195	const ext::shared_ptr<PricingEngine> analyticEngine(
2196	ext::make_shared<AnalyticHestonEngine>(args: model, args: 178));
2197
2198	const Date maturityDates[] = {
2199	settlementDate + Period(1, Days),
2200	settlementDate + Period(1, Weeks),
2201	settlementDate + Period(1, Years),
2202	settlementDate + Period(10, Years)
2203	};
2204
2205	const Option::Type optionTypes[] = { Option::Call, Option::Put };
2206	const Real strikes[] = { 50, 75, 90, 100, 110, 130, 150, 200};
2207
2208	const Real tol = 1e-7;
2209
2210	for (auto maturityDate : maturityDates) {
2211	const ext::shared_ptr<Exercise> exercise = ext::make_shared<EuropeanExercise>(args&: maturityDate);
2212
2213	for (auto optionType : optionTypes) {
2214	for (Real strike : strikes) {
2215	VanillaOption option(ext::make_shared<PlainVanillaPayoff>(args&: optionType, args&: strike),
2216	exercise);
2217
2218	option.setPricingEngine(analyticEngine);
2219	const Real expected = option.NPV();
2220
2221	for (Size k=0; k < LENGTH(engines); ++k) {
2222	option.setPricingEngine(engines[k]);
2223	const Real calculated = option.NPV();
2224
2225	const Real error = std::fabs(x: calculated-expected);
2226
2227	if (error > tol) {
2228	BOOST_ERROR(" failed to reproduce prices with Andersen-"
2229	"Piterbarg control variate"
2230	<< "\n algorithm : " << algos[k]
2231	<< "\n strike : " << strike
2232	<< "\n control variate: " << calculated
2233	<< "\n classic engine : " << expected
2234	<< "\n difference: " << error);
2235	}
2236	}
2237	}
2238	}
2239	}
2240	}
2241
2242
2243	void HestonModelTest::testAndersenPiterbargControlVariateIntegrand() {
2244	BOOST_TEST_MESSAGE("Testing Andersen-Piterbarg Integrand "
2245	"with control variate...");
2246
2247	const Date settlementDate(17, April, 2017);
2248	Settings::instance().evaluationDate() = settlementDate;
2249	const Date maturityDate = settlementDate + Period(2, Years);
2250
2251	const DayCounter dayCounter = Actual365Fixed();
2252	const Rate r = 0.075;
2253	const Rate q = 0.05;
2254	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc: dayCounter));
2255	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc: dayCounter));
2256
2257	const Time maturity = dayCounter.yearFraction(d1: settlementDate, d2: maturityDate);
2258	const DiscountFactor df = rTS->discount(t: maturity);
2259
2260	const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
2261	const Real fwd = s0->value()*qTS->discount(t: maturity)/df;
2262
2263	const Real strike = 150;
2264	const Real sx = std::log(x: strike);
2265	const Real dd = std::log(x: s0->value()*qTS->discount(t: maturity)/df);
2266
2267	const Real v0 = 0.08;
2268	const Real rho = -0.8;
2269	const Real sigma = 0.5;
2270	const Real kappa = 4.0;
2271	const Real theta = 0.05;
2272
2273	const ext::shared_ptr<HestonModel> hestonModel(
2274	ext::make_shared<HestonModel>(
2275	args: ext::make_shared<HestonProcess>(
2276	args: rTS, args: qTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
2277
2278	const ext::shared_ptr<COSHestonEngine> cosEngine(
2279	ext::make_shared<COSHestonEngine>(args: hestonModel));
2280
2281	const ext::shared_ptr<AnalyticHestonEngine> engine(
2282	ext::make_shared<AnalyticHestonEngine>(
2283	args: hestonModel,
2284	args: AnalyticHestonEngine::AndersenPiterbarg,
2285	args: AnalyticHestonEngine::Integration::gaussLaguerre()));
2286
2287	VanillaOption option(
2288	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
2289	ext::make_shared<EuropeanExercise>(args: maturityDate));
2290	option.setPricingEngine(engine);
2291
2292	const Real refNPV = option.NPV();
2293
2294	const Volatility implStdDev = blackFormulaImpliedStdDev(
2295	optionType: Option::Call, strike, forward: fwd, blackPrice: refNPV, discount: df);
2296
2297	const Real var = cosEngine->var(t: maturity);
2298	const Real stdDev = std::sqrt(x: var);
2299
2300	const Real d = (std::log(x: s0->value()/strike)
2301	+ (r-q)*maturity+ 0.5*var)/stdDev;
2302
2303	const Real skew = cosEngine->skew(t: maturity);
2304	const Real kurt = cosEngine->kurtosis(t: maturity);
2305
2306	const NormalDistribution n;
2307
2308	const Real q3 = 1/6.*s0->value()*stdDev*(2*stdDev - d)*n(d);
2309	const Real q4 = 1/24.*s0->value()*stdDev*(d*d - 3*d*stdDev - 1)*n(d);
2310	const Real q5 = 1/72.*s0->value()*stdDev*(
2311	d*d*d*d - 5*d*d*d*stdDev - 6*d*d + 15*d*stdDev + 3)*n(d);
2312
2313	const Real bsNPV = blackFormula(optionType: Option::Call, strike, forward: fwd, stdDev, discount: df);
2314
2315	// different variance values for the control variate
2316	const Real variances[] = {
2317	v0*maturity,
2318	((1-std::exp(x: -kappa*maturity))*(v0-theta)/(kappa*maturity) + theta)
2319	*maturity,
2320	// second moment as control variate
2321	var,
2322	// third and fourth moment pricing based on
2323	// Corrado C. and T. Su, (1996-b),
2324	// “Skewness and Kurtosis in S&P 500 IndexReturns Implied by Option Prices”,
2325	// Journal of Financial Research 19 (2), 175-192.
2326	squared(x: blackFormulaImpliedStdDev(
2327	optionType: Option::Call, strike, forward: fwd, blackPrice: bsNPV + skew*q3, discount: df)),
2328	squared(x: blackFormulaImpliedStdDev(
2329	optionType: Option::Call, strike, forward: fwd, blackPrice: bsNPV + skew*q3 + kurt*q4, discount: df)),
2330	// Moment matching based on
2331	// Rubinstein M., (1998), “Edgeworth Binomial Trees”,
2332	// Journal of Derivatives 5 (3), 20-27.
2333	squared(x: blackFormulaImpliedStdDev(
2334	optionType: Option::Call, strike, forward: fwd,
2335	blackPrice: bsNPV + skew*q3 + kurt*q4 + skew*skew*q5, discount: df)),
2336	// implied vol as control variate
2337	squared(x: implStdDev),
2338	// remaining function becomes zero for u -> 0
2339	-8.0*std::log(x: engine->chF(z: std::complex<Real>(0, -0.5), t: maturity).real())
2340	};
2341
2342	for (Size i=0; i < LENGTH(variances); ++i) {
2343	const Real sigmaBS = std::sqrt(x: variances[i]/maturity);
2344
2345	for (Real u =0.001; u < 15; u*=1.05) {
2346	const std::complex<Real> z(u, -0.5);
2347
2348	const std::complex<Real> phiBS
2349	= std::exp(z: -0.5*sigmaBS*sigmaBS*maturity
2350	*(z*z + std::complex<Real>(-z.imag(), z.real())));
2351
2352	const std::complex<Real> ex
2353	= std::exp(z: std::complex<Real>(0.0, u*(dd-sx)));
2354
2355	const std::complex<Real> chf = engine->chF(z, t: maturity);
2356
2357	const Real orig = (-ex*chf / (u*u + 0.25)).real();
2358	const Real cv = (ex*(phiBS - chf) / (u*u + 0.25)).real();
2359
2360	if (std::fabs(x: cv) > 0.03) {
2361	BOOST_ERROR(" Control variate function is greater "
2362	"than original function"
2363	<< "\n control variate method : " << i
2364	<< "\n z value : " << u
2365	<< "\n control variate function: " << cv
2366	<< "\n original function : " << orig);
2367	}
2368	}
2369	}
2370	}
2371
2372	void HestonModelTest::testAndersenPiterbargConvergence() {
2373	BOOST_TEST_MESSAGE("Testing Andersen-Piterbarg pricing convergence...");
2374
2375	const Date settlementDate(5, July, 2002);
2376	Settings::instance().evaluationDate() = settlementDate;
2377	const Date maturityDate(5, July, 2003);
2378
2379	const DayCounter dayCounter = Actual365Fixed();
2380	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.01, dc: dayCounter));
2381	const Handle<YieldTermStructure> qTS(flatRate(forward: 0.02, dc: dayCounter));
2382
2383	const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
2384
2385	const Real v0 = 0.04;
2386	const Real rho = -0.5;
2387	const Real sigma = 1.0;
2388	const Real kappa = 4.0;
2389	const Real theta = 0.25;
2390
2391	const ext::shared_ptr<HestonModel> hestonModel(
2392	ext::make_shared<HestonModel>(
2393	args: ext::make_shared<HestonProcess>(
2394	args: rTS, args: qTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
2395
2396	VanillaOption option(
2397	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0->value()),
2398	ext::make_shared<EuropeanExercise>(args: maturityDate));
2399
2400
2401	// Alan Lewis reference prices posted in
2402	// http://wilmott.com/messageview.cfm?catid=34&threadid=90957
2403	const Real reference = 16.070154917028834278213466703938231827658768230714;
2404
2405	const Real diffs[] = {
2406	0.0892433814611486298, 0.00013096156482816923,
2407	1.34107015270501506e-07, 1.22913235145460931e-10,
2408	1.24344978758017533e-13 };
2409
2410	for (Size n=10; n <= 50; n+=10) {
2411	option.setPricingEngine(ext::make_shared<AnalyticHestonEngine>(
2412	args: hestonModel, args: AnalyticHestonEngine::AndersenPiterbarg,
2413	args: AnalyticHestonEngine::Integration::discreteTrapezoid(evaluation: n), args: 1e-13));
2414
2415	const Real calculatedDiff = std::fabs(x: option.NPV()-reference);
2416	if (calculatedDiff > 1.25*diffs[n/10-1])
2417	BOOST_ERROR("failed to prove convergence for trapezoid rule "
2418	<< "\n calculated difference: " << calculatedDiff
2419	<< "\n expected difference: " << diffs[n/10-1]);
2420	}
2421	}
2422
2423
2424	void HestonModelTest::testPiecewiseTimeDependentChFvsHestonChF() {
2425	BOOST_TEST_MESSAGE("Testing piecewise time dependent "
2426	"ChF vs Heston ChF...");
2427
2428	const Date settlementDate(5, July, 2017);
2429	Settings::instance().evaluationDate() = settlementDate;
2430	const Date maturityDate(5, July, 2018);
2431
2432	const DayCounter dayCounter = Actual365Fixed();
2433	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.01, dc: dayCounter));
2434	const Handle<YieldTermStructure> qTS(flatRate(forward: 0.02, dc: dayCounter));
2435
2436	const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
2437
2438	const Real v0 = 0.04;
2439	const Real rho = -0.5;
2440	const Real sigma = 1.0;
2441	const Real kappa = 4.0;
2442	const Real theta = 0.25;
2443
2444	const ConstantParameter thetaP(theta, PositiveConstraint());
2445	const ConstantParameter kappaP(kappa, PositiveConstraint());
2446	const ConstantParameter sigmaP(sigma, PositiveConstraint());
2447	const ConstantParameter rhoP (rho, BoundaryConstraint(-1.0, 1.0));
2448
2449	const ext::shared_ptr<AnalyticHestonEngine> analyticEngine(
2450	ext::make_shared<AnalyticHestonEngine>(
2451	args: ext::make_shared<HestonModel>(
2452	args: ext::make_shared<HestonProcess>(
2453	args: rTS, args: qTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho))));
2454
2455	const ext::shared_ptr<AnalyticPTDHestonEngine> ptdHestonEngine(
2456	ext::make_shared<AnalyticPTDHestonEngine>(
2457	args: ext::make_shared<PiecewiseTimeDependentHestonModel>(
2458	args: rTS, args: qTS, args: s0, args: v0, args: thetaP, args: kappaP, args: sigmaP, args: rhoP,
2459	args: TimeGrid(dayCounter.yearFraction(d1: settlementDate, d2: maturityDate),
2460	10))));
2461
2462	constexpr double tol = 100 * QL_EPSILON;
2463	for (Real r = 0.1; r < 4; r+=0.25) {
2464	for (Real phi = 0; phi < 360; phi+=60) {
2465	for (Time t=0.1; t <= 1.0; t+=0.3) {
2466	const std::complex<Real> z
2467	= r*std::exp(z: std::complex<Real>(0, phi));
2468
2469	const std::complex<Real> a = analyticEngine->chF(z, t);
2470	const std::complex<Real> b = ptdHestonEngine->chF(z, t);
2471
2472	if (std::abs(z: a-b) > tol)
2473	BOOST_ERROR("failed to compare characteristic function "
2474	<< "\n time dependent model: " << b
2475	<< "\n Heston model : " << a
2476	<< "\n Difference : " << std::abs(a-b));
2477	}
2478	}
2479	}
2480	}
2481
2482
2483	void HestonModelTest::testPiecewiseTimeDependentComparison() {
2484	BOOST_TEST_MESSAGE("Testing piecewise time dependent "
2485	"ChF vs Heston ChF...");
2486
2487	const Date settlementDate(5, July, 2017);
2488	Settings::instance().evaluationDate() = settlementDate;
2489
2490	const DayCounter dc = Actual365Fixed();
2491	const Date maturityDate(5, July, 2018);
2492	const Time maturity = dc.yearFraction(d1: settlementDate, d2: maturityDate);
2493
2494	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc));
2495	const Handle<YieldTermStructure> qTS(flatRate(forward: 0.08, dc));
2496
2497	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
2498
2499	std::vector<Time> modelTimes = {0.25, 0.75, 10.0};
2500	const TimeGrid modelGrid(modelTimes.begin(), modelTimes.end());
2501
2502	const Real v0 = 0.1;
2503	ConstantParameter theta( 0.1, PositiveConstraint());
2504	ConstantParameter kappa( 1.0, PositiveConstraint());
2505	ConstantParameter rho( -0.75, BoundaryConstraint(-1.0, 1.0));
2506
2507	std::vector<Time> pTimes(2);
2508	pTimes[0] = 0.25;
2509	pTimes[1] = 0.75;
2510	PiecewiseConstantParameter sigma(pTimes, PositiveConstraint());
2511
2512	sigma.setParam(i: 0, x: 0.30);
2513	sigma.setParam(i: 1, x: 0.15);
2514	sigma.setParam(i: 2, x: 1.25);
2515
2516	VanillaOption option(
2517	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0->value()),
2518	ext::make_shared<EuropeanExercise>(args: maturityDate));
2519
2520	const ext::shared_ptr<PiecewiseTimeDependentHestonModel> ptdModel(
2521	ext::make_shared<PiecewiseTimeDependentHestonModel>(
2522	args: rTS, args: qTS, args: s0, args: v0, args&: theta, args&: kappa, args&: sigma, args&: rho, args: modelGrid));
2523
2524	const ext::shared_ptr<AnalyticPTDHestonEngine> ptdHestonEngine(
2525	ext::make_shared<AnalyticPTDHestonEngine>(args: ptdModel));
2526
2527	option.setPricingEngine(ptdHestonEngine);
2528	const Real calculatedGatheral = option.NPV();
2529
2530	const ext::shared_ptr<AnalyticPTDHestonEngine> ptdAPEngine(
2531	ext::make_shared<AnalyticPTDHestonEngine>(
2532	args: ptdModel,
2533	args: AnalyticPTDHestonEngine::AndersenPiterbarg,
2534	args: AnalyticPTDHestonEngine::Integration::discreteTrapezoid(evaluation: 128),
2535	args: 1e-12));
2536	option.setPricingEngine(ptdAPEngine);
2537	const Real calculatedAndersenPiterbarg = option.NPV();
2538
2539	if (std::fabs(x: calculatedGatheral - calculatedAndersenPiterbarg) > 1e-10)
2540	BOOST_ERROR("failed to reproduce npv for time dependent Heston model "
2541	<< "\n Gatheral ChF : " << calculatedGatheral
2542	<< "\n AndersenPiterbarg ChF: " << calculatedAndersenPiterbarg
2543	<< "\n Difference : "
2544	<< std::fabs(calculatedGatheral - calculatedAndersenPiterbarg));
2545
2546	const ext::shared_ptr<HestonProcess> firstPartProcess(
2547	ext::make_shared<HestonProcess>(
2548	args: rTS, args: qTS, args: s0, args: v0, args: 1.0, args: 0.1, args: 0.30, args: -0.75,
2549	args: HestonProcess::QuadraticExponentialMartingale));
2550
2551	typedef PseudoRandom::rsg_type rsg_type;
2552	typedef PseudoRandom::urng_type urng_type;
2553	typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
2554
2555	const MultiPathGenerator<rsg_type> firstPathGen(
2556	firstPartProcess,
2557	TimeGrid(pTimes.front(), 6),
2558	PseudoRandom::make_sequence_generator(dimension: 12, seed: 1234));
2559
2560	const urng_type urng(5678);
2561
2562	Statistics stat;
2563	const DiscountFactor df = rTS->discount(d: maturityDate);
2564
2565	const Size nSims = 10000;
2566	for (Size i=0; i < nSims; ++i) {
2567	Real priceS = 0.0;
2568
2569	for (Size j=0; j < 2; ++j) {
2570	const sample_type& path1 =
2571	(j & 1) != 0U ? firstPathGen.antithetic() : firstPathGen.next();
2572	const Real spot1 = path1.value[0].back();
2573	const Real v1 = path1.value[1].back();
2574
2575	const MultiPathGenerator<rsg_type> secondPathGen(
2576	ext::make_shared<HestonProcess>(
2577	args: rTS, args: qTS,
2578	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: spot1)),
2579	args: v1, args: 1.0, args: 0.1, args: 0.15, args: -0.75,
2580	args: HestonProcess::QuadraticExponentialMartingale),
2581	TimeGrid(pTimes[1]-pTimes[0], 12),
2582	PseudoRandom::make_sequence_generator(dimension: 24, seed: urng.nextInt32()));
2583
2584	const sample_type& path2 = secondPathGen.next();
2585	const Real spot2 = path2.value[0].back();
2586	const Real v2 = path2.value[1].back();
2587
2588	const MultiPathGenerator<rsg_type> thirdPathGen(
2589	ext::make_shared<HestonProcess>(
2590	args: rTS, args: qTS,
2591	args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: spot2)),
2592	args: v2, args: 1.0, args: 0.1, args: 1.25, args: -0.75,
2593	args: HestonProcess::QuadraticExponentialMartingale),
2594	TimeGrid(maturity-pTimes[1], 6),
2595	PseudoRandom::make_sequence_generator(dimension: 12, seed: urng.nextInt32()));
2596	const sample_type& path3 = thirdPathGen.next();
2597	const Real spot3 = path3.value[0].back();
2598
2599	priceS += 0.5*(*option.payoff())(spot3);
2600	}
2601
2602	stat.add(value: priceS*df);
2603	}
2604
2605	const Real calculatedMC = stat.mean();
2606	const Real errorEstimate = stat.errorEstimate();
2607
2608	if (std::fabs(x: calculatedMC - calculatedGatheral) > 3.0*errorEstimate)
2609	BOOST_ERROR("failed to reproduce npv for time dependent Heston model"
2610	<< "\n Gatheral ChF : " << calculatedGatheral
2611	<< "\n Monte-Carlo : " << calculatedMC
2612	<< "\n Monte-Carlo error: " << errorEstimate
2613	<< "\n Difference : "
2614	<< std::fabs(calculatedGatheral - calculatedMC));
2615	}
2616
2617	void HestonModelTest::testPiecewiseTimeDependentChFAsymtotic() {
2618	BOOST_TEST_MESSAGE("Testing piecewise time dependent "
2619	"ChF vs Heston ChF...");
2620
2621	const Date settlementDate(5, July, 2017);
2622	Settings::instance().evaluationDate() = settlementDate;
2623	const Date maturityDate = settlementDate + Period(13, Months);
2624
2625	const DayCounter dc = Actual365Fixed();
2626	const Time maturity = dc.yearFraction(d1: settlementDate, d2: maturityDate);
2627	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.0, dc));
2628
2629	std::vector<Time> modelTimes = {0.01, 0.5, 2.0};
2630
2631	const TimeGrid modelGrid(modelTimes.begin(), modelTimes.end());
2632
2633	const Real v0 = 0.1;
2634	const std::vector<Time> pTimes(modelTimes.begin(), modelTimes.end()-1);
2635
2636	PiecewiseConstantParameter sigma(pTimes, PositiveConstraint());
2637	PiecewiseConstantParameter theta(pTimes, PositiveConstraint());
2638	PiecewiseConstantParameter kappa(pTimes, PositiveConstraint());
2639	PiecewiseConstantParameter rho(pTimes, BoundaryConstraint(-1.0, 1.0));
2640
2641	const Real sigmas[] = { 0.01, 0.2, 0.6 };
2642	const Real thetas[] = { 0.16, 0.06, 0.36 };
2643	const Real kappas[] = { 1.0, 0.3, 4.0 };
2644	const Real rhos[] = { 0.5, -0.75, -0.25 };
2645
2646	for (Size i=0; i < 3; ++i) {
2647	sigma.setParam(i, x: sigmas[i]);
2648	theta.setParam(i, x: thetas[i]);
2649	kappa.setParam(i, x: kappas[i]);
2650	rho.setParam(i, x: rhos[i]);
2651	}
2652
2653	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
2654	const ext::shared_ptr<PiecewiseTimeDependentHestonModel> ptdModel(
2655	ext::make_shared<PiecewiseTimeDependentHestonModel>(
2656	args: rTS, args: rTS, args: s0, args: v0, args&: theta, args&: kappa, args&: sigma, args&: rho, args: modelGrid));
2657
2658	const Real eps = 1e-8;
2659
2660	const ext::shared_ptr<AnalyticPTDHestonEngine> ptdHestonEngine(
2661	ext::make_shared<AnalyticPTDHestonEngine>(
2662	args: ptdModel,
2663	args: AnalyticPTDHestonEngine::AndersenPiterbarg,
2664	args: AnalyticPTDHestonEngine::Integration::discreteTrapezoid(evaluation: 128),
2665	args: eps));
2666
2667	const std::complex<Real> D_u_inf = -
2668	std::complex<Real>(std::sqrt(x: 1-rhos[0]*rhos[0]),rhos[0])/sigmas[0];
2669
2670	const std::complex<Real> dd = std::complex<Real>(kappas[0],
2671	(2*kappas[0]*rhos[0]-sigmas[0])
2672	/(2*std::sqrt(x: 1-rhos[0]*rhos[0])))/(sigmas[0]*sigmas[0]);
2673
2674	std::complex<Real> C_u_inf(0.0, 0.0), cc(0.0, 0.0), clog(0.0, 0.0);
2675
2676	for (Size i=0; i < 3; ++i) {
2677	const Real kappa = kappas[i];
2678	const Real theta = thetas[i];
2679	const Real sigma = sigmas[i];
2680	const Real rho = rhos[i];
2681	const Time tau = std::min(a: maturity, b: modelGrid[i+1]) - modelGrid[i];
2682
2683	C_u_inf += -kappa*theta*tau / sigma
2684	*std::complex<Real>(std::sqrt(x: 1-rho*rho), rho);
2685
2686	cc += kappa*std::complex<Real>(2*kappa,(2*kappa*rho-sigma)
2687	/sqrt(x: 1-rho*rho))*tau*theta/(2*sigma*sigma);
2688
2689	const std::complex<Real> Di =
2690	(i < 2) ? sigma/sigmas[i+1]
2691	*std::complex<Real>(std::sqrt(x: 1-rhos[i+1]*rhos[i+1]), rhos[i+1])
2692	: std::complex<Real>(0.0, 0.0);
2693
2694	clog += 2*kappa*theta/(sigma*sigma)*std::log(z: 1.0 -
2695	( Di - std::complex<Real>(std::sqrt(x: 1-rho*rho), rho)) /
2696	( Di + std::complex<Real>(std::sqrt(x: 1-rho*rho), -rho)));
2697	}
2698
2699	const Real epsilon = eps*M_PI/s0->value();
2700
2701	const Real uM =
2702	AnalyticHestonEngine::Integration::andersenPiterbargIntegrationLimit(
2703	c_inf: -(C_u_inf + D_u_inf*v0).real(), epsilon, v0, t: maturity);
2704
2705	const Real expectedUM = 18.6918883427;
2706	if (std::fabs(x: uM - expectedUM) > 1e-5) {
2707	BOOST_ERROR("failed to reproduce Andersen-Piterbarg "
2708	"Integration bounds for piecewise constant "
2709	"time dependent Heston Model"
2710	<< "\n calculated : " << uM
2711	<< "\n expected : " << expectedUM
2712	<< "\n diff : " << std::fabs(uM - expectedUM)
2713	<< "\n tolerance : " << 1e-5);
2714	}
2715
2716	const Real u = 1e8;
2717	const std::complex<Real> expectedlnChF = ptdHestonEngine->lnChF(z: u, t: maturity);
2718	const std::complex<Real> calculatedAsympotic =
2719	(D_u_inf*u + dd)*v0 + C_u_inf*u + cc + clog;
2720
2721	if (std::abs(z: expectedlnChF - calculatedAsympotic) > 0.01) {
2722	BOOST_ERROR("failed to reproduce asymptotic of characteristic function"
2723	<< "\n ln(ChF) : " << expectedlnChF
2724	<< "\n asymptotic: " << calculatedAsympotic
2725	<< "\n diff : "
2726	<< std::abs(expectedlnChF - calculatedAsympotic)
2727	<< "\n tolerance : " << 0.01);
2728	}
2729
2730	VanillaOption option(
2731	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0->value()),
2732	ext::make_shared<EuropeanExercise>(args: maturityDate));
2733	option.setPricingEngine(ptdHestonEngine);
2734
2735	const Real expectedNPV = 17.43851162589377;
2736	const Real calculatedNPV = option.NPV();
2737	const Real diffNPV = std::fabs(x: expectedNPV - calculatedNPV);
2738	if (diffNPV > 1e-9) {
2739	BOOST_ERROR("failed to reproduce high precision prices for "
2740	"piecewise constant time dependent Heston model"
2741	<< "\n expeceted : " << expectedNPV
2742	<< "\n calclated : " << calculatedNPV
2743	<< "\n diff : " << diffNPV
2744	<< "\n tolerance : " << 1e-9);
2745	}
2746	}
2747
2748	void HestonModelTest::testSmallSigmaExpansion() {
2749	BOOST_TEST_MESSAGE("Testing small sigma expansion of "
2750	"the characteristic function...");
2751
2752	const Date settlementDate(20, March, 2020);
2753	Settings::instance().evaluationDate() = settlementDate;
2754	const Date maturityDate = settlementDate + Period(2, Years);
2755
2756	const DayCounter dc = Actual365Fixed();
2757	const Time t = dc.yearFraction(d1: settlementDate, d2: maturityDate);
2758	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.0, dc));
2759
2760	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: 100));
2761
2762	const Real theta = 0.1 * 0.1;
2763	const Real v0 = theta + 0.02;
2764	const Real kappa = 1.25;
2765	const Real sigma = 1e-9;
2766	const Real rho = -0.9;
2767
2768	const ext::shared_ptr<HestonModel> hestonModel =
2769	ext::make_shared<HestonModel>(
2770	args: ext::make_shared<HestonProcess>(
2771	args: rTS, args: rTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
2772
2773	const ext::shared_ptr<AnalyticHestonEngine> engine =
2774	ext::make_shared<AnalyticHestonEngine>(args: hestonModel);
2775
2776	const std::complex<Real> expectedChF(
2777	0.990463578538352651,2.60693475987521132e-12);
2778
2779	const std::complex<Real> calculatedChF = engine->chF(
2780	z: std::complex<Real>(0.55, -0.5), t);
2781
2782	const Real diffChF = std::abs(z: expectedChF - calculatedChF);
2783	const Real tolChF = 1e-12;
2784	if (diffChF > tolChF) {
2785	BOOST_ERROR("failed to reproduce normalized characteristic function "
2786	"value for small sigma"
2787	<< "\n expeceted : " << expectedChF
2788	<< "\n calclated : " << calculatedChF
2789	<< "\n diff : " << diffChF
2790	<< "\n tolerance : " << tolChF);
2791	}
2792
2793	VanillaOption option(
2794	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 120.0),
2795	ext::make_shared<EuropeanExercise>(args: maturityDate));
2796
2797	option.setPricingEngine(
2798	ext::make_shared<AnalyticHestonEngine>(
2799	args: hestonModel,
2800	args: AnalyticHestonEngine::AndersenPiterbarg,
2801	args: AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder: 192)));
2802
2803	const Real calculatedNPV = option.NPV();
2804
2805	const Real stdDev =
2806	std::sqrt(x: ((1-std::exp(x: -kappa*t))*(v0-theta)/(kappa*t) + theta)*t);
2807
2808	const Real expectedNPV =
2809	blackFormula(optionType: Option::Call, strike: 120.0, forward: 100.0, stdDev);
2810
2811	const Real diffNPV =std::fabs(x: calculatedNPV - expectedNPV);
2812	const Real tolNPV = 50*sigma;
2813
2814	if (diffNPV > tolNPV) {
2815	BOOST_ERROR("failed to reproduce Black Scholes prices "
2816	"for Heston model with very small sigma"
2817	<< "\n expeceted : " << expectedNPV
2818	<< "\n calclated : " << calculatedNPV
2819	<< "\n diff : " << diffNPV
2820	<< "\n tolerance : " << tolNPV);
2821	}
2822	}
2823
2824	void HestonModelTest::testSmallSigmaExpansion4ExpFitting() {
2825	BOOST_TEST_MESSAGE("Testing small sigma expansion for the "
2826	"exponential fitting Heston engine...");
2827
2828	const Date todaysDate(13, March, 2020);
2829	Settings::instance().evaluationDate() = todaysDate;
2830
2831	const DayCounter dc = Actual365Fixed();
2832	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc));
2833	const Handle<YieldTermStructure> qTS(flatRate(forward: 0.075, dc));
2834
2835	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: 100.0));
2836
2837	// special case: reduce sigma
2838	const Date maturityDate = Date(14, March, 2021);
2839	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
2840	const Real fwd =
2841	spot->value()*qTS->discount(t: maturity)/rTS->discount(t: maturity);
2842
2843	const Real v0 = 0.04;
2844	const Real rho = -0.5;
2845	const Real kappa = 4.0;
2846	const Real theta = 0.04;
2847
2848	const Real moneyness = 0.1;
2849	const Real strike = std::exp(x: -moneyness*std::sqrt(x: theta*maturity))*fwd;
2850
2851	const Real expected = blackFormula(
2852	optionType: Option::Call, strike, forward: fwd,
2853	stdDev: std::sqrt(x: v0*maturity), discount: rTS->discount(t: maturity));
2854
2855	VanillaOption option(
2856	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
2857	ext::make_shared<EuropeanExercise>(args: maturityDate));
2858
2859	for (Real sigma = 1e-4; sigma > 1e-12; sigma*=0.1) {
2860	option.setPricingEngine(
2861	ext::make_shared<ExponentialFittingHestonEngine>(
2862	args: ext::make_shared<HestonModel>(
2863	args: ext::make_shared<HestonProcess>(
2864	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args&: sigma, args: rho))));
2865	const Real calculated = option.NPV();
2866
2867	const Real diff = std::fabs(x: expected - calculated);
2868
2869	if (diff > 0.01*sigma) {
2870	BOOST_ERROR("failed to reproduce Black Scholes prices "
2871	"for Heston model with very small sigma"
2872	<< "\n expeceted : " << expected
2873	<< "\n calclated : " << calculated
2874	<< "\n sigma : " << sigma
2875	<< "\n diff : " << diff
2876	<< "\n tolerance : " << 10*sigma);
2877	}
2878	}
2879
2880
2881	// generic cases
2882	const Real kappas[] = { 0.5, 1.0, 4.0 };
2883	const Real thetas[] = { 0.04, 0.09};
2884	const Real v0s[] = { 0.025, 0.20 };
2885	const Integer maturities[] = { 1, 31, 182, 1850 };
2886
2887	for (int maturitie : maturities) {
2888	const Date maturityDate = todaysDate + Period(maturitie, Days);
2889	const DiscountFactor df = rTS->discount(d: maturityDate);
2890	const Real fwd = spot->value() * qTS->discount(d: maturityDate)/df;
2891
2892	const ext::shared_ptr<Exercise> exercise =
2893	ext::make_shared<EuropeanExercise>(args: maturityDate);
2894
2895	const Time t = dc.yearFraction(d1: todaysDate, d2: maturityDate);
2896
2897	Option::Type optionType = Option::Call;
2898
2899	for (Real kappa : kappas) {
2900	for (Real theta : thetas) {
2901	for (Real v0 : v0s) {
2902	const ext::shared_ptr<PricingEngine> engine =
2903	ext::make_shared<ExponentialFittingHestonEngine>(
2904	args: ext::make_shared<HestonModel>(args: ext::make_shared<HestonProcess>(
2905	args: rTS, args: qTS, args: spot, args&: v0, args&: kappa, args&: theta, args: 1e-13, args: -0.8)));
2906
2907	const Real stdDev =
2908	std::sqrt(x: ((1-std::exp(x: -kappa*t))*(v0-theta)/(kappa*t) + theta)*t);
2909
2910	for (Real strike = spot->value()*exp(x: -10*stdDev);
2911	strike < spot->value()*exp(x: 10*stdDev); strike*= 1.2) {
2912
2913	VanillaOption option(
2914	ext::make_shared<PlainVanillaPayoff>(
2915	args&: optionType, args&: strike), exercise);
2916
2917	option.setPricingEngine(engine);
2918	const Real calculated = option.NPV();
2919
2920	const Real expected =
2921	blackFormula(optionType, strike, forward: fwd, stdDev, discount: df);
2922
2923	const Real diff = std::fabs(x: expected - calculated);
2924	if (diff > 1e-10) {
2925	BOOST_ERROR("failed to reproduce Black Scholes prices "
2926	"for Heston model with very small sigma"
2927	<< "\n expceted : " << expected
2928	<< "\n calculated: " << calculated
2929	<< "\n diff : " << diff
2930	<< "\n tolerance : " << 1e-10);
2931	}
2932
2933	optionType = (optionType == Option::Call)
2934	? Option::Put : Option::Call;
2935	}
2936	}
2937	}
2938	}
2939	}
2940	}
2941
2942	void HestonModelTest::testExponentialFitting4StrikesAndMaturities() {
2943	BOOST_TEST_MESSAGE("Testing exponential fitting Heston engine "
2944	"with high precision results for large moneyness...");
2945
2946	const Date todaysDate = Date(13, May, 2020);
2947	Settings::instance().evaluationDate() = todaysDate;
2948
2949	const DayCounter dc = Actual365Fixed();
2950
2951	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.0507, dc));
2952	const Handle<YieldTermStructure> qTS(flatRate(forward: 0.0469, dc));
2953
2954	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.0));
2955
2956	const Real moneyness[] = { -20, -10, -5, 2.5, 1, 0, 1, 2.5, 5, 10, 20 };
2957	const Period maturities[] = {
2958	Period(1, Days),
2959	Period(1, Months),
2960	Period(1, Years),
2961	Period(10, Years)
2962	};
2963
2964	const Real v0 = 0.04;
2965	const Real rho = -0.6;
2966	const Real sigma = 0.75;
2967	const Real kappa = 2.5;
2968	const Real theta = 0.06;
2969
2970	// Reference prices are calculated using a boost multi-precision
2971	// implementation of the AnalyticHestonEngine,
2972	// https://github.com/klausspanderen/HestonExponentialFitting
2973
2974	const Real referenceValues[] = {
2975	1.1631865252540813e-58,
2976	1.06426822273258466e-49,
2977	6.92896489110422086e-16,
2978	8.19515526286263236e-06,
2979	0.000625608178476390504,
2980	0.00417261379371945684,
2981	0.000625608178476390504,
2982	8.19515526286263236e-06,
2983	1.92308901296741414e-10,
2984	1.57327901822368115e-23,
2985	5.7830515043285098e-58,
2986	3.56081886910098813e-48,
2987	2.9489071194212509e-23,
2988	1.54181757781090727e-11,
2989	0.000367960011879847279,
2990	0.00493886106106039818,
2991	0.0227152343265593776,
2992	0.00493886106106039818,
2993	0.000367960011879847279,
2994	3.06653474407784574e-06,
2995	8.86665241279348934e-11,
2996	1.51206812371708868e-20,
2997	4.18506719865401643e-29,
2998	2.46637786897559908e-15,
2999	1.75338784910563671e-08,
3000	0.00284789176080218294,
3001	0.0199133097064688458,
3002	0.0776848755698912041,
3003	0.0199133097064688458,
3004	0.00284789176080218294,
3005	0.00012462190796343504,
3006	2.59755319566692257e-07,
3007	1.13853114743124721e-12,
3008	4.27612073892114211e-39,
3009	1.08387452075906664e-25,
3010	4.15179522944463802e-11,
3011	0.00134157732880653131,
3012	0.029018582813884912,
3013	0.176405213088554197,
3014	0.029018582813884912,
3015	0.00134157732880653131,
3016	5.43674074281991917e-06,
3017	6.51443921040230507e-11,
3018	9.25756999394709285e-21
3019	};
3020
3021	const ext::shared_ptr<HestonModel> model =
3022	ext::make_shared<HestonModel>(
3023	args: ext::make_shared<HestonProcess>(
3024	args: rTS, args: qTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
3025
3026	const ext::shared_ptr<PricingEngine> engine =
3027	ext::make_shared<ExponentialFittingHestonEngine>(args: model);
3028
3029	Size idx = 0;
3030	for (auto maturitie : maturities) {
3031	const Date maturityDate = todaysDate + maturitie;
3032	const Time t = dc.yearFraction(d1: todaysDate, d2: maturityDate);
3033
3034	const ext::shared_ptr<Exercise> exercise =
3035	ext::make_shared<EuropeanExercise>(args: maturityDate);
3036
3037	const DiscountFactor df = rTS->discount(t);
3038	const Real fwd = s0->value()*qTS->discount(t)/df;
3039
3040	for (Size j=0; j < LENGTH(moneyness); ++j, ++idx) {
3041	const Real strike =
3042	std::exp(x: -moneyness[j]*std::sqrt(x: theta*t))*fwd;
3043
3044	for (Size k=0; k < 2; ++k) {
3045	const ext::shared_ptr<PlainVanillaPayoff> payoff =
3046	ext::make_shared<PlainVanillaPayoff>(args: (k) != 0U ? Option::Put : Option::Call,
3047	args: strike);
3048
3049	VanillaOption option(payoff, exercise);
3050	option.setPricingEngine(engine);
3051
3052	const Real calculated = option.NPV();
3053
3054	Real expected;
3055	if (payoff->optionType() == Option::Call)
3056	if (fwd < strike)
3057	expected = referenceValues[idx];
3058	else
3059	expected = (fwd - strike)*df + referenceValues[idx];
3060	else
3061	if (fwd > strike)
3062	expected = referenceValues[idx];
3063	else
3064	expected = referenceValues[idx] - (fwd - strike)*df;
3065
3066	const Real diff = std::fabs(x: calculated - expected);
3067	if (diff > 1e-12) {
3068	BOOST_ERROR("failed to reproduce cached extreme "
3069	"Heston model prices with exponential fitted "
3070	"Gauss-Laguerre quadrature rule"
3071	<< "\n forward : " << fwd
3072	<< "\n strike : " << strike
3073	<< "\n expected : " << expected
3074	<< "\n calculated: " << calculated
3075	<< "\n diff : " << diff
3076	<< "\n tolerance : " << 1e-12);
3077	}
3078	}
3079	}
3080	}
3081	}
3082
3083	namespace {
3084	class HestonIntegrationMaxBoundTestFct {
3085	public:
3086	explicit HestonIntegrationMaxBoundTestFct(Real maxBound)
3087	: maxBound_(maxBound),
3088	callCounter_(ext::make_shared<Size>(args: Size(0))) {}
3089
3090	Real operator()() {
3091	++(*callCounter_);
3092	return maxBound_;
3093	}
3094
3095	Size getCallCounter() const {
3096	return *callCounter_;
3097	}
3098	private:
3099	const Real maxBound_;
3100	const ext::shared_ptr<Size> callCounter_;
3101	};
3102	}
3103
3104	void HestonModelTest::testHestonEngineIntegration() {
3105	BOOST_TEST_MESSAGE("Testing Heston engine integration signature...");
3106
3107	auto square = [](Real x) -> Real { return x * x; };
3108
3109	const AnalyticHestonEngine::Integration integration =
3110	AnalyticHestonEngine::Integration::gaussLobatto(relTolerance: 1e-12, absTolerance: 1e-12);
3111
3112	const Real c1 = integration.calculate(c_inf: 1.0, f: square, maxBound: Real(1.0));
3113
3114	HestonIntegrationMaxBoundTestFct testFct(1.0);
3115	const Real c2 = integration.calculate(c_inf: 1.0, f: square, maxBound: testFct);
3116
3117	if (testFct.getCallCounter() == 0 ||
3118	std::fabs(x: c1 - 1/3.) > 1e-10 || std::fabs(x: c2 - 1/3.) > 1e-10) {
3119	BOOST_ERROR("failed to test Heston engine integration signature");
3120	}
3121	}
3122
3123
3124	void HestonModelTest::testOptimalControlVariateChoice() {
3125	BOOST_TEST_MESSAGE(
3126	"Testing optimal control variate choice for the Heston model...");
3127
3128	Real v0 = 0.0225;
3129	Real rho = 0.5;
3130	Real sigma = 2.0;
3131	Real kappa = 0.1;
3132	Real theta = 0.01;
3133	Time t = 2.0;
3134
3135	AnalyticHestonEngine::ComplexLogFormula calculated =
3136	AnalyticHestonEngine::optimalControlVariate(
3137	t, v0, kappa, theta, sigma, rho);
3138
3139	if (calculated != AnalyticHestonEngine::AsymptoticChF) {
3140	BOOST_ERROR("failed to reproduce optimal control variate choice");
3141	}
3142
3143	calculated = AnalyticHestonEngine::optimalControlVariate(
3144	t, v0, kappa, theta, sigma: 0.2, rho);
3145	if (calculated != AnalyticHestonEngine::AndersenPiterbargOptCV) {
3146	BOOST_ERROR("failed to reproduce optimal control variate choice");
3147	}
3148
3149	calculated = AnalyticHestonEngine::optimalControlVariate(
3150	t, v0: 0.2, kappa, theta, sigma, rho);
3151	if (calculated != AnalyticHestonEngine::AndersenPiterbargOptCV) {
3152	BOOST_ERROR("failed to reproduce optimal control variate choice");
3153	}
3154
3155	}
3156
3157	void HestonModelTest::testAsymptoticControlVariate() {
3158	BOOST_TEST_MESSAGE("Testing Heston asymptotic control variate...");
3159
3160	const Date todaysDate = Date(4, August, 2020);
3161	Settings::instance().evaluationDate() = todaysDate;
3162
3163	const DayCounter dc = Actual365Fixed();
3164
3165	const Handle<YieldTermStructure> rTS(flatRate(forward: 0.0, dc));
3166
3167	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 1.0));
3168
3169	const Real v0 = 0.0225;
3170	const Real rho = 0.5;
3171	const Real sigma = 2.0;
3172	const Real kappa = 0.1;
3173	const Real theta = 0.01;
3174
3175	const ext::shared_ptr<HestonModel> model =
3176	ext::make_shared<HestonModel>(
3177	args: ext::make_shared<HestonProcess>(
3178	args: rTS, args: rTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho));
3179
3180	const Date maturityDate = todaysDate + Period(2, Years);
3181	const Time t = dc.yearFraction(d1: todaysDate, d2: maturityDate);
3182	const ext::shared_ptr<Exercise> exercise =
3183	ext::make_shared<EuropeanExercise>(args: maturityDate);
3184
3185	const Real moneynesses[] = { -15, -10, -5, 0, 5, 10, 15 };
3186
3187	const Real expected[] = {
3188	0.0074676425640918,
3189	0.008680823863233695,
3190	0.010479611906112223,
3191	0.023590088942038945,
3192	0.0019575784806211706,
3193	0.0005490310253748906,
3194	0.0001657118753134695
3195	};
3196
3197	const ext::shared_ptr<PricingEngine> engines[] = {
3198	ext::make_shared<AnalyticHestonEngine>(
3199	args: model,
3200	args: AnalyticHestonEngine::OptimalCV,
3201	args: AnalyticHestonEngine::Integration::gaussLobatto(relTolerance: 1e-10, absTolerance: 1e-10, maxEvaluations: 100000)),
3202	ext::make_shared<AnalyticHestonEngine>(
3203	args: model,
3204	args: AnalyticHestonEngine::OptimalCV,
3205	args: AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder: 96)),
3206	ext::make_shared<ExponentialFittingHestonEngine>(args: model)
3207	};
3208
3209	for (Size j=0; j < LENGTH(engines); ++j) {
3210	for (Size i=0; i < LENGTH(moneynesses); ++i) {
3211	const Real moneyness = moneynesses[i];
3212
3213	const Real strike = std::exp(x: -moneyness*std::sqrt(x: theta*t));
3214
3215	const ext::shared_ptr<PlainVanillaPayoff> payoff =
3216	ext::make_shared<PlainVanillaPayoff>(
3217	args: strike > 1.0 ? Option::Call : Option::Put, args: strike);
3218
3219	const ext::shared_ptr<PricingEngine> engine = engines[j];
3220
3221	VanillaOption option(payoff, exercise);
3222	option.setPricingEngine(engine);
3223
3224	const Real calculated = option.NPV();
3225
3226	const ext::shared_ptr<AnalyticHestonEngine> analyticHestonEngine =
3227	ext::dynamic_pointer_cast<AnalyticHestonEngine>(r: engine);
3228
3229	if ((analyticHestonEngine != nullptr) &&
3230	analyticHestonEngine->numberOfEvaluations() > 5000) {
3231	BOOST_ERROR("too many function valuation needed "
3232	<< "\n moneyness : " << moneyness
3233	<< "\n evaluations : "
3234	<< analyticHestonEngine->numberOfEvaluations()
3235	<< "\n max evaluations: " << 5000);
3236	}
3237
3238	const Real diff = std::fabs(x: calculated - expected[i]);
3239	if (diff > 5e-8) {
3240	BOOST_ERROR("failed to reproduce extreme Heston model values for"
3241	<< "\n moneyness : " << moneyness
3242	<< "\n #engine : " << j
3243	<< "\n calculated: " << calculated
3244	<< "\n expected : " << expected[i]
3245	<< "\n difference: " << diff
3246	<< "\n tolerance : " << 1e-8);
3247	}
3248	}
3249	}
3250	}
3251
3252	void HestonModelTest::testLocalVolFromHestonModel() {
3253	BOOST_TEST_MESSAGE("Testing Local Volatility pricing from Heston Model...");
3254
3255	const auto todaysDate = Date(28, June, 2021);
3256	Settings::instance().evaluationDate() = todaysDate;
3257
3258	const auto dc = Actual365Fixed();
3259
3260	const Handle<YieldTermStructure> rTS(
3261	ext::make_shared<ZeroCurve>(
3262	args: std::vector<Date>{
3263	todaysDate, todaysDate + Period(90, Days),
3264	todaysDate + Period(180, Days), todaysDate + Period(1, Years)
3265	},
3266	args: std::vector<Rate>{0.075, 0.05, 0.075, 0.1},
3267	args: dc
3268	)
3269	);
3270
3271	const Handle<YieldTermStructure> qTS(
3272	ext::make_shared<ZeroCurve>(
3273	args: std::vector<Date>{
3274	todaysDate, todaysDate + Period(90, Days), todaysDate + Period(1, Years)
3275	},
3276	args: std::vector<Rate>{0.06, 0.04, 0.12},
3277	args: dc
3278	)
3279	);
3280
3281	const Handle<Quote> s0(ext::make_shared<SimpleQuote>(args: 100.0));
3282
3283	const Real v0 = 0.1;
3284	const Real rho = -0.75;
3285	const Real sigma = 0.8;
3286	const Real kappa = 1.0;
3287	const Real theta = 0.16;
3288
3289	const auto hestonModel = ext::make_shared<HestonModel>(
3290	args: ext::make_shared<HestonProcess>(
3291	args: rTS, args: qTS, args: s0, args: v0, args: kappa, args: theta, args: sigma, args: rho)
3292	);
3293
3294	VanillaOption option(
3295	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: 120.0),
3296	ext::make_shared<EuropeanExercise>(args: todaysDate + Period(1, Years))
3297	);
3298
3299	option.setPricingEngine(
3300	ext::make_shared<AnalyticHestonEngine>(
3301	args: hestonModel,
3302	args: AnalyticHestonEngine::OptimalCV,
3303	args: AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder: 192)
3304	)
3305	);
3306
3307	const Real expected = option.NPV();
3308
3309	option.setPricingEngine(
3310	ext::make_shared<FdBlackScholesVanillaEngine>(
3311	args: ext::make_shared<BlackScholesMertonProcess>(
3312	args: s0, args: qTS, args: rTS,
3313	args: Handle<BlackVolTermStructure>(
3314	ext::make_shared<HestonBlackVolSurface>(
3315	args: Handle<HestonModel>(hestonModel),
3316	args: AnalyticHestonEngine::OptimalCV,
3317	args: AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder: 24)
3318	)
3319	)
3320	),
3321	args: 25, args: 125, args: 1, args: FdmSchemeDesc::Douglas(), args: true, args: 0.4
3322	)
3323	);
3324
3325	const Real calculated = option.NPV();
3326
3327	const Real tol = 0.002;
3328	const Real diff = std::fabs(x: calculated - expected);
3329	if (diff > tol) {
3330	BOOST_ERROR("failed to reproduce Heston model values with "
3331	"local volatility pricing"
3332	<< "\n calculated: " << calculated
3333	<< "\n expected : " << expected
3334	<< "\n difference: " << diff
3335	<< "\n tolerance : " << tol);
3336	}
3337	}
3338
3339
3340	test_suite* HestonModelTest::suite(SpeedLevel speed) {
3341	auto* suite = BOOST_TEST_SUITE("Heston model tests");
3342
3343	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testBlackCalibration));
3344	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testDAXCalibration));
3345	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAnalyticVsBlack));
3346	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAnalyticVsCached));
3347	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testMultipleStrikesEngine));
3348	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testMcVsCached));
3349	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testFdVanillaVsCached));
3350	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testFdAmerican));
3351	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAnalyticPiecewiseTimeDependent));
3352	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testDAXCalibrationOfTimeDependentModel));
3353	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAlanLewisReferencePrices));
3354	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testExpansionOnAlanLewisReference));
3355	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testExpansionOnFordeReference));
3356	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAllIntegrationMethods));
3357	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testCosHestonCumulants));
3358	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testCosHestonEngine));
3359	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testCosHestonEngineTruncation));
3360	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testCharacteristicFct));
3361	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAndersenPiterbargPricing));
3362	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAndersenPiterbargControlVariateIntegrand));
3363	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAndersenPiterbargConvergence));
3364	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testPiecewiseTimeDependentChFvsHestonChF));
3365	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testPiecewiseTimeDependentComparison));
3366	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testPiecewiseTimeDependentChFAsymtotic));
3367	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testSmallSigmaExpansion));
3368	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testSmallSigmaExpansion4ExpFitting));
3369	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testExponentialFitting4StrikesAndMaturities));
3370	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testHestonEngineIntegration));
3371	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testOptimalControlVariateChoice));
3372	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAsymptoticControlVariate));
3373	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testLocalVolFromHestonModel));
3374
3375	if (speed <= Fast) {
3376	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testDifferentIntegrals));
3377	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testFdVanillaWithDividendsVsCached));
3378	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testFdBarrierVsCached));
3379	}
3380
3381	if (speed == Slow) {
3382	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testKahlJaeckelCase));
3383	}
3384
3385	return suite;
3386	}
3387
3388	test_suite* HestonModelTest::experimental() {
3389	auto* suite = BOOST_TEST_SUITE("Heston model experimental tests");
3390	suite->add(QUANTLIB_TEST_CASE(&HestonModelTest::testAnalyticPDFHestonEngine));
3391	return suite;
3392	}
3393