1	/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3	/*
4	Copyright (C) 2014, 2015 Johannes Göttker-Schnetmann
5	Copyright (C) 2014, 2015 Klaus Spanderen
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 "hestonslvmodel.hpp"
22	#include "utilities.hpp"
23	#include <ql/quotes/simplequote.hpp>
24	#include <ql/time/calendars/target.hpp>
25	#include <ql/time/daycounters/actualactual.hpp>
26	#include <ql/math/functional.hpp>
27	#include <ql/math/solvers1d/brent.hpp>
28	#include <ql/instruments/barrieroption.hpp>
29	#include <ql/instruments/forwardvanillaoption.hpp>
30	#include <ql/instruments/impliedvolatility.hpp>
31	#include <ql/math/distributions/gammadistribution.hpp>
32	#include <ql/math/distributions/normaldistribution.hpp>
33	#include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
34	#include <ql/math/interpolations/bilinearinterpolation.hpp>
35	#include <ql/math/integrals/gausslobattointegral.hpp>
36	#include <ql/math/integrals/discreteintegrals.hpp>
37	#include <ql/math/randomnumbers/rngtraits.hpp>
38	#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
39	#include <ql/math/optimization/levenbergmarquardt.hpp>
40	#include <ql/models/equity/hestonmodel.hpp>
41	#include <ql/models/equity/hestonmodelhelper.hpp>
42	#include <ql/processes/hestonprocess.hpp>
43	#include <ql/termstructures/yield/zerocurve.hpp>
44	#include <ql/termstructures/volatility/equityfx/blackvariancesurface.hpp>
45	#include <ql/termstructures/volatility/equityfx/noexceptlocalvolsurface.hpp>
46	#include <ql/termstructures/volatility/equityfx/fixedlocalvolsurface.hpp>
47	#include <ql/termstructures/volatility/equityfx/gridmodellocalvolsurface.hpp>
48	#include <ql/termstructures/volatility/equityfx/localconstantvol.hpp>
49	#include <ql/termstructures/volatility/equityfx/localvolsurface.hpp>
50	#include <ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp>
51	#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
52	#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
53	#include <ql/pricingengines/vanilla/fdhestonvanillaengine.hpp>
54	#include <ql/pricingengines/barrier/fdhestonbarrierengine.hpp>
55	#include <ql/pricingengines/barrier/fdblackscholesbarrierengine.hpp>
56	#include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>
57	#include <ql/pricingengines/vanilla/mceuropeanhestonengine.hpp>
58	#include <ql/pricingengines/forward/forwardengine.hpp>
59	#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
60	#include <ql/methods/finitedifferences/meshers/fdmmeshercomposite.hpp>
61	#include <ql/methods/finitedifferences/meshers/fdmblackscholesmesher.hpp>
62	#include <ql/methods/finitedifferences/meshers/predefined1dmesher.hpp>
63	#include <ql/methods/finitedifferences/meshers/uniform1dmesher.hpp>
64	#include <ql/methods/finitedifferences/meshers/concentrating1dmesher.hpp>
65	#include <ql/methods/finitedifferences/schemes/douglasscheme.hpp>
66	#include <ql/methods/finitedifferences/schemes/hundsdorferscheme.hpp>
67	#include <ql/methods/finitedifferences/schemes/craigsneydscheme.hpp>
68	#include <ql/methods/finitedifferences/schemes/modifiedcraigsneydscheme.hpp>
69	#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
70	#include <ql/methods/finitedifferences/utilities/fdmmesherintegral.hpp>
71	#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
72	#include <ql/methods/finitedifferences/operators/fdmlocalvolfwdop.hpp>
73	#include <ql/models/marketmodels/browniangenerators/mtbrowniangenerator.hpp>
74	#include <ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp>
75	#include <ql/models/equity/hestonslvfdmmodel.hpp>
76	#include <ql/models/equity/hestonslvmcmodel.hpp>
77	#include <ql/methods/finitedifferences/operators/fdmhestonfwdop.hpp>
78	#include <ql/methods/finitedifferences/operators/fdmsquarerootfwdop.hpp>
79	#include <ql/methods/finitedifferences/operators/fdmblackscholesfwdop.hpp>
80	#include <ql/methods/finitedifferences/utilities/fdmhestongreensfct.hpp>
81	#include <ql/methods/finitedifferences/utilities/localvolrndcalculator.hpp>
82	#include <ql/methods/finitedifferences/utilities/squarerootprocessrndcalculator.hpp>
83	#include <ql/pricingengines/barrier/fdhestondoublebarrierengine.hpp>
84	#include <ql/experimental/exoticoptions/analyticpdfhestonengine.hpp>
85	#include <ql/processes/hestonslvprocess.hpp>
86	#include <ql/instruments/doublebarrieroption.hpp>
87	#include <ql/pricingengines/barrier/analyticdoublebarrierbinaryengine.hpp>
88	#include <boost/math/special_functions/gamma.hpp>
89	#include <boost/multi_array.hpp>
90	#include <iomanip>
91
92	using namespace QuantLib;
93	using boost::unit_test_framework::test_suite;
94
95	namespace {
96	Real fokkerPlanckPrice1D(const ext::shared_ptr<FdmMesher>& mesher,
97	const ext::shared_ptr<FdmLinearOpComposite>& op,
98	const ext::shared_ptr<StrikedTypePayoff>& payoff,
99	Real x0, Time maturity, Size tGrid) {
100
101	const Array x = mesher->locations(direction: 0);
102	Array p(x.size(), 0.0);
103
104	QL_REQUIRE(x.size() > 3 && x[1] <= x0 && x[x.size()-2] >= x0,
105	"insufficient mesher");
106
107	const Array::const_iterator upperb
108	= std::upper_bound(first: x.begin(), last: x.end(), val: x0);
109	const Array::const_iterator lowerb = upperb-1;
110
111	if (close_enough(x: *upperb, y: x0)) {
112	const Size idx = std::distance(first: x.begin(), last: upperb);
113	const Real dx = (x[idx+1]-x[idx-1])/2.0;
114	p[idx] = 1.0/dx;
115	}
116	else if (close_enough(x: *lowerb, y: x0)) {
117	const Size idx = std::distance(first: x.begin(), last: lowerb);
118	const Real dx = (x[idx+1]-x[idx-1])/2.0;
119	p[idx] = 1.0/dx;
120	} else {
121	const Real dx = *upperb - *lowerb;
122	const Real lowerP = (*upperb - x0)/dx;
123	const Real upperP = (x0 - *lowerb)/dx;
124
125	const Size lowerIdx = std::distance(first: x.begin(), last: lowerb);
126	const Size upperIdx = std::distance(first: x.begin(), last: upperb);
127
128	const Real lowerDx = (x[lowerIdx+1]-x[lowerIdx-1])/2.0;
129	const Real upperDx = (x[upperIdx+1]-x[upperIdx-1])/2.0;
130
131	p[lowerIdx] = lowerP/lowerDx;
132	p[upperIdx] = upperP/upperDx;
133	}
134
135	DouglasScheme evolver(FdmSchemeDesc::Douglas().theta, op);
136	const Time dt = maturity/tGrid;
137	evolver.setStep(dt);
138
139	for (Time t=dt; t <= maturity+20*QL_EPSILON; t+=dt) {
140	evolver.step(a&: p, t);
141	}
142
143	Array payoffTimesDensity(x.size());
144	for (Size i=0; i < x.size(); ++i) {
145	payoffTimesDensity[i] = (*payoff)(std::exp(x: x[i]))*p[i];
146	}
147
148	CubicNaturalSpline f(x.begin(), x.end(), payoffTimesDensity.begin());
149	f.enableExtrapolation();
150	return GaussLobattoIntegral(1000, 1e-6)(f, x.front(), x.back());
151	}
152	}
153
154	void HestonSLVModelTest::testBlackScholesFokkerPlanckFwdEquation() {
155	BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation for BS process...");
156
157	const DayCounter dc = ActualActual(ActualActual::ISDA);
158	const Date todaysDate = Date(28, Dec, 2012);
159	Settings::instance().evaluationDate() = todaysDate;
160
161	const Date maturityDate = todaysDate + Period(2, Years);
162	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
163
164	const Real s0 = 100;
165	const Real x0 = std::log(x: s0);
166	const Rate r = 0.035;
167	const Rate q = 0.01;
168	const Volatility v = 0.35;
169
170	const Size xGrid = 2*100+1;
171	const Size tGrid = 400;
172
173	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
174	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
175	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
176	const Handle<BlackVolTermStructure> vTS(flatVol(volatility: v, dc));
177
178	const ext::shared_ptr<GeneralizedBlackScholesProcess> process(
179	ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS, args: vTS));
180
181	const ext::shared_ptr<PricingEngine> engine(
182	ext::make_shared<AnalyticEuropeanEngine>(args: process));
183
184	const ext::shared_ptr<FdmMesher> uniformMesher(
185	ext::make_shared<FdmMesherComposite>(
186	args: ext::make_shared<FdmBlackScholesMesher>(args: xGrid, args: process, args: maturity, args: s0)));
187
188	const ext::shared_ptr<FdmLinearOpComposite> uniformBSFwdOp(
189	ext::make_shared<FdmBlackScholesFwdOp>(args: uniformMesher, args: process, args: s0, args: false));
190
191	const ext::shared_ptr<FdmMesher> concentratedMesher(
192	ext::make_shared<FdmMesherComposite>(
193	args: ext::make_shared<FdmBlackScholesMesher>(args: xGrid, args: process, args: maturity, args: s0,
194	args: Null<Real>(), args: Null<Real>(), args: 0.0001, args: 1.5,
195	args: std::pair<Real, Real>(s0, 0.1))));
196
197	const ext::shared_ptr<FdmLinearOpComposite> concentratedBSFwdOp(
198	ext::make_shared<FdmBlackScholesFwdOp>(args: concentratedMesher, args: process, args: s0, args: false));
199
200	const ext::shared_ptr<FdmMesher> shiftedMesher(
201	ext::make_shared<FdmMesherComposite>(
202	args: ext::make_shared<FdmBlackScholesMesher>(args: xGrid, args: process, args: maturity, args: s0,
203	args: Null<Real>(), args: Null<Real>(), args: 0.0001, args: 1.5,
204	args: std::pair<Real, Real>(s0*1.1, 0.2))));
205
206	const ext::shared_ptr<FdmLinearOpComposite> shiftedBSFwdOp(
207	ext::make_shared<FdmBlackScholesFwdOp>(args: shiftedMesher, args: process, args: s0, args: false));
208
209	const ext::shared_ptr<Exercise> exercise(
210	ext::make_shared<EuropeanExercise>(args: maturityDate));
211	const Real strikes[] = { 50, 80, 100, 130, 150 };
212
213	for (Real strike : strikes) {
214	const ext::shared_ptr<StrikedTypePayoff> payoff(
215	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args&: strike));
216
217	VanillaOption option(payoff, exercise);
218	option.setPricingEngine(engine);
219
220	const Real expected = option.NPV()/rTS->discount(d: maturityDate);
221	const Real calcUniform
222	= fokkerPlanckPrice1D(mesher: uniformMesher, op: uniformBSFwdOp,
223	payoff, x0, maturity, tGrid);
224	const Real calcConcentrated
225	= fokkerPlanckPrice1D(mesher: concentratedMesher, op: concentratedBSFwdOp,
226	payoff, x0, maturity, tGrid);
227	const Real calcShifted
228	= fokkerPlanckPrice1D(mesher: shiftedMesher, op: shiftedBSFwdOp,
229	payoff, x0, maturity, tGrid);
230	const Real tol = 0.02;
231
232	if (std::fabs(x: expected - calcUniform) > tol) {
233	BOOST_FAIL("failed to reproduce european option price "
234	<< "with an uniform mesher"
235	<< "\n strike: " << strike << std::fixed << std::setprecision(8)
236	<< "\n calculated: " << calcUniform << "\n expected: " << expected
237	<< "\n tolerance: " << tol);
238	}
239	if (std::fabs(x: expected - calcConcentrated) > tol) {
240	BOOST_FAIL("failed to reproduce european option price "
241	<< "with a concentrated mesher"
242	<< "\n strike: " << strike << std::fixed << std::setprecision(8)
243	<< "\n calculated: " << calcConcentrated << "\n expected: " << expected
244	<< "\n tolerance: " << tol);
245	}
246	if (std::fabs(x: expected - calcShifted) > tol) {
247	BOOST_FAIL("failed to reproduce european option price "
248	<< "with a shifted mesher"
249	<< "\n strike: " << strike << std::fixed << std::setprecision(8)
250	<< "\n calculated: " << calcShifted << "\n expected: " << expected
251	<< "\n tolerance: " << tol);
252	}
253	}
254	}
255
256
257	namespace {
258	Real stationaryLogProbabilityFct(Real kappa, Real theta,
259	Real sigma, Real z) {
260	const Real alpha = 2*kappa*theta/(sigma*sigma);
261	const Real beta = alpha/theta;
262
263	return std::pow(x: beta, y: alpha)*std::exp(x: z*alpha)
264	*std::exp(x: -beta*std::exp(x: z)-GammaFunction().logValue(x: alpha));
265	}
266	}
267
268	void HestonSLVModelTest::testSquareRootZeroFlowBC() {
269	BOOST_TEST_MESSAGE("Testing zero-flow BC for the square root process...");
270
271	const Real kappa = 1.0;
272	const Real theta = 0.4;
273	const Real sigma = 0.8;
274	const Real v_0 = 0.1;
275	const Time t = 1.0;
276
277	const Real vmin = 0.0005;
278	const Real h = 0.0001;
279
280	const Real expected[5][5]
281	= {{ 0.000548, -0.000245, -0.005657, -0.001167, -0.000024},
282	{-0.000595, -0.000701, -0.003296, -0.000883, -0.000691},
283	{-0.001277, -0.001320, -0.003128, -0.001399, -0.001318},
284	{-0.001979, -0.002002, -0.003425, -0.002047, -0.002001},
285	{-0.002715, -0.002730, -0.003920, -0.002760, -0.002730} };
286
287	for (Size i=0; i < 5; ++i) {
288	const Real v = vmin + i*0.001;
289	const Real vm2 = v - 2*h;
290	const Real vm1 = v - h;
291	const Real v0 = v;
292	const Real v1 = v + h;
293	const Real v2 = v + 2*h;
294
295	const SquareRootProcessRNDCalculator rndCalculator(
296	v_0, kappa, theta, sigma);
297
298	const Real pm2 = rndCalculator.pdf(v: vm2, t);
299	const Real pm1 = rndCalculator.pdf(v: vm1, t);
300	const Real p0 = rndCalculator.pdf(v: v0 , t);
301	const Real p1 = rndCalculator.pdf(v: v1 , t);
302	const Real p2 = rndCalculator.pdf(v: v2 , t);
303
304	// test derivatives
305	const Real flowSym2Order = sigma*sigma*v0/(4*h)*(p1-pm1)
306	+ (kappa*(v0-theta)+sigma*sigma/2)*p0;
307
308	const Real flowSym4Order
309	= sigma*sigma*v0/(24*h)*(-p2 + 8*p1 - 8*pm1 + pm2)
310	+ (kappa*(v0-theta)+sigma*sigma/2)*p0;
311
312	const Real fwd1Order = sigma*sigma*v0/(2*h)*(p1-p0)
313	+ (kappa*(v0-theta)+sigma*sigma/2)*p0;
314
315	const Real fwd2Order = sigma*sigma*v0/(4*h)*(4*p1-3*p0-p2)
316	+ (kappa*(v0-theta)+sigma*sigma/2)*p0;
317
318	const Real fwd3Order
319	= sigma*sigma*v0/(12*h)*(-p2 + 6*p1 - 3*p0 - 2*pm1)
320	+ (kappa*(v0-theta)+sigma*sigma/2)*p0;
321
322	const Real tol = 0.000002;
323	if ( std::fabs(x: expected[i][0] - flowSym2Order) > tol
324	|| std::fabs(x: expected[i][1] - flowSym4Order) > tol
325	|| std::fabs(x: expected[i][2] - fwd1Order) > tol
326	|| std::fabs(x: expected[i][3] - fwd2Order) > tol
327	|| std::fabs(x: expected[i][4] - fwd3Order) > tol ) {
328	BOOST_ERROR("failed to reproduce Zero Flow BC at"
329	<< "\n v: " << v
330	<< "\n tolerance: " << tol);
331	}
332	}
333	}
334
335
336	namespace {
337	ext::shared_ptr<FdmMesher> createStationaryDistributionMesher(
338	Real kappa, Real theta, Real sigma, Size vGrid) {
339
340	const Real qMin = 0.01;
341	const Real qMax = 0.99;
342	const Real dq = (qMax-qMin)/(vGrid-1);
343
344	const SquareRootProcessRNDCalculator rnd(theta, kappa, theta, sigma);
345	std::vector<Real> v(vGrid);
346	for (Size i=0; i < vGrid; ++i) {
347	v[i] = rnd.stationary_invcdf(q: qMin + i*dq);
348	}
349
350	return ext::shared_ptr<FdmMesher>(
351	ext::make_shared<FdmMesherComposite>(
352	args: ext::make_shared<Predefined1dMesher>(args&: v)));
353	}
354	}
355
356
357	void HestonSLVModelTest::testTransformedZeroFlowBC() {
358	BOOST_TEST_MESSAGE("Testing zero-flow BC for transformed "
359	"Fokker-Planck forward equation...");
360
361	const Real kappa = 1.0;
362	const Real theta = 0.4;
363	const Real sigma = 2.0;
364	const Size vGrid = 100;
365
366	const ext::shared_ptr<FdmMesher> mesher
367	= createStationaryDistributionMesher(kappa, theta, sigma, vGrid);
368	const Array v = mesher->locations(direction: 0);
369
370	Array p(vGrid);
371	const SquareRootProcessRNDCalculator rnd(theta, kappa, theta, sigma);
372	for (Size i=0; i < v.size(); ++i)
373	p[i] = rnd.stationary_pdf(v: v[i]);
374
375
376	const Real alpha = 1.0 - 2*kappa*theta/(sigma*sigma);
377	const Array q = Pow(v, alpha)*p;
378
379	for (Size i=0; i < vGrid/2; ++i) {
380	const Real hm = v[i+1] - v[i];
381	const Real hp = v[i+2] - v[i+1];
382
383	const Real eta=1.0/(hm*(hm+hp)*hp);
384	const Real a = -eta*(squared(x: hm+hp) - hm*hm);
385	const Real b = eta*squared(x: hm+hp);
386	const Real c = -eta*hm*hm;
387
388	const Real df = a*q[i] + b*q[i+1] + c*q[i+2];
389	const Real flow = 0.5*sigma*sigma*v[i]*df + kappa*v[i]*q[i];
390
391	const Real tol = 1e-6;
392	if (std::fabs(x: flow) > tol) {
393	BOOST_ERROR("failed to reproduce Zero Flow BC at"
394	<< "\n v: " << v
395	<< "\n flow: " << flow
396	<< "\n tolerance: " << tol);
397	}
398	}
399	}
400
401	namespace {
402	class q_fct {
403	public:
404	q_fct(const Array& v, const Array& p, const Real alpha)
405	: v_(v), q_(Pow(v, alpha)*p), alpha_(alpha),
406	spline_(ext::make_shared<CubicNaturalSpline>(args: v_.begin(), args: v_.end(), args: q_.begin())) {
407	}
408
409	Real operator()(Real v) const {
410	return (*spline_)(v, true)*std::pow(x: v, y: -alpha_);
411	}
412	private:
413	const Array v_, q_;
414	const Real alpha_;
415	const ext::shared_ptr<CubicNaturalSpline> spline_;
416	};
417	}
418
419	void HestonSLVModelTest::testSquareRootEvolveWithStationaryDensity() {
420	BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation "
421	"for the square root process with stationary density...");
422
423	// Documentation for this test case:
424	// http://www.spanderen.de/2013/05/04/fokker-planck-equation-feller-constraint-and-boundary-conditions/
425	const Real kappa = 2.5;
426	const Real theta = 0.2;
427	const Size vGrid = 100;
428	const Real eps = 1e-2;
429
430	for (Real sigma = 0.2; sigma < 2.01; sigma+=0.1) {
431	const Real alpha = (1.0 - 2*kappa*theta/(sigma*sigma));
432
433	const SquareRootProcessRNDCalculator rnd(theta, kappa, theta, sigma);
434	const Real vMin = rnd.stationary_invcdf(q: eps);
435	const Real vMax = rnd.stationary_invcdf(q: 1-eps);
436
437	const ext::shared_ptr<FdmMesher> mesher(
438	ext::make_shared<FdmMesherComposite>(
439	args: ext::make_shared<Uniform1dMesher>(args: vMin, args: vMax, args: vGrid)));
440
441	const Array v = mesher->locations(direction: 0);
442	const FdmSquareRootFwdOp::TransformationType transform =
443	(sigma < 0.75) ?
444	FdmSquareRootFwdOp::Plain :
445	FdmSquareRootFwdOp::Power;
446
447	Array vq (v.size());
448	Array vmq(v.size());
449	for (Size i=0; i < v.size(); ++i) {
450	vmq[i] = 1.0/(vq[i] = std::pow(x: v[i], y: alpha));
451	}
452
453	Array p(vGrid);
454	for (Size i=0; i < v.size(); ++i) {
455	p[i] = rnd.stationary_pdf(v: v[i]);
456	if (transform == FdmSquareRootFwdOp::Power)
457	p[i] *= vq[i];
458	}
459
460	const ext::shared_ptr<FdmSquareRootFwdOp> op(
461	ext::make_shared<FdmSquareRootFwdOp>(args: mesher, args: kappa, args: theta,
462	args&: sigma, args: 0, args: transform));
463
464	const Size n = 100;
465	const Time dt = 0.01;
466	DouglasScheme evolver(0.5, op);
467	evolver.setStep(dt);
468
469	for (Size i=1; i <= n; ++i) {
470	evolver.step(a&: p, t: i*dt);
471	}
472
473	const Real expected = 1-2*eps;
474
475	if (transform == FdmSquareRootFwdOp::Power)
476	for (Size i=0; i < v.size(); ++i) {
477	p[i] *= vmq[i];
478	}
479
480	const q_fct f(v, p, alpha);
481
482	const Real calculated = GaussLobattoIntegral(1000000, 1e-6)(
483	f, v.front(), v.back());
484
485	const Real tol = 0.005;
486	if (std::fabs(x: calculated-expected) > tol) {
487	BOOST_ERROR("failed to reproduce stationary probability function"
488	<< "\n calculated: " << calculated
489	<< "\n expected: " << expected
490	<< "\n tolerance: " << tol);
491	}
492	}
493	}
494
495	void HestonSLVModelTest::testSquareRootLogEvolveWithStationaryDensity() {
496	BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation "
497	"for the square root log process with stationary density...");
498
499	// Documentation for this test case:
500	// nowhere yet :)
501	const Real kappa = 2.5;
502	const Real theta = 0.2;
503	const Size vGrid = 1000;
504	const Real eps = 1e-2;
505
506	for (Real sigma = 0.2; sigma < 2.01; sigma+=0.1) {
507	const Real lowerLimit = 0.001;
508	// should not go to very large negative values, distributions flattens with sigma
509	// causing numerical instabilities log/exp evaluations
510
511	const SquareRootProcessRNDCalculator rnd(theta, kappa, theta, sigma);
512
513	const Real vMin = std::max(a: lowerLimit, b: rnd.stationary_invcdf(q: eps));
514	const Real lowEps = std::max(a: eps, b: rnd.stationary_cdf(v: lowerLimit));
515
516	const Real expected = 1-eps-lowEps;
517	const Real vMax = rnd.stationary_invcdf(q: 1-eps);
518
519	const ext::shared_ptr<FdmMesherComposite> mesher(
520	ext::make_shared<FdmMesherComposite>(
521	args: ext::make_shared<Uniform1dMesher>(args: std::log(x: vMin), args: std::log(x: vMax), args: vGrid)));
522
523	const Array v = mesher->locations(direction: 0);
524
525	Array p(vGrid);
526	for (Size i=0; i < v.size(); ++i)
527	p[i] = stationaryLogProbabilityFct(kappa, theta, sigma, z: v[i]);
528
529	const ext::shared_ptr<FdmSquareRootFwdOp> op(
530	ext::make_shared<FdmSquareRootFwdOp>(args: mesher, args: kappa, args: theta,
531	args&: sigma, args: 0,
532	args: FdmSquareRootFwdOp::Log));
533
534	const Size n = 100;
535	const Time dt = 0.01;
536	FdmBoundaryConditionSet bcSet;
537	DouglasScheme evolver(0.5, op);
538	evolver.setStep(dt);
539
540	for (Size i=1; i <= n; ++i) {
541	evolver.step(a&: p, t: i*dt);
542	}
543
544	const Real calculated
545	= FdmMesherIntegral(mesher, DiscreteSimpsonIntegral()).integrate(f: p);
546
547	const Real tol = 0.005;
548	if (std::fabs(x: calculated-expected) > tol) {
549	BOOST_ERROR("failed to reproduce stationary probability function for "
550	<< "\n sigma: " << sigma
551	<< "\n calculated: " << calculated
552	<< "\n expected: " << expected
553	<< "\n tolerance: " << tol);
554	}
555	}
556	}
557
558	void HestonSLVModelTest::testSquareRootFokkerPlanckFwdEquation() {
559	BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation "
560	"for the square root process with Dirac start...");
561
562	const Real kappa = 1.2;
563	const Real theta = 0.4;
564	const Real sigma = 0.7;
565	const Real v0 = theta;
566	const Real alpha = 1.0 - 2*kappa*theta/(sigma*sigma);
567
568	const Time maturity = 1.0;
569
570	const Size xGrid = 1001;
571	const Size tGrid = 500;
572
573	const Real vol = sigma*std::sqrt(x: theta/(2*kappa));
574	const Real upperBound = theta+6*vol;
575	const Real lowerBound = std::max(a: 0.0002, b: theta-6*vol);
576
577	const ext::shared_ptr<FdmMesher> mesher(
578	ext::make_shared<FdmMesherComposite>(
579	args: ext::make_shared<Uniform1dMesher>(args: lowerBound, args: upperBound, args: xGrid)));
580
581	const Array x(mesher->locations(direction: 0));
582
583	const ext::shared_ptr<FdmSquareRootFwdOp> op(
584	ext::make_shared<FdmSquareRootFwdOp>(args: mesher, args: kappa, args: theta, args: sigma, args: 0)); //!
585
586	const Time dt = maturity/tGrid;
587	const Size n = 5;
588
589	Array p(xGrid);
590	SquareRootProcessRNDCalculator rndCalculator(v0, kappa, theta, sigma);
591	for (Size i=0; i < p.size(); ++i) {
592	p[i] = rndCalculator.pdf(v: x[i], t: n*dt);
593	}
594	Array q = Pow(v: x, alpha)*p;
595
596	DouglasScheme evolver(0.5, op);
597	evolver.setStep(dt);
598
599	for (Time t=(n+1)*dt; t <= maturity+20*QL_EPSILON; t+=dt) {
600	evolver.step(a&: p, t);
601	evolver.step(a&: q, t);
602	}
603
604	const Real tol = 0.002;
605
606	Array y(x.size());
607	for (Size i=0; i < x.size(); ++i) {
608	const Real expected = rndCalculator.pdf(v: x[i], t: maturity);
609
610	const Real calculated = p[i];
611	if (std::fabs(x: expected - calculated) > tol) {
612	BOOST_FAIL("failed to reproduce pdf at"
613	<< std::fixed << std::setprecision(5)
614	<< "\n x: " << x[i]
615	<< "\n calculated: " << calculated
616	<< "\n expected: " << expected
617	<< "\n tolerance: " << tol);
618	}
619	}
620	}
621
622
623
624	namespace {
625	Real fokkerPlanckPrice2D(const Array& p,
626	const ext::shared_ptr<FdmMesherComposite>& mesher) {
627
628	std::vector<Real> x, y;
629
630	x.reserve(n: mesher->layout()->dim()[0]);
631	y.reserve(n: mesher->layout()->dim()[1]);
632
633	for (const auto& iter : *mesher->layout()) {
634	if (iter.coordinates()[1] == 0U) {
635	x.push_back(x: mesher->location(iter, direction: 0));
636	}
637	if (iter.coordinates()[0] == 0U) {
638	y.push_back(x: mesher->location(iter, direction: 1));
639	}
640	}
641
642	return FdmMesherIntegral(mesher,
643	DiscreteSimpsonIntegral()).integrate(f: p);
644	}
645
646	Real hestonPxBoundary(
647	Time maturity, Real eps,
648	const ext::shared_ptr<HestonModel>& model) {
649
650	const AnalyticPDFHestonEngine pdfEngine(model);
651	const Real sInit = model->process()->s0()->value();
652	const Real xMin = Brent().solve(
653	f: [&](Real x) -> Real { return pdfEngine.cdf(X: x, t: maturity) - eps; },
654	accuracy: sInit*1e-3, guess: sInit, xMin: sInit*0.001, xMax: 1000*sInit);
655
656	return xMin;
657	}
658
659	struct FokkerPlanckFwdTestCase {
660	const Real s0, r, q, v0, kappa, theta, rho, sigma;
661	const Size xGrid, vGrid, tGridPerYear, tMinGridPerYear;
662	const Real avgEps, eps;
663	const FdmSquareRootFwdOp::TransformationType trafoType;
664	const FdmHestonGreensFct::Algorithm greensAlgorithm;
665	const FdmSchemeDesc::FdmSchemeType schemeType;
666	};
667
668	void hestonFokkerPlanckFwdEquationTest(
669	const FokkerPlanckFwdTestCase& testCase) {
670
671	const DayCounter dc = ActualActual(ActualActual::ISDA);
672	const Date todaysDate = Date(28, Dec, 2014);
673	Settings::instance().evaluationDate() = todaysDate;
674
675	std::vector<Period> maturities = {
676	Period(1, Months), Period(3, Months), Period(6, Months), Period(9, Months),
677	Period(1, Years), Period(2, Years), Period(3, Years)
678	};
679
680	const Date maturityDate = todaysDate + maturities.back();
681	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
682
683	const Real s0 = testCase.s0;
684	const Real x0 = std::log(x: s0);
685	const Rate r = testCase.r;
686	const Rate q = testCase.q;
687
688	const Real kappa = testCase.kappa;
689	const Real theta = testCase.theta;
690	const Real rho = testCase.rho;
691	const Real sigma = testCase.sigma;
692	const Real v0 = testCase.v0;
693	const Real alpha = 1.0 - 2*kappa*theta/(sigma*sigma);
694
695	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
696	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
697	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
698
699	const ext::shared_ptr<HestonProcess> process(
700	ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
701
702	const ext::shared_ptr<HestonModel> model(ext::make_shared<HestonModel>(args: process));
703
704	const ext::shared_ptr<PricingEngine> engine(
705	ext::make_shared<AnalyticHestonEngine>(args: model));
706
707	const Size xGrid = testCase.xGrid;
708	const Size vGrid = testCase.vGrid;
709	const Size tGridPerYear = testCase.tGridPerYear;
710
711	const FdmSquareRootFwdOp::TransformationType transformationType
712	= testCase.trafoType;
713	Real lowerBound, upperBound;
714	std::vector<ext::tuple<Real, Real, bool> > cPoints;
715
716	const SquareRootProcessRNDCalculator rnd(v0, kappa, theta, sigma);
717	switch (transformationType) {
718	case FdmSquareRootFwdOp::Log:
719	{
720	upperBound = std::log(x: rnd.stationary_invcdf(q: 0.9995));
721	lowerBound = std::log(x: 0.00001);
722
723	const Real v0Center = std::log(x: v0);
724	const Real v0Density = 10.0;
725	const Real upperBoundDensity = 100;
726	const Real lowerBoundDensity = 1.0;
727	cPoints = {
728	ext::make_tuple(args&: lowerBound, args: lowerBoundDensity, args: false),
729	ext::make_tuple(args: v0Center, args: v0Density, args: true),
730	ext::make_tuple(args&: upperBound, args: upperBoundDensity, args: false)
731	};
732	}
733	break;
734	case FdmSquareRootFwdOp::Plain:
735	{
736	upperBound = rnd.stationary_invcdf(q: 0.9995);
737	lowerBound = rnd.stationary_invcdf(q: 1e-5);
738
739	const Real v0Center = v0;
740	const Real v0Density = 0.1;
741	const Real lowerBoundDensity = 0.0001;
742	cPoints = {
743	ext::make_tuple(args&: lowerBound, args: lowerBoundDensity, args: false),
744	ext::make_tuple(args: v0Center, args: v0Density, args: true)
745	};
746	}
747	break;
748	case FdmSquareRootFwdOp::Power:
749	{
750	upperBound = rnd.stationary_invcdf(q: 0.9995);
751	lowerBound = 0.000075;
752
753	const Real v0Center = v0;
754	const Real v0Density = 1.0;
755	const Real lowerBoundDensity = 0.005;
756	cPoints = {
757	ext::make_tuple(args&: lowerBound, args: lowerBoundDensity, args: false),
758	ext::make_tuple(args: v0Center, args: v0Density, args: true)
759	};
760	}
761	break;
762	default:
763	QL_FAIL("unknown transformation type");
764	}
765
766	const ext::shared_ptr<Fdm1dMesher> varianceMesher(
767	ext::make_shared<Concentrating1dMesher>(args&: lowerBound, args&: upperBound,
768	args: vGrid, args&: cPoints, args: 1e-12));
769
770	const Real sEps = 1e-4;
771	const Real sLowerBound
772	= std::log(x: hestonPxBoundary(maturity, eps: sEps, model));
773	const Real sUpperBound
774	= std::log(x: hestonPxBoundary(maturity, eps: 1-sEps,model));
775
776	const ext::shared_ptr<Fdm1dMesher> spotMesher(
777	ext::make_shared<Concentrating1dMesher>(args: sLowerBound, args: sUpperBound, args: xGrid,
778	args: std::make_pair(x: x0, y: 0.1), args: true));
779
780	const ext::shared_ptr<FdmMesherComposite>
781	mesher(ext::make_shared<FdmMesherComposite>(args: spotMesher, args: varianceMesher));
782
783	const ext::shared_ptr<FdmLinearOpComposite> hestonFwdOp(
784	ext::make_shared<FdmHestonFwdOp>(args: mesher, args: process, args: transformationType));
785
786	ModifiedCraigSneydScheme evolver(
787	FdmSchemeDesc::ModifiedCraigSneyd().theta,
788	FdmSchemeDesc::ModifiedCraigSneyd().mu, hestonFwdOp);
789
790	// step one days using non-correlated process
791	const Time eT = 1.0/365;
792	Array p = FdmHestonGreensFct(mesher, process, testCase.trafoType)
793	.get(t: eT, algorithm: testCase.greensAlgorithm);
794
795	const Real strikes[] = { 50, 80, 90, 100, 110, 120, 150, 200 };
796
797	Time t=eT;
798	for (auto maturitie : maturities) {
799
800	// calculate step size
801	const Date nextMaturityDate = todaysDate + maturitie;
802	const Time nextMaturityTime
803	= dc.yearFraction(d1: todaysDate, d2: nextMaturityDate);
804
805	Time dt = (nextMaturityTime - t)/tGridPerYear;
806	evolver.setStep(dt);
807
808	for (Size i=0; i < tGridPerYear; ++i, t+=dt) {
809	evolver.step(a&: p, t: t+dt);
810	}
811
812	Real avg=0, min=QL_MAX_REAL, max=0;
813	for (Real strike : strikes) {
814	const ext::shared_ptr<StrikedTypePayoff> payoff(
815	ext::make_shared<PlainVanillaPayoff>(args: (strike > s0) ? Option::Call : Option::Put,
816	args&: strike));
817
818	Array pd(p.size());
819	for (const auto& iter : *mesher->layout()) {
820	const Size idx = iter.index();
821	const Real s = std::exp(x: mesher->location(iter, direction: 0));
822
823	pd[idx] = (*payoff)(s)*p[idx];
824	if (transformationType == FdmSquareRootFwdOp::Power) {
825	const Real v = mesher->location(iter, direction: 1);
826	pd[idx] *= std::pow(x: v, y: -alpha);
827	}
828	}
829
830	const Real calculated = fokkerPlanckPrice2D(p: pd, mesher)
831	* rTS->discount(d: nextMaturityDate);
832
833	const ext::shared_ptr<Exercise> exercise(
834	ext::make_shared<EuropeanExercise>(args: nextMaturityDate));
835
836	VanillaOption option(payoff, exercise);
837	option.setPricingEngine(engine);
838
839	const Real expected = option.NPV();
840	const Real absDiff = std::fabs(x: expected - calculated);
841	const Real relDiff = absDiff / std::max(QL_EPSILON, b: expected);
842	const Real diff = std::min(a: absDiff, b: relDiff);
843
844	avg += diff;
845	min = std::min(a: diff, b: min);
846	max = std::max(a: diff, b: max);
847
848	if (diff > testCase.eps) {
849	BOOST_FAIL("failed to reproduce Heston SLV prices at"
850	<< "\n strike " << strike
851	<< "\n kappa " << kappa
852	<< "\n theta " << theta
853	<< "\n rho " << rho
854	<< "\n sigma " << sigma
855	<< "\n v0 " << v0
856	<< "\n transform " << transformationType
857	<< std::fixed << std::setprecision(5)
858	<< "\n calculated: " << calculated
859	<< "\n expected: " << expected
860	<< "\n tolerance: " << testCase.eps);
861	}
862	}
863
864	avg/=LENGTH(strikes); // NOLINT(bugprone-integer-division)
865
866	if (avg > testCase.avgEps) {
867	BOOST_FAIL("failed to reproduce Heston SLV prices"
868	" on average at"
869	<< "\n kappa " << kappa
870	<< "\n theta " << theta
871	<< "\n rho " << rho
872	<< "\n sigma " << sigma
873	<< "\n v0 " << v0
874	<< "\n transform " << transformationType
875	<< std::fixed << std::setprecision(5)
876	<< "\n average diff: " << avg
877	<< "\n tolerance: " << testCase.avgEps);
878	}
879	}
880	}
881	}
882
883	void HestonSLVModelTest::testHestonFokkerPlanckFwdEquation() {
884	BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation "
885	"for the Heston process...");
886
887	FokkerPlanckFwdTestCase testCases[] = {
888	{
889	.s0: 100.0, .r: 0.01, .q: 0.02,
890	// Feller constraint violated, high vol-of-vol case
891	// \frac{2\kappa\theta}{\sigma^2} = 2.0 * 1.0 * 0.05 / 0.2 = 0.5 < 1
892	.v0: 0.05, .kappa: 1.0, .theta: 0.05, .rho: -0.75, .sigma: std::sqrt(x: 0.2),
893	.xGrid: 101, .vGrid: 401, .tGridPerYear: 25, .tMinGridPerYear: 25,
894	.avgEps: 0.02, .eps: 0.05,
895	.trafoType: FdmSquareRootFwdOp::Power,
896	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
897	.schemeType: FdmSchemeDesc::DouglasType
898	},
899	{
900	.s0: 100.0, .r: 0.01, .q: 0.02,
901	// Feller constraint violated, high vol-of-vol case
902	// \frac{2\kappa\theta}{\sigma^2} = 2.0 * 1.0 * 0.05 / 0.2 = 0.5 < 1
903	.v0: 0.05, .kappa: 1.0, .theta: 0.05, .rho: -0.75, .sigma: std::sqrt(x: 0.2),
904	.xGrid: 201, .vGrid: 501, .tGridPerYear: 10, .tMinGridPerYear: 10,
905	.avgEps: 0.005, .eps: 0.02,
906	.trafoType: FdmSquareRootFwdOp::Log,
907	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
908	.schemeType: FdmSchemeDesc::HundsdorferType
909	},
910	{
911	.s0: 100.0, .r: 0.01, .q: 0.02,
912	// Feller constraint violated, high vol-of-vol case
913	// \frac{2\kappa\theta}{\sigma^2} = 2.0 * 1.0 * 0.05 / 0.2 = 0.5 < 1
914	.v0: 0.05, .kappa: 1.0, .theta: 0.05, .rho: -0.75, .sigma: std::sqrt(x: 0.2),
915	.xGrid: 201, .vGrid: 501, .tGridPerYear: 25, .tMinGridPerYear: 25,
916	.avgEps: 0.01, .eps: 0.03,
917	.trafoType: FdmSquareRootFwdOp::Log,
918	.greensAlgorithm: FdmHestonGreensFct::ZeroCorrelation,
919	.schemeType: FdmSchemeDesc::HundsdorferType
920	},
921	{
922	.s0: 100.0, .r: 0.01, .q: 0.02,
923	// Feller constraint fulfilled, low vol-of-vol case
924	// \frac{2\kappa\theta}{\sigma^2} = 2.0 * 1.0 * 0.05 / 0.05 = 2.0 > 1
925	.v0: 0.05, .kappa: 1.0, .theta: 0.05, .rho: -0.75, .sigma: std::sqrt(x: 0.05),
926	.xGrid: 201, .vGrid: 401, .tGridPerYear: 5, .tMinGridPerYear: 5,
927	.avgEps: 0.01, .eps: 0.02,
928	.trafoType: FdmSquareRootFwdOp::Plain,
929	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
930	.schemeType: FdmSchemeDesc::HundsdorferType
931	}
932	};
933
934	for (const auto& testCase : testCases) {
935	hestonFokkerPlanckFwdEquationTest(testCase);
936	}
937	}
938
939
940	namespace {
941	ext::shared_ptr<Matrix>
942	createLocalVolMatrixFromProcess(const ext::shared_ptr<BlackScholesMertonProcess>& lvProcess,
943	const std::vector<Real>& strikes,
944	const std::vector<Date>& dates,
945	std::vector<Time>& times) {
946
947	const ext::shared_ptr<LocalVolTermStructure> localVol =
948	lvProcess->localVolatility().currentLink();
949
950	const DayCounter dc = localVol->dayCounter();
951	const Date todaysDate = Settings::instance().evaluationDate();
952
953	QL_REQUIRE(times.size() == dates.size(), "mismatch");
954
955	for (Size i=0; i < times.size(); ++i) {
956	times[i] = dc.yearFraction(d1: todaysDate, d2: dates[i]);
957	}
958
959	ext::shared_ptr<Matrix> surface(
960	ext::make_shared<Matrix>(args: strikes.size(), args: dates.size()));
961
962	for (Size i=0; i < strikes.size(); ++i) {
963	for (Size j=0; j < dates.size(); ++j) {
964	try {
965	(*surface)[i][j] = localVol->localVol(d: dates[j], underlyingLevel: strikes[i], extrapolate: true);
966	} catch (Error&) {
967	(*surface)[i][j] = 0.2;
968	}
969	}
970	}
971
972	return surface;
973	}
974
975	ext::tuple<std::vector<Real>, std::vector<Date>,
976	ext::shared_ptr<BlackVarianceSurface> >
977	createSmoothImpliedVol(const DayCounter& dc, const Calendar& cal) {
978
979	const Date todaysDate = Settings::instance().evaluationDate();
980
981	Integer times[] = { 13, 41, 75, 165, 256, 345, 524, 703 };
982	std::vector<Date> dates;
983	for (int time : times) {
984	Date date = todaysDate + time;
985	dates.push_back(x: date);
986	}
987
988	const std::vector<Real> surfaceStrikes = {
989	2.222222222, 11.11111111, 44.44444444, 75.55555556, 80, 84.44444444, 88.88888889, 93.33333333, 97.77777778, 100,
990	102.2222222, 106.6666667, 111.1111111, 115.5555556, 120, 124.4444444, 166.6666667, 222.2222222, 444.4444444, 666.6666667
991	};
992
993	Volatility v[] =
994	{ 1.015873, 1.015873, 0.915873, 0.89729, 0.796493, 0.730914, 0.631335, 0.568895,
995	0.851309, 0.821309, 0.781309, 0.641309, 0.635593, 0.583653, 0.508045, 0.463182,
996	0.686034, 0.630534, 0.590534, 0.500534, 0.448706, 0.416661, 0.375470, 0.353442,
997	0.526034, 0.482263, 0.447713, 0.387703, 0.355064, 0.337438, 0.316966, 0.306859,
998	0.497587, 0.464373, 0.430764, 0.374052, 0.344336, 0.328607, 0.310619, 0.301865,
999	0.479511, 0.446815, 0.414194, 0.361010, 0.334204, 0.320301, 0.304664, 0.297180,
1000	0.461866, 0.429645, 0.398092, 0.348638, 0.324680, 0.312512, 0.299082, 0.292785,
1001	0.444801, 0.413014, 0.382634, 0.337026, 0.315788, 0.305239, 0.293855, 0.288660,
1002	0.428604, 0.397219, 0.368109, 0.326282, 0.307555, 0.298483, 0.288972, 0.284791,
1003	0.420971, 0.389782, 0.361317, 0.321274, 0.303697, 0.295302, 0.286655, 0.282948,
1004	0.413749, 0.382754, 0.354917, 0.316532, 0.300016, 0.292251, 0.284420, 0.281164,
1005	0.400889, 0.370272, 0.343525, 0.307904, 0.293204, 0.286549, 0.280189, 0.277767,
1006	0.390685, 0.360399, 0.334344, 0.300507, 0.287149, 0.281380, 0.276271, 0.274588,
1007	0.383477, 0.353434, 0.327580, 0.294408, 0.281867, 0.276746, 0.272655, 0.271617,
1008	0.379106, 0.349214, 0.323160, 0.289618, 0.277362, 0.272641, 0.269332, 0.268846,
1009	0.377073, 0.347258, 0.320776, 0.286077, 0.273617, 0.269057, 0.266293, 0.266265,
1010	0.399925, 0.369232, 0.338895, 0.289042, 0.265509, 0.255589, 0.249308, 0.249665,
1011	0.423432, 0.406891, 0.373720, 0.314667, 0.281009, 0.263281, 0.246451, 0.242166,
1012	0.453704, 0.453704, 0.453704, 0.381255, 0.334578, 0.305527, 0.268909, 0.251367,
1013	0.517748, 0.517748, 0.517748, 0.416577, 0.364770, 0.331595, 0.287423, 0.264285 };
1014
1015	Matrix blackVolMatrix(surfaceStrikes.size(), dates.size());
1016	for (Size i=0; i < surfaceStrikes.size(); ++i)
1017	for (Size j=0; j < dates.size(); ++j) {
1018	blackVolMatrix[i][j] = v[i*(dates.size())+j];
1019	}
1020
1021	const ext::shared_ptr<BlackVarianceSurface> volTS(
1022	ext::make_shared<BlackVarianceSurface>(args: todaysDate, args: cal,
1023	args&: dates,
1024	args: surfaceStrikes, args&: blackVolMatrix,
1025	args: dc,
1026	args: BlackVarianceSurface::ConstantExtrapolation,
1027	args: BlackVarianceSurface::ConstantExtrapolation));
1028	volTS->setInterpolation<Bicubic>();
1029
1030	return ext::make_tuple(args: surfaceStrikes, args&: dates, args: volTS);
1031	}
1032	}
1033
1034	void HestonSLVModelTest::testHestonFokkerPlanckFwdEquationLogLVLeverage() {
1035	BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation "
1036	"for the Heston process Log Transformation with leverage LV limiting case...");
1037
1038	const DayCounter dc = ActualActual(ActualActual::ISDA);
1039	const Date todaysDate = Date(28, Dec, 2012);
1040	Settings::instance().evaluationDate() = todaysDate;
1041
1042	const Date maturityDate = todaysDate + Period(1, Years);
1043	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
1044
1045	const Real s0 = 100;
1046	const Real x0 = std::log(x: s0);
1047	const Rate r = 0.0;
1048	const Rate q = 0.0;
1049
1050	const Real kappa = 1.0;
1051	const Real theta = 1.0;
1052	const Real rho = -0.75;
1053	const Real sigma = 0.02;
1054	const Real v0 = theta;
1055
1056	const FdmSquareRootFwdOp::TransformationType transform
1057	= FdmSquareRootFwdOp::Plain;
1058
1059	const DayCounter dayCounter = Actual365Fixed();
1060	const Calendar calendar = TARGET();
1061
1062	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1063	const Handle<YieldTermStructure> rTS(flatRate(today: todaysDate, forward: r, dc: dayCounter));
1064	const Handle<YieldTermStructure> qTS(flatRate(today: todaysDate, forward: q, dc: dayCounter));
1065
1066	ext::shared_ptr<HestonProcess> hestonProcess(
1067	ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1068
1069	const Size xGrid = 201;
1070	const Size vGrid = 401;
1071	const Size tGrid = 25;
1072
1073	const SquareRootProcessRNDCalculator rnd(v0, kappa, theta, sigma);
1074
1075	const Real upperBound = rnd.stationary_invcdf(q: 0.99);
1076	const Real lowerBound = rnd.stationary_invcdf(q: 0.01);
1077
1078	const Real beta = 10.0;
1079	std::vector<ext::tuple<Real, Real, bool>> critPoints = {
1080	ext::make_tuple(args: lowerBound, args: beta, args: true),
1081	ext::make_tuple(args: v0, args: beta/100, args: true),
1082	ext::make_tuple(args: upperBound, args: beta, args: true)
1083	};
1084	const ext::shared_ptr<Fdm1dMesher> varianceMesher(
1085	ext::make_shared<Concentrating1dMesher>(args: lowerBound, args: upperBound, args: vGrid, args&: critPoints));
1086
1087	const ext::shared_ptr<Fdm1dMesher> equityMesher(
1088	ext::make_shared<Concentrating1dMesher>(args: std::log(x: 2.0), args: std::log(x: 600.0), args: xGrid,
1089	args: std::make_pair(x: x0+0.005, y: 0.1), args: true));
1090
1091	const ext::shared_ptr<FdmMesherComposite>
1092	mesher(ext::make_shared<FdmMesherComposite>(args: equityMesher, args: varianceMesher));
1093
1094	const ext::tuple<std::vector<Real>, std::vector<Date>,
1095	ext::shared_ptr<BlackVarianceSurface> > smoothSurface =
1096	createSmoothImpliedVol(dc: dayCounter, cal: calendar);
1097	const ext::shared_ptr<BlackScholesMertonProcess> lvProcess(
1098	ext::make_shared<BlackScholesMertonProcess>(args: spot, args: qTS, args: rTS,
1099	args: Handle<BlackVolTermStructure>(ext::get<2>(t: smoothSurface))));
1100
1101	// step two days using non-correlated process
1102	const Time eT = 2.0/365;
1103
1104	Real v=-Null<Real>(), p_v(0.0);
1105	Array p(mesher->layout()->size(), 0.0);
1106	const Real bsV0 = squared(x: lvProcess->blackVolatility()->blackVol(t: 0.0, strike: s0, extrapolate: true));
1107
1108	SquareRootProcessRNDCalculator rndCalculator(v0, kappa, theta, sigma);
1109	for (const auto& iter : *mesher->layout()) {
1110	const Real x = mesher->location(iter, direction: 0);
1111	if (v != mesher->location(iter, direction: 1)) {
1112	v = mesher->location(iter, direction: 1);
1113	// the extreme tail probabilities of the non central
1114	// chi square distribution lead to numerical exceptions
1115	// on some platforms
1116	if (std::fabs(x: v - v0) < 5*sigma*std::sqrt(x: v0*eT))
1117	p_v = rndCalculator.pdf(v, t: eT);
1118	else
1119	p_v = 0.0;
1120	}
1121	const Real p_x = 1.0/(std::sqrt(M_TWOPI*bsV0*eT))
1122	* std::exp(x: -0.5*squared(x: x - x0)/(bsV0*eT));
1123	p[iter.index()] = p_v*p_x;
1124	}
1125	const Time dt = (maturity-eT)/tGrid;
1126
1127
1128	const std::vector<Real> denseStrikes =
1129	{ 2.222222222, 11.11111111, 20, 25, 30, 35, 40,
1130	44.44444444, 50, 55, 60, 65, 70, 75.55555556,
1131	80, 84.44444444, 88.88888889, 93.33333333, 97.77777778, 100,
1132	102.2222222, 106.6666667, 111.1111111, 115.5555556, 120,
1133	124.4444444, 166.6666667, 222.2222222, 444.4444444, 666.6666667 };
1134
1135	Matrix surface(denseStrikes.size(), ext::get<1>(t: smoothSurface).size());
1136	std::vector<Time> times(surface.columns());
1137
1138	const std::vector<Date>& dates = ext::get<1>(t: smoothSurface);
1139	ext::shared_ptr<Matrix> m = createLocalVolMatrixFromProcess(
1140	lvProcess, strikes: denseStrikes, dates, times);
1141
1142	const ext::shared_ptr<FixedLocalVolSurface> leverage(
1143	ext::make_shared<FixedLocalVolSurface>(args: todaysDate, args: dates, args: denseStrikes, args&: m, args: dc));
1144
1145	const ext::shared_ptr<PricingEngine> lvEngine(
1146	ext::make_shared<AnalyticEuropeanEngine>(args: lvProcess));
1147
1148	const ext::shared_ptr<FdmLinearOpComposite> hestonFwdOp(
1149	ext::make_shared<FdmHestonFwdOp>(args: mesher, args&: hestonProcess, args: transform, args: leverage));
1150
1151	HundsdorferScheme evolver(FdmSchemeDesc::Hundsdorfer().theta,
1152	FdmSchemeDesc::Hundsdorfer().mu,
1153	hestonFwdOp);
1154
1155	Time t=dt;
1156	evolver.setStep(dt);
1157
1158	for (Size i=0; i < tGrid; ++i, t+=dt) {
1159	evolver.step(a&: p, t);
1160	}
1161
1162	const ext::shared_ptr<Exercise> exercise(
1163	ext::make_shared<EuropeanExercise>(args: maturityDate));
1164
1165	const ext::shared_ptr<PricingEngine> fdmEngine(
1166	ext::make_shared<FdBlackScholesVanillaEngine>(args: lvProcess, args: 50, args: 201, args: 0,
1167	args: FdmSchemeDesc::Douglas(), args: true,args: 0.2));
1168
1169	for (Size strike=5; strike < 200; strike+=10) {
1170	const ext::shared_ptr<StrikedTypePayoff> payoff(
1171	ext::make_shared<CashOrNothingPayoff>(args: Option::Put, args: Real(strike), args: 1.0));
1172
1173	Array pd(p.size());
1174	for (const auto& iter : *mesher->layout()) {
1175	const Size idx = iter.index();
1176	const Real s = std::exp(x: mesher->location(iter, direction: 0));
1177
1178	pd[idx] = (*payoff)(s)*p[idx];
1179	}
1180
1181	const Real calculated
1182	= fokkerPlanckPrice2D(p: pd, mesher)*rTS->discount(d: maturityDate);
1183
1184	VanillaOption option(payoff, exercise);
1185	option.setPricingEngine(fdmEngine);
1186	const Real expected = option.NPV();
1187
1188	const Real tol = 0.015;
1189	if (std::fabs(x: expected - calculated ) > tol) {
1190	BOOST_FAIL("failed to reproduce Heston prices at"
1191	<< "\n strike " << strike
1192	<< std::fixed << std::setprecision(5)
1193	<< "\n calculated: " << calculated
1194	<< "\n expected: " << expected
1195	<< "\n tolerance: " << tol);
1196	}
1197	}
1198	}
1199
1200
1201	void HestonSLVModelTest::testBlackScholesFokkerPlanckFwdEquationLocalVol() {
1202	BOOST_TEST_MESSAGE(
1203	"Testing Fokker-Planck forward equation for BS Local Vol process...");
1204
1205	const DayCounter dc = ActualActual(ActualActual::ISDA);
1206	const Date todaysDate(5, July, 2014);
1207	Settings::instance().evaluationDate() = todaysDate;
1208
1209	const Real s0 = 100;
1210	const Real x0 = std::log(x: s0);
1211	const Rate r = 0.035;
1212	const Rate q = 0.01;
1213
1214	const Calendar calendar = TARGET();
1215	const DayCounter dayCounter = Actual365Fixed();
1216
1217	const Handle<YieldTermStructure> rTS(
1218	flatRate(today: todaysDate, forward: r, dc: dayCounter));
1219	const Handle<YieldTermStructure> qTS(
1220	flatRate(today: todaysDate, forward: q, dc: dayCounter));
1221
1222	ext::tuple<std::vector<Real>, std::vector<Date>,
1223	ext::shared_ptr<BlackVarianceSurface> > smoothImpliedVol =
1224	createSmoothImpliedVol(dc: dayCounter, cal: calendar);
1225
1226	const std::vector<Real>& strikes = ext::get<0>(t&: smoothImpliedVol);
1227	const std::vector<Date>& dates = ext::get<1>(t&: smoothImpliedVol);
1228	const Handle<BlackVolTermStructure> vTS =
1229	Handle<BlackVolTermStructure>(ext::get<2>(t&: smoothImpliedVol));
1230
1231	const Size xGrid = 101;
1232	const Size tGrid = 51;
1233
1234	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1235	const ext::shared_ptr<BlackScholesMertonProcess> process(
1236	ext::make_shared<BlackScholesMertonProcess>(args: spot, args: qTS, args: rTS, args: vTS));
1237
1238	const ext::shared_ptr<LocalVolTermStructure> localVol(
1239	ext::make_shared<NoExceptLocalVolSurface>(args: vTS, args: rTS, args: qTS, args: spot, args: 0.2));
1240
1241	const ext::shared_ptr<PricingEngine> engine(
1242	ext::make_shared<AnalyticEuropeanEngine>(args: process));
1243
1244	for (Size i = 1; i < dates.size(); i += 2) {
1245	for (Size j = 3; j < strikes.size() - 3; j += 2) {
1246	const Date& exDate = dates[i];
1247	const Date maturityDate = exDate;
1248	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
1249	const ext::shared_ptr<Exercise> exercise(
1250	ext::make_shared<EuropeanExercise>(args: exDate));
1251
1252	const ext::shared_ptr<FdmMesher> uniformMesher(
1253	ext::make_shared<FdmMesherComposite>(
1254	args: ext::make_shared<FdmBlackScholesMesher>(args: xGrid, args: process,
1255	args: maturity, args: s0)));
1256
1257	//-- operator --
1258	const ext::shared_ptr<FdmLinearOpComposite> uniformBSFwdOp(
1259	ext::make_shared<FdmLocalVolFwdOp>(
1260	args: uniformMesher, args: *spot, args: *rTS, args: *qTS, args: localVol));
1261
1262	const ext::shared_ptr<FdmMesher> concentratedMesher(
1263	ext::make_shared<FdmMesherComposite>(
1264	args: ext::make_shared<FdmBlackScholesMesher>(args: xGrid, args: process,
1265	args: maturity, args: s0, args: Null<Real>(),
1266	args: Null<Real>(), args: 0.0001, args: 1.5,
1267	args: std::pair<Real, Real>(s0, 0.1))));
1268
1269	//-- operator --
1270	const ext::shared_ptr<FdmLinearOpComposite> concentratedBSFwdOp(
1271	ext::make_shared<FdmLocalVolFwdOp>(
1272	args: concentratedMesher, args: *spot, args: *rTS, args: *qTS, args: localVol));
1273
1274	const ext::shared_ptr<FdmMesher> shiftedMesher(
1275	ext::make_shared<FdmMesherComposite>(
1276	args: ext::make_shared<FdmBlackScholesMesher>(args: xGrid, args: process,
1277	args: maturity, args: s0, args: Null<Real>(),
1278	args: Null<Real>(), args: 0.0001, args: 1.5,
1279	args: std::pair<Real, Real>(s0 * 1.1,
1280	0.2))));
1281
1282	//-- operator --
1283	const ext::shared_ptr<FdmLinearOpComposite> shiftedBSFwdOp(
1284	ext::make_shared<FdmLocalVolFwdOp>(
1285	args: shiftedMesher, args: *spot, args: *rTS, args: *qTS, args: localVol));
1286
1287	const ext::shared_ptr<StrikedTypePayoff> payoff(
1288	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strikes[j]));
1289
1290	VanillaOption option(payoff, exercise);
1291	option.setPricingEngine(engine);
1292
1293	const Real expected = option.NPV();
1294	const Real calcUniform = fokkerPlanckPrice1D(mesher: uniformMesher,
1295	op: uniformBSFwdOp, payoff, x0, maturity, tGrid)
1296	* rTS->discount(d: maturityDate);
1297	const Real calcConcentrated = fokkerPlanckPrice1D(
1298	mesher: concentratedMesher, op: concentratedBSFwdOp, payoff, x0,
1299	maturity, tGrid) * rTS->discount(d: maturityDate);
1300	const Real calcShifted = fokkerPlanckPrice1D(mesher: shiftedMesher,
1301	op: shiftedBSFwdOp, payoff, x0, maturity, tGrid)
1302	* rTS->discount(d: maturityDate);
1303	const Real tol = 0.05;
1304
1305	if (std::fabs(x: expected - calcUniform) > tol) {
1306	BOOST_FAIL("failed to reproduce european option price "
1307	<< "with an uniform mesher"
1308	<< "\n strike: " << strikes[i]
1309	<< std::fixed << std::setprecision(8)
1310	<< "\n calculated: " << calcUniform
1311	<< "\n expected: " << expected
1312	<< "\n tolerance: " << tol);
1313	}
1314	if (std::fabs(x: expected - calcConcentrated) > tol) {
1315	BOOST_FAIL("failed to reproduce european option price "
1316	<< "with a concentrated mesher"
1317	<< "\n strike: " << strikes[i]
1318	<< std::fixed << std::setprecision(8)
1319	<< "\n calculated: " << calcConcentrated
1320	<< "\n expected: " << expected
1321	<< "\n tolerance: " << tol);
1322	}
1323	if (std::fabs(x: expected - calcShifted) > tol) {
1324	BOOST_FAIL("failed to reproduce european option price "
1325	<< "with a shifted mesher"
1326	<< "\n strike: " << strikes[i]
1327	<< std::fixed << std::setprecision(8)
1328	<< "\n calculated: " << calcShifted
1329	<< "\n expected: " << expected
1330	<< "\n tolerance: " << tol);
1331	}
1332	}
1333	}
1334	}
1335
1336
1337	namespace {
1338	struct HestonModelParams {
1339	const Rate r, q;
1340	const Real kappa, theta, rho, sigma, v0;
1341	};
1342
1343	struct HestonSLVTestCase {
1344	const HestonModelParams hestonParams;
1345	const HestonSLVFokkerPlanckFdmParams fdmParams;
1346	};
1347
1348
1349	void lsvCalibrationTest(const HestonSLVTestCase& testCase) {
1350	const Date todaysDate(2, June, 2015);
1351	Settings::instance().evaluationDate() = todaysDate;
1352	const Date finalDate(2, June, 2020);
1353
1354	const DayCounter dc = Actual365Fixed();
1355
1356	const Real s0 = 100;
1357	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1358
1359	const Rate r = testCase.hestonParams.r;
1360	const Rate q = testCase.hestonParams.q;
1361
1362	const Real kappa = testCase.hestonParams.kappa;
1363	const Real theta = testCase.hestonParams.theta;
1364	const Real rho = testCase.hestonParams.rho;
1365	const Real sigma = testCase.hestonParams.sigma;
1366	const Real v0 = testCase.hestonParams.v0;
1367	const Volatility lv = 0.3;
1368
1369	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
1370	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
1371
1372	const ext::shared_ptr<HestonProcess> hestonProcess(
1373	ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1374
1375	const Handle<HestonModel> hestonModel(
1376	ext::make_shared<HestonModel>(args: hestonProcess));
1377
1378	const Handle<LocalVolTermStructure> localVol(
1379	ext::make_shared<LocalConstantVol>(args: todaysDate, args: lv, args: dc));
1380
1381	const HestonSLVFDMModel slvModel(
1382	localVol, hestonModel, finalDate, testCase.fdmParams);
1383
1384	// this includes a calibration of the leverage function!
1385	ext::shared_ptr<LocalVolTermStructure>
1386	l = slvModel.leverageFunction();
1387
1388	const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
1389	ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS,
1390	args: Handle<BlackVolTermStructure>(flatVol(volatility: lv, dc))));
1391
1392	const ext::shared_ptr<PricingEngine> analyticEngine(
1393	ext::make_shared<AnalyticEuropeanEngine>(args: bsProcess));
1394
1395	const Real strikes[] = { 50, 75, 80, 90, 100, 110, 125, 150 };
1396	const Size times[] = { 3, 6, 9, 12, 24, 36, 60 };
1397
1398	for (unsigned long time : times) {
1399	const Date expiry = todaysDate + Period(time, Months);
1400	const ext::shared_ptr<Exercise> exercise(
1401	ext::make_shared<EuropeanExercise>(args: expiry));
1402
1403	const ext::shared_ptr<PricingEngine> slvEngine(
1404	(time <= 3) ?
1405	ext::make_shared<FdHestonVanillaEngine>(
1406	args: hestonModel.currentLink(), args: Size(std::max(a: 101.0, b: 51 * time / 12.0)), args: 401,
1407	args: 101, args: 0, args: FdmSchemeDesc::ModifiedCraigSneyd(), args&: l) :
1408	ext::make_shared<FdHestonVanillaEngine>(
1409	args: hestonModel.currentLink(), args: Size(std::max(a: 51.0, b: 51 * time / 12.0)), args: 201, args: 101,
1410	args: 0, args: FdmSchemeDesc::ModifiedCraigSneyd(), args&: l));
1411
1412	for (Real strike : strikes) {
1413	const ext::shared_ptr<StrikedTypePayoff> payoff(
1414	ext::make_shared<PlainVanillaPayoff>(args: (strike > s0) ? Option::Call : Option::Put,
1415	args&: strike));
1416
1417	VanillaOption option(payoff, exercise);
1418
1419	option.setPricingEngine(slvEngine);
1420	const Real calculated = option.NPV();
1421
1422	option.setPricingEngine(analyticEngine);
1423	const Real expected = option.NPV();
1424	const Real vega = option.vega();
1425
1426	const ext::shared_ptr<GeneralizedBlackScholesProcess> bp(
1427	ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS,
1428	args: Handle<BlackVolTermStructure>(flatVol(volatility: lv,
1429	dc))));
1430
1431	const Real tol = 0.001;//testCase.eps;
1432	if (std::fabs(x: (calculated-expected)/vega) > tol) {
1433	BOOST_FAIL("failed to reproduce round trip vola "
1434	<< "\n strike " << strike << "\n time " << time
1435	<< "\n expected NPV " << expected << "\n calculated NPV "
1436	<< calculated << "\n vega " << vega << std::fixed
1437	<< std::setprecision(5)
1438	<< "\n calculated: " << lv + (calculated - expected) / vega
1439	<< "\n expected: " << lv << "\n diff (in bp) "
1440	<< (calculated - expected) / vega * 1e4
1441	<< "\n tolerance: " << tol);
1442	}
1443	}
1444	}
1445	}
1446	}
1447
1448
1449
1450	void HestonSLVModelTest::testFDMCalibration() {
1451	const HestonSLVFokkerPlanckFdmParams plainParams =
1452	{ .xGrid: 201, .vGrid: 301, .tMaxStepsPerYear: 1000, .tMinStepsPerYear: 25, .tStepNumberDecay: 3.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
1453	.x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
1454	.vLowerEps: 1e-8, .vUpperEps: 1e-8, .vMin: 0.0, .v0Density: 1.0, .vLowerBoundDensity: 1.0, .vUpperBoundDensity: 1.0, .leverageFctPropEps: 1e-6,
1455	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
1456	.trafoType: FdmSquareRootFwdOp::Plain,
1457	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
1458	};
1459
1460	const HestonSLVFokkerPlanckFdmParams logParams =
1461	{ .xGrid: 301, .vGrid: 601, .tMaxStepsPerYear: 2000, .tMinStepsPerYear: 30, .tStepNumberDecay: 2.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
1462	.x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
1463	.vLowerEps: 1e-5, .vUpperEps: 1e-5, .vMin: 0.0000025, .v0Density: 1.0, .vLowerBoundDensity: 0.1, .vUpperBoundDensity: 0.9, .leverageFctPropEps: 1e-5,
1464	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
1465	.trafoType: FdmSquareRootFwdOp::Log,
1466	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
1467	};
1468
1469	const HestonSLVFokkerPlanckFdmParams powerParams =
1470	{ .xGrid: 401, .vGrid: 801, .tMaxStepsPerYear: 2000, .tMinStepsPerYear: 30, .tStepNumberDecay: 2.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
1471	.x0Density: 0.1, .localVolEpsProb: 1e-3, .maxIntegrationIterations: 10000,
1472	.vLowerEps: 1e-6, .vUpperEps: 1e-6, .vMin: 0.001, .v0Density: 1.0, .vLowerBoundDensity: 0.001, .vUpperBoundDensity: 1.0, .leverageFctPropEps: 1e-5,
1473	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
1474	.trafoType: FdmSquareRootFwdOp::Power,
1475	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
1476	};
1477
1478
1479	const HestonSLVTestCase testCases[] = {
1480	{ .hestonParams: {.r: 0.035, .q: 0.01, .kappa: 1.0, .theta: 0.06, .rho: -0.75, .sigma: 0.1, .v0: 0.09}, .fdmParams: plainParams},
1481	{ .hestonParams: {.r: 0.035, .q: 0.01, .kappa: 1.0, .theta: 0.06, .rho: -0.75, .sigma: std::sqrt(x: 0.2), .v0: 0.09}, .fdmParams: logParams },
1482	{ .hestonParams: {.r: 0.035, .q: 0.01, .kappa: 1.0, .theta: 0.09, .rho: -0.75, .sigma: std::sqrt(x: 0.2), .v0: 0.06}, .fdmParams: logParams },
1483	{ .hestonParams: {.r: 0.035, .q: 0.01, .kappa: 1.0, .theta: 0.06, .rho: -0.75, .sigma: 0.2, .v0: 0.09}, .fdmParams: powerParams }
1484	};
1485
1486
1487	for (Size i=0; i < LENGTH(testCases); ++i) {
1488	BOOST_TEST_MESSAGE("Testing stochastic local volatility calibration case " << i << " ...");
1489	lsvCalibrationTest(testCase: testCases[i]);
1490	}
1491
1492	}
1493
1494	void HestonSLVModelTest::testLocalVolsvSLVPropDensity() {
1495	BOOST_TEST_MESSAGE("Testing local volatility vs SLV model...");
1496
1497	const Date todaysDate(5, Oct, 2015);
1498	const Date finalDate = todaysDate + Period(1, Years);
1499	Settings::instance().evaluationDate() = todaysDate;
1500
1501	const Real s0 = 100;
1502	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1503	const Rate r = 0.01;
1504	const Rate q = 0.02;
1505
1506	const Calendar calendar = TARGET();
1507	const DayCounter dayCounter = Actual365Fixed();
1508
1509	const Handle<YieldTermStructure> rTS(
1510	flatRate(today: todaysDate, forward: r, dc: dayCounter));
1511	const Handle<YieldTermStructure> qTS(
1512	flatRate(today: todaysDate, forward: q, dc: dayCounter));
1513
1514	const Handle<BlackVolTermStructure> vTS = Handle<BlackVolTermStructure>(
1515	ext::get<2>(t: createSmoothImpliedVol(dc: dayCounter, cal: calendar)));
1516
1517	// Heston parameter from implied calibration
1518	const Real kappa = 2.0;
1519	const Real theta = 0.074;
1520	const Real rho = -0.51;
1521	const Real sigma = 0.8;
1522	const Real v0 = 0.1974;
1523
1524	const ext::shared_ptr<HestonProcess> hestonProcess(
1525	ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1526
1527	const Handle<HestonModel> hestonModel(
1528	ext::make_shared<HestonModel>(args: hestonProcess));
1529
1530	const Handle<LocalVolTermStructure> localVol(
1531	ext::make_shared<NoExceptLocalVolSurface>(args: vTS, args: rTS, args: qTS, args: spot, args: 0.3));
1532	localVol->enableExtrapolation(b: true);
1533
1534	const Size vGrid = 151;
1535	const Size xGrid = 51;
1536
1537	const HestonSLVFokkerPlanckFdmParams fdmParams = {
1538	.xGrid: xGrid, .vGrid: vGrid, .tMaxStepsPerYear: 500, .tMinStepsPerYear: 50, .tStepNumberDecay: 100.0, .nRannacherTimeSteps: 5, .predictionCorretionSteps: 2,
1539	.x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
1540	.vLowerEps: 1e-5, .vUpperEps: 1e-5, .vMin: 0.0000025,
1541	.v0Density: 1.0, .vLowerBoundDensity: 0.1, .vUpperBoundDensity: 0.9, .leverageFctPropEps: 1e-5,
1542	.greensAlgorithm: FdmHestonGreensFct::ZeroCorrelation,
1543	.trafoType: FdmSquareRootFwdOp::Log,
1544	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
1545	};
1546
1547	const HestonSLVFDMModel slvModel(
1548	localVol, hestonModel, finalDate, fdmParams, true);
1549
1550	const std::list<HestonSLVFDMModel::LogEntry>& logEntries
1551	= slvModel.logEntries();
1552
1553	const SquareRootProcessRNDCalculator squareRootRndCalculator(
1554	v0, kappa, theta, sigma);
1555
1556	for (const auto& logEntrie : logEntries) {
1557
1558	const Time t = logEntrie.t;
1559	if (t > 0.2) {
1560	const Array x(logEntrie.mesher->getFdm1dMeshers().at(n: 0)->locations().begin(),
1561	logEntrie.mesher->getFdm1dMeshers().at(n: 0)->locations().end());
1562	const std::vector<Real>& z = logEntrie.mesher->getFdm1dMeshers().at(n: 1)->locations();
1563
1564	const ext::shared_ptr<Array>& prob = logEntrie.prob;
1565
1566	for (Size i=0; i < z.size(); ++i) {
1567	const Real pCalc = DiscreteSimpsonIntegral()(
1568	x, Array(prob->begin()+i*xGrid,
1569	prob->begin()+(i+1)*xGrid));
1570
1571	const Real expected
1572	= squareRootRndCalculator.pdf(v: std::exp(x: z[i]), t);
1573	const Real calculated = pCalc/std::exp(x: z[i]);
1574
1575	if ( std::fabs(x: expected-calculated) > 0.01
1576	&& std::fabs(x: (expected-calculated)/expected) > 0.04) {
1577	BOOST_ERROR("failed to reproduce probability at "
1578	<< "\n v : " << std::exp(z[i])
1579	<< "\n t : " << t
1580	<< "\n expected : " << expected
1581	<< "\n calculated : " << calculated);
1582	}
1583	}
1584	}
1585	}
1586	}
1587
1588
1589	void HestonSLVModelTest::testBarrierPricingViaHestonLocalVol() {
1590	BOOST_TEST_MESSAGE("Testing calibration via vanilla options...");
1591
1592	const DayCounter dc = ActualActual(ActualActual::ISDA);
1593	const Date todaysDate(5, Nov, 2015);
1594	Settings::instance().evaluationDate() = todaysDate;
1595
1596	const Real s0 = 100;
1597	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1598	const Rate r = 0.1;
1599	const Rate q = 0.025;
1600
1601	const Real kappa = 2.0;
1602	const Real theta = 0.09;
1603	const Real rho = -0.75;
1604	const Real sigma = 0.8;
1605	const Real v0 = 0.19;
1606
1607	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
1608	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
1609
1610	const ext::shared_ptr<HestonProcess> hestonProcess(
1611	ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1612
1613	const Handle<HestonModel> hestonModel(
1614	ext::make_shared<HestonModel>(args: hestonProcess));
1615
1616	const Handle<BlackVolTermStructure> surf(
1617	ext::make_shared<HestonBlackVolSurface>(args: hestonModel));
1618
1619	const Real strikeValues[] = { 50, 75, 100, 125, 150, 200, 400 };
1620	const Period maturities[] = {
1621	Period(1, Months), Period(2,Months),
1622	Period(3, Months), Period(4, Months),
1623	Period(5, Months), Period(6, Months),
1624	Period(9, Months), Period(1, Years),
1625	Period(18, Months), Period(2, Years),
1626	Period(3, Years), Period(5, Years) };
1627
1628	const ext::shared_ptr<LocalVolSurface> localVolSurface(
1629	ext::make_shared<LocalVolSurface>(args: surf, args: rTS, args: qTS, args: spot));
1630
1631	const ext::shared_ptr<PricingEngine> hestonEngine(
1632	ext::make_shared<AnalyticHestonEngine>(args: hestonModel.currentLink(), args: 164));
1633
1634	for (Real strike : strikeValues) {
1635	for (auto maturitie : maturities) {
1636	const Date exerciseDate = todaysDate + maturitie;
1637	const Time t = dc.yearFraction(d1: todaysDate, d2: exerciseDate);
1638
1639	const Volatility impliedVol = surf->blackVol(t, strike, extrapolate: true);
1640
1641	const ext::shared_ptr<GeneralizedBlackScholesProcess>
1642	bsProcess(ext::make_shared<GeneralizedBlackScholesProcess>(
1643	args: spot, args: qTS, args: rTS,
1644	args: Handle<BlackVolTermStructure>(flatVol(volatility: impliedVol, dc))));
1645
1646	const ext::shared_ptr<PricingEngine> analyticEngine(
1647	ext::make_shared<AnalyticEuropeanEngine>(args: bsProcess));
1648
1649	const ext::shared_ptr<Exercise> exercise(
1650	ext::make_shared<EuropeanExercise>(args: exerciseDate));
1651	const ext::shared_ptr<StrikedTypePayoff> payoff(
1652	ext::make_shared<PlainVanillaPayoff>(
1653	args: spot->value() < strike ? Option::Call : Option::Put,
1654	args&: strike));
1655
1656	const ext::shared_ptr<PricingEngine> localVolEngine(
1657	ext::make_shared<FdBlackScholesVanillaEngine>(args: bsProcess, args: 201, args: 801, args: 0,
1658	args: FdmSchemeDesc::Douglas(), args: true));
1659
1660	VanillaOption option(payoff, exercise);
1661
1662	option.setPricingEngine(analyticEngine);
1663	const Real analyticNPV = option.NPV();
1664
1665	option.setPricingEngine(hestonEngine);
1666	const Real hestonNPV = option.NPV();
1667
1668	option.setPricingEngine(localVolEngine);
1669	const Real localVolNPV = option.NPV();
1670
1671	const Real tol = 1e-3;
1672	if (std::fabs(x: analyticNPV - hestonNPV) > tol)
1673	BOOST_ERROR("Heston and BS price do not match "
1674	<< "\n Heston : " << hestonNPV
1675	<< "\n Black-Scholes: " << analyticNPV
1676	<< "\n diff : "
1677	<< std::fabs(analyticNPV - hestonNPV));
1678
1679	if (std::fabs(x: analyticNPV - localVolNPV) > tol)
1680	BOOST_ERROR("LocalVol and BS price do not match "
1681	<< "\n LocalVol : " << localVolNPV
1682	<< "\n Black-Scholes: " << analyticNPV
1683	<< "\n diff : "
1684	<< std::fabs(analyticNPV - localVolNPV));
1685	}
1686	}
1687	}
1688
1689	void HestonSLVModelTest::testBarrierPricingMixedModels() {
1690	BOOST_TEST_MESSAGE("Testing Barrier pricing with mixed models...");
1691
1692	const DayCounter dc = ActualActual(ActualActual::ISDA);
1693	const Date todaysDate(5, Nov, 2015);
1694	const Date exerciseDate = todaysDate + Period(1, Years);
1695	Settings::instance().evaluationDate() = todaysDate;
1696
1697	const Real s0 = 100;
1698	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1699	const Rate r = 0.05;
1700	const Rate q = 0.02;
1701
1702	const Real kappa = 2.0;
1703	const Real theta = 0.09;
1704	const Real rho = -0.75;
1705	const Real sigma = 0.4;
1706	const Real v0 = 0.19;
1707
1708	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
1709	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
1710
1711	const ext::shared_ptr<HestonProcess> hestonProcess(
1712	ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
1713
1714	const Handle<HestonModel> hestonModel(
1715	ext::make_shared<HestonModel>(args: hestonProcess));
1716
1717	const Handle<BlackVolTermStructure> impliedVolSurf(
1718	ext::make_shared<HestonBlackVolSurface>(args: hestonModel));
1719
1720	const Handle<LocalVolTermStructure> localVolSurf(
1721	ext::make_shared<NoExceptLocalVolSurface>(
1722	args: impliedVolSurf, args: rTS, args: qTS, args: spot, args: 0.3));
1723
1724	const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
1725	ext::make_shared<GeneralizedBlackScholesProcess>(
1726	args: spot, args: qTS, args: rTS, args: impliedVolSurf));
1727
1728	const ext::shared_ptr<Exercise> exercise(
1729	ext::make_shared<EuropeanExercise>(args: exerciseDate));
1730	const ext::shared_ptr<StrikedTypePayoff> payoff(
1731	ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: s0));
1732
1733	const ext::shared_ptr<PricingEngine> hestonEngine(
1734	ext::make_shared<FdHestonBarrierEngine>(
1735	args: hestonModel.currentLink(), args: 26, args: 101, args: 51));
1736
1737	const ext::shared_ptr<PricingEngine> localEngine(
1738	ext::make_shared<FdBlackScholesBarrierEngine>(
1739	args: bsProcess, args: 26, args: 101, args: 0, args: FdmSchemeDesc::Douglas(), args: true, args: 0.3));
1740
1741	const Real barrier = 10.0;
1742	BarrierOption barrierOption(
1743	Barrier::DownOut, barrier, 0.0, payoff, exercise);
1744
1745	barrierOption.setPricingEngine(hestonEngine);
1746	const Real hestonDeltaCalculated = barrierOption.delta();
1747
1748	barrierOption.setPricingEngine(localEngine);
1749	const Real localDeltaCalculated = barrierOption.delta();
1750
1751	const Real localDeltaExpected = -0.439068;
1752	const Real hestonDeltaExpected = -0.342059;
1753	const Real tol = 0.001;
1754	if (std::fabs(x: hestonDeltaExpected - hestonDeltaCalculated) > tol) {
1755	BOOST_ERROR("Heston Delta does not match"
1756	<< "\n calculated : " << hestonDeltaCalculated
1757	<< "\n expected : " << hestonDeltaExpected);
1758	}
1759	if (std::fabs(x: localDeltaExpected - localDeltaCalculated) > tol) {
1760	BOOST_ERROR("Local Vol Delta does not match"
1761	<< "\n calculated : " << localDeltaCalculated
1762	<< "\n expected : " << localDeltaExpected);
1763	}
1764
1765	const HestonSLVFokkerPlanckFdmParams params =
1766	{ .xGrid: 51, .vGrid: 201, .tMaxStepsPerYear: 1000, .tMinStepsPerYear: 100, .tStepNumberDecay: 3.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
1767	.x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
1768	.vLowerEps: 1e-8, .vUpperEps: 1e-8, .vMin: 0.0, .v0Density: 1.0, .vLowerBoundDensity: 1.0, .vUpperBoundDensity: 1.0, .leverageFctPropEps: 1e-6,
1769	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
1770	.trafoType: FdmSquareRootFwdOp::Plain,
1771	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
1772	};
1773
1774	const Real eta[] = { 0.1, 0.2, 0.3, 0.4,
1775	0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
1776
1777	const Real slvDeltaExpected[] = {
1778	-0.429475, -0.419749, -0.410055, -0.400339,
1779	-0.390616, -0.380888, -0.371156, -0.361425,
1780	-0.351699, -0.341995 };
1781
1782	for (Size i=0; i < LENGTH(eta); ++i) {
1783	const Handle<HestonModel> modHestonModel(
1784	ext::make_shared<HestonModel>(
1785	args: ext::make_shared<HestonProcess>(
1786	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: eta[i]*sigma, args: rho)));
1787
1788	const HestonSLVFDMModel slvModel(
1789	localVolSurf, modHestonModel, exerciseDate, params);
1790
1791	const ext::shared_ptr<LocalVolTermStructure>
1792	leverageFct = slvModel.leverageFunction();
1793
1794	const ext::shared_ptr<PricingEngine> slvEngine(
1795	ext::make_shared<FdHestonBarrierEngine>(
1796	args: modHestonModel.currentLink(), args: 201, args: 801, args: 201, args: 0,
1797	args: FdmSchemeDesc::Hundsdorfer(), args: leverageFct));
1798
1799	BarrierOption barrierOption(
1800	Barrier::DownOut, barrier, 0.0, payoff, exercise);
1801
1802	barrierOption.setPricingEngine(slvEngine);
1803	const Real slvDeltaCalculated = barrierOption.delta();
1804
1805	if (std::fabs(x: slvDeltaExpected[i] - slvDeltaCalculated) > tol) {
1806	BOOST_ERROR("Stochastic Local Vol Delta does not match"
1807	<< "\n calculated : " << slvDeltaCalculated
1808	<< "\n expected : " << slvDeltaExpected);
1809	}
1810	}
1811	}
1812
1813	void HestonSLVModelTest::testMonteCarloVsFdmPricing() {
1814	BOOST_TEST_MESSAGE(
1815	"Testing Monte-Carlo vs FDM Pricing for "
1816	"Heston SLV models...");
1817
1818	const DayCounter dc = ActualActual(ActualActual::ISDA);
1819	const Date todaysDate(5, Dec, 2015);
1820	const Date exerciseDate = todaysDate + Period(1, Years);
1821	Settings::instance().evaluationDate() = todaysDate;
1822
1823	const Real s0 = 100;
1824	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1825	const Rate r = 0.05;
1826	const Rate q = 0.02;
1827
1828	const Real kappa = 2.0;
1829	const Real theta = 0.18;
1830	const Real rho = -0.75;
1831	const Real sigma = 0.8;
1832	const Real v0 = 0.19;
1833
1834	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
1835	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
1836
1837	const ext::shared_ptr<HestonProcess> hestonProcess
1838	= ext::make_shared<HestonProcess>(
1839	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho);
1840
1841	const ext::shared_ptr<HestonModel> hestonModel
1842	= ext::make_shared<HestonModel>(args: hestonProcess);
1843
1844	const ext::shared_ptr<LocalVolTermStructure> leverageFct
1845	= ext::make_shared<LocalConstantVol>(args: todaysDate, args: 0.25, args: dc);
1846
1847	const ext::shared_ptr<HestonSLVProcess> slvProcess
1848	= ext::make_shared<HestonSLVProcess>(args: hestonProcess, args: leverageFct);
1849
1850	const ext::shared_ptr<PricingEngine> mcEngine
1851	= MakeMCEuropeanHestonEngine<
1852	PseudoRandom, GeneralStatistics, HestonSLVProcess>(slvProcess)
1853	.withStepsPerYear(steps: 100)
1854	.withAntitheticVariate()
1855	.withSamples(samples: 10000)
1856	.withSeed(seed: 1234);
1857
1858	const ext::shared_ptr<PricingEngine> fdEngine
1859	= ext::make_shared<FdHestonVanillaEngine>(
1860	args: hestonModel, args: 51, args: 401, args: 101, args: 0,
1861	args: FdmSchemeDesc::ModifiedCraigSneyd(), args: leverageFct);
1862
1863	const ext::shared_ptr<Exercise> exercise
1864	= ext::make_shared<EuropeanExercise>(args: exerciseDate);
1865
1866	const ext::shared_ptr<HestonProcess> mixingProcess
1867	= ext::make_shared<HestonProcess>(args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma * 10, args: rho,
1868	args: HestonProcess::QuadraticExponentialMartingale);
1869	const ext::shared_ptr<HestonModel> mixingModel
1870	= ext::make_shared<HestonModel>(args: mixingProcess);
1871
1872	const ext::shared_ptr<PricingEngine> fdEngineWithMixingFactor
1873	= ext::make_shared<FdHestonVanillaEngine>(
1874	args: mixingModel, args: 51, args: 401, args: 101, args: 0,
1875	args: FdmSchemeDesc::ModifiedCraigSneyd(), args: leverageFct, args: 0.1);
1876
1877	const Real strikes[] = { s0, 1.1*s0 };
1878	for (Real strike : strikes) {
1879	const ext::shared_ptr<StrikedTypePayoff> payoff =
1880	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args&: strike);
1881
1882	VanillaOption option(payoff, exercise);
1883
1884	option.setPricingEngine(fdEngine);
1885
1886	const Real priceFDM = option.NPV();
1887
1888	option.setPricingEngine(fdEngineWithMixingFactor);
1889
1890	const Real priceFDMWithMix = option.NPV();
1891
1892	option.setPricingEngine(mcEngine);
1893
1894	const Real priceMC = option.NPV();
1895	const Real priceError = option.errorEstimate();
1896
1897	if (priceError > 0.1) {
1898	BOOST_ERROR("Heston Monte-Carlo error is too large"
1899	<< "\n MC Error: " << priceError
1900	<< "\n Limit : " << 0.1);
1901	}
1902
1903	if (std::fabs(x: priceFDM - priceMC) > 2.3*priceError) {
1904	BOOST_ERROR("Heston Monte-Carlo price does not match with FDM"
1905	<< "\n MC Price : " << priceMC
1906	<< "\n MC Error : " << priceError
1907	<< "\n FDM Price: " << priceFDM);
1908	}
1909
1910	if (priceFDM != priceFDMWithMix) {
1911	BOOST_ERROR("Heston mixing FDM price does not match with non-mixing FDM"
1912	<< "\n Mixing FDM Price : " << priceFDMWithMix
1913	<< "\n Non Mixing FDM Price : " << priceFDM);
1914	}
1915	}
1916	}
1917
1918	void HestonSLVModelTest::testMonteCarloCalibration() {
1919	BOOST_TEST_MESSAGE(
1920	"Testing Monte-Carlo Calibration...");
1921
1922	const DayCounter dc = ActualActual(ActualActual::ISDA);
1923	const Date todaysDate(5, Jan, 2016);
1924	const Date maturityDate = todaysDate + Period(1, Years);
1925	Settings::instance().evaluationDate() = todaysDate;
1926
1927	const Real s0 = 100;
1928	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
1929	const Rate r = 0.05;
1930	const Rate q = 0.02;
1931
1932	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
1933	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
1934
1935	const ext::shared_ptr<LocalVolTermStructure> localVol
1936	= ext::make_shared<LocalConstantVol>(args: todaysDate, args: 0.3, args: dc);
1937
1938	// parameter of the "calibrated" Heston model
1939	const Real kappa = 1.0;
1940	const Real theta = 0.06;
1941	const Real rho = -0.75;
1942	const Real sigma = 0.4;
1943	const Real v0 = 0.09;
1944
1945	const ext::shared_ptr<HestonProcess> hestonProcess
1946	= ext::make_shared<HestonProcess>(
1947	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho);
1948
1949	const ext::shared_ptr<HestonModel> hestonModel
1950	= ext::make_shared<HestonModel>(args: hestonProcess);
1951
1952	const Size xGrid = 400;
1953	const Size nSims[] = { 40000 };
1954
1955	for (unsigned long nSim : nSims) {
1956	const bool sobol = true;
1957
1958	const ext::shared_ptr<LocalVolTermStructure> leverageFct =
1959	HestonSLVMCModel(
1960	Handle<LocalVolTermStructure>(localVol), Handle<HestonModel>(hestonModel),
1961	sobol ? ext::shared_ptr<BrownianGeneratorFactory>(new SobolBrownianGeneratorFactory(
1962	SobolBrownianGenerator::Diagonal, 1234UL, SobolRsg::JoeKuoD7)) :
1963	ext::shared_ptr<BrownianGeneratorFactory>(
1964	new MTBrownianGeneratorFactory(1234UL)),
1965	maturityDate, 91, xGrid, nSim)
1966	.leverageFunction();
1967
1968	const ext::shared_ptr<PricingEngine> bsEngine(
1969	ext::make_shared<AnalyticEuropeanEngine>(
1970	args: ext::make_shared<GeneralizedBlackScholesProcess>(
1971	args: spot, args: qTS, args: rTS,
1972	args: Handle<BlackVolTermStructure>(flatVol(volatility: 0.3, dc)))));
1973
1974	const Real strikes[] = { 50, 80, 100, 120, 150, 200 };
1975	const Date maturities[] = {
1976	todaysDate + Period(3, Months), todaysDate + Period(6, Months),
1977	todaysDate + Period(12, Months)
1978	};
1979
1980	Real qualityFactor = 0.0;
1981	Real maxQualityFactor = 0.0;
1982	Size nValues = 0U;
1983
1984	for (auto maturity : maturities) {
1985	const Time maturityTime = dc.yearFraction(d1: todaysDate, d2: maturity);
1986
1987	const ext::shared_ptr<PricingEngine> fdEngine
1988	= ext::make_shared<FdHestonVanillaEngine>(
1989	args: hestonModel, args: std::max(a: Size(26), b: Size(maturityTime*51)),
1990	args: 201, args: 51, args: 0,
1991	args: FdmSchemeDesc::ModifiedCraigSneyd(), args: leverageFct);
1992
1993	const ext::shared_ptr<Exercise> exercise
1994	= ext::make_shared<EuropeanExercise>(args&: maturity);
1995
1996	for (Real strike : strikes) {
1997	const ext::shared_ptr<StrikedTypePayoff> payoff =
1998	ext::make_shared<PlainVanillaPayoff>(args: strike < s0 ? Option::Put : Option::Call,
1999	args&: strike);
2000
2001	VanillaOption option(payoff, exercise);
2002
2003	option.setPricingEngine(bsEngine);
2004	const Real bsNPV = option.NPV();
2005	const Real bsVega = option.vega();
2006
2007	if (bsNPV > 0.02) {
2008	option.setPricingEngine(fdEngine);
2009	const Real fdmNPV = option.NPV();
2010
2011	const Real diff = std::fabs(x: fdmNPV-bsNPV)/bsVega*1e4;
2012
2013	qualityFactor+=diff;
2014	maxQualityFactor = std::max(a: maxQualityFactor, b: diff);
2015	++nValues;
2016	}
2017	}
2018	}
2019
2020	if (qualityFactor/nValues > 7.5) {
2021	BOOST_ERROR(
2022	"Failed to reproduce average calibration quality"
2023	<< "\n average calibration quality : "
2024	<< qualityFactor/nValues << "bp"
2025	<< "\n tolerance : 5.0bp");
2026	}
2027
2028	if (qualityFactor/nValues > 15.0) {
2029	BOOST_ERROR(
2030	"Failed to reproduce maximum calibration error"
2031	<< "\n maximum calibration error : "
2032	<< maxQualityFactor << "bp"
2033	<< "\n tolerance : 15.0bp");
2034	}
2035	}
2036	}
2037
2038
2039	void HestonSLVModelTest::testForwardSkewSLV() {
2040	BOOST_TEST_MESSAGE("Testing the implied volatility skew of "
2041	"forward starting options in SLV model...");
2042
2043	const DayCounter dc = ActualActual(ActualActual::ISDA);
2044	const Date todaysDate(5, Jan, 2017);
2045	const Date maturityDate = todaysDate + Period(2, Years);
2046	Settings::instance().evaluationDate() = todaysDate;
2047
2048	const Real s0 = 100;
2049	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
2050	const Rate r = 0.05;
2051	const Rate q = 0.02;
2052	const Volatility flatLocalVol = 0.3;
2053
2054	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
2055	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
2056
2057	const Handle<LocalVolTermStructure> localVol(
2058	ext::make_shared<LocalConstantVol>(args: todaysDate, args: flatLocalVol, args: dc));
2059
2060	// parameter of the "calibrated" Heston model
2061	const Real kappa = 2.0;
2062	const Real theta = 0.06;
2063	const Real rho = -0.75;
2064	const Real sigma = 0.6;
2065	const Real v0 = 0.09;
2066
2067	const ext::shared_ptr<HestonProcess> hestonProcess
2068	= ext::make_shared<HestonProcess>(
2069	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho);
2070
2071	const Handle<HestonModel> hestonModel(
2072	ext::make_shared<HestonModel>(args: hestonProcess));
2073
2074
2075	// Monte-Carlo calibration
2076	const Size nSim = 40000;
2077	const Size xGrid = 200;
2078
2079	const ext::shared_ptr<LocalVolTermStructure> leverageFctMC =
2080	HestonSLVMCModel(
2081	localVol, hestonModel,
2082	ext::shared_ptr<BrownianGeneratorFactory>(new MTBrownianGeneratorFactory(1234UL)),
2083	maturityDate, 182, xGrid, nSim)
2084	.leverageFunction();
2085
2086	const ext::shared_ptr<HestonSLVProcess> mcSlvProcess(
2087	ext::make_shared<HestonSLVProcess>(args: hestonProcess, args: leverageFctMC));
2088
2089	// finite difference calibration
2090	const HestonSLVFokkerPlanckFdmParams logParams = {
2091	.xGrid: 201, .vGrid: 401, .tMaxStepsPerYear: 1000, .tMinStepsPerYear: 30, .tStepNumberDecay: 2.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
2092	.x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
2093	.vLowerEps: 1e-5, .vUpperEps: 1e-5, .vMin: 0.0000025, .v0Density: 1.0, .vLowerBoundDensity: 0.1, .vUpperBoundDensity: 0.9, .leverageFctPropEps: 1e-5,
2094	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
2095	.trafoType: FdmSquareRootFwdOp::Log,
2096	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
2097	};
2098
2099	const ext::shared_ptr<LocalVolTermStructure> leverageFctFDM =
2100	HestonSLVFDMModel(
2101	localVol, hestonModel, maturityDate, logParams).leverageFunction();
2102
2103	const ext::shared_ptr<HestonSLVProcess> fdmSlvProcess(
2104	ext::make_shared<HestonSLVProcess>(args: hestonProcess, args: leverageFctFDM));
2105
2106	const Date resetDate = todaysDate + Period(12, Months);
2107	const Time resetTime = dc.yearFraction(d1: todaysDate, d2: resetDate);
2108	const Time maturityTime = dc.yearFraction(d1: todaysDate, d2: maturityDate);
2109	std::vector<Time> mandatoryTimes = {resetTime, maturityTime};
2110
2111	const Size tSteps = 100;
2112	const TimeGrid grid(mandatoryTimes.begin(), mandatoryTimes.end(), tSteps);
2113	const Size resetIndex = grid.closestIndex(t: resetTime);
2114
2115	typedef SobolBrownianBridgeRsg rsg_type;
2116	typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
2117
2118	const Size factors = mcSlvProcess->factors();
2119
2120	std::vector<ext::shared_ptr<MultiPathGenerator<rsg_type> > > pathGen = {
2121	ext::make_shared<MultiPathGenerator<rsg_type> >(
2122	args: mcSlvProcess, args: grid, args: rsg_type(factors, tSteps), args: false),
2123	ext::make_shared<MultiPathGenerator<rsg_type> >(
2124	args: fdmSlvProcess, args: grid, args: rsg_type(factors, tSteps), args: false)
2125	};
2126
2127	const Real strikes[] = {
2128	0.5, 0.7, 0.8, 0.9, 1.0, 1.1, 1.25, 1.5, 1.75, 2.0
2129	};
2130
2131	std::vector<std::vector<GeneralStatistics> >
2132	stats(2, std::vector<GeneralStatistics>(LENGTH(strikes)));
2133
2134	for (Size i=0; i < 5*nSim; ++i) {
2135	for (Size k= 0; k < 2; ++k) {
2136	const sample_type& path = pathGen[k]->next();
2137
2138	const Real S_t1 = path.value[0][resetIndex-1];
2139	const Real S_T1 = path.value[0][tSteps-1];
2140
2141	const sample_type& antiPath = pathGen[k]->antithetic();
2142	const Real S_t2 = antiPath.value[0][resetIndex-1];
2143	const Real S_T2 = antiPath.value[0][tSteps-1];
2144
2145	for (Size j=0; j < LENGTH(strikes); ++j) {
2146	const Real strike = strikes[j];
2147	if (strike < 1.0)
2148	stats[k][j].add(value: 0.5*(
2149	S_t1 * std::max(a: 0.0, b: strike - S_T1/S_t1)
2150	+ S_t2 * std::max(a: 0.0, b: strike - S_T2/S_t2)));
2151	else
2152	stats[k][j].add(value: 0.5*(
2153	S_t1 * std::max(a: 0.0, b: S_T1/S_t1 - strike)
2154	+ S_t2 * std::max(a: 0.0, b: S_T2/S_t2 - strike)));
2155	}
2156	}
2157	}
2158
2159	const ext::shared_ptr<Exercise> exercise(
2160	ext::make_shared<EuropeanExercise>(args: maturityDate));
2161
2162	const ext::shared_ptr<SimpleQuote> vol(
2163	ext::make_shared<SimpleQuote>(args: flatLocalVol));
2164	const Handle<BlackVolTermStructure> volTS(flatVol(today: todaysDate, volatility: vol, dc));
2165
2166	const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
2167	ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS, args: volTS));
2168
2169	const ext::shared_ptr<PricingEngine> fwdEngine(
2170	ext::make_shared<ForwardVanillaEngine<AnalyticEuropeanEngine> >(
2171	args: bsProcess));
2172
2173	const Volatility expected[] = {
2174	0.37804, 0.346608, 0.330682, 0.314978, 0.300399,
2175	0.287273, 0.272916, 0.26518, 0.268663, 0.277052
2176	};
2177
2178	const DiscountFactor df = rTS->discount(t: grid.back());
2179
2180	for (Size j=0; j < LENGTH(strikes); ++j) {
2181	for (Size k= 0; k < 2; ++k) {
2182	const Real strike = strikes[j];
2183	const Real npv = stats[k][j].mean() * df;
2184
2185	const ext::shared_ptr<StrikedTypePayoff> payoff(
2186	ext::make_shared<PlainVanillaPayoff>(
2187	args: (strike < 1.0) ? Option::Put : Option::Call, args: strike));
2188
2189	const ext::shared_ptr<ForwardVanillaOption> fwdOption(
2190	ext::make_shared<ForwardVanillaOption>(
2191	args: strike, args: resetDate, args: payoff, args: exercise));
2192
2193	const Volatility implVol =
2194	QuantLib::detail::ImpliedVolatilityHelper::calculate(
2195	instrument: *fwdOption, engine: *fwdEngine, volQuote&: *vol, targetValue: npv, accuracy: 1e-8, maxEvaluations: 200, minVol: 1e-4, maxVol: 2.0);
2196
2197	const Real tol = 0.002;
2198	const Volatility volError = std::fabs(x: implVol - expected[j]);
2199
2200	if (volError > tol) {
2201	BOOST_ERROR("Implied forward volatility error is too large"
2202	<< "\n expected forward volatility: " << expected[j]
2203	<< "\n SLV forward volatility : " << implVol
2204	<< "\n difference : " << volError
2205	<< "\n tolerance : " << tol
2206	<< "\n calibration method : " <<
2207	((k) ? "Monte-Carlo" : "Finite Difference"));
2208	}
2209	}
2210	}
2211	}
2212
2213	namespace {
2214	ext::shared_ptr<LocalVolTermStructure> getFixedLocalVolFromHeston(
2215	const ext::shared_ptr<HestonModel>& hestonModel,
2216	const ext::shared_ptr<TimeGrid>& timeGrid) {
2217
2218	const Handle<BlackVolTermStructure> trueImpliedVolSurf(
2219	ext::make_shared<HestonBlackVolSurface>(
2220	args: Handle<HestonModel>(hestonModel),
2221	args: AnalyticHestonEngine::AndersenPiterbarg,
2222	args: AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder: 32)));
2223
2224	const ext::shared_ptr<HestonProcess> hestonProcess
2225	= hestonModel->process();
2226
2227	const ext::shared_ptr<LocalVolTermStructure> localVol(
2228	ext::make_shared<NoExceptLocalVolSurface>(
2229	args: trueImpliedVolSurf,
2230	args: hestonProcess->riskFreeRate(),
2231	args: hestonProcess->dividendYield(),
2232	args: hestonProcess->s0(),
2233	args: std::sqrt(x: hestonProcess->theta())));
2234
2235
2236	const ext::shared_ptr<LocalVolRNDCalculator> localVolRND(
2237	ext::make_shared<LocalVolRNDCalculator>(
2238	args: hestonProcess->s0().currentLink(),
2239	args: hestonProcess->riskFreeRate().currentLink(),
2240	args: hestonProcess->dividendYield().currentLink(),
2241	args: localVol,
2242	args: timeGrid));
2243
2244	std::vector<ext::shared_ptr<std::vector<Real> > > strikes;
2245	for (Size i=1; i < timeGrid->size(); ++i) {
2246	const Time t = timeGrid->at(i);
2247	const ext::shared_ptr<Fdm1dMesher> fdm1dMesher
2248	= localVolRND->mesher(t);
2249
2250	const std::vector<Real>& logStrikes = fdm1dMesher->locations();
2251	const ext::shared_ptr<std::vector<Real> > strikeSlice(
2252	ext::make_shared<std::vector<Real> >(args: logStrikes.size()));
2253
2254	for (Size j=0; j < logStrikes.size(); ++j) {
2255	(*strikeSlice)[j] = std::exp(x: logStrikes[j]);
2256	}
2257
2258	strikes.push_back(x: strikeSlice);
2259	}
2260
2261	const Size nStrikes = strikes.front()->size();
2262	const ext::shared_ptr<Matrix> localVolMatrix(
2263	ext::make_shared<Matrix>(args: nStrikes, args: timeGrid->size()-1));
2264	for (Size i=1; i < timeGrid->size(); ++i) {
2265	const Time t = timeGrid->at(i);
2266	const ext::shared_ptr<std::vector<Real> > strikeSlice
2267	= strikes[i-1];
2268
2269	for (Size j=0; j < nStrikes; ++j) {
2270	const Real s = (*strikeSlice)[j];
2271	(*localVolMatrix)[j][i-1] = localVol->localVol(t, underlyingLevel: s, extrapolate: true);
2272	}
2273	}
2274
2275	const Date todaysDate
2276	= hestonProcess->riskFreeRate()->referenceDate();
2277	const DayCounter dc = hestonProcess->riskFreeRate()->dayCounter();
2278
2279	const std::vector<Time> expiries(
2280	timeGrid->begin()+1, timeGrid->end());
2281
2282	return ext::make_shared<FixedLocalVolSurface>(
2283	args: todaysDate, args: expiries, args&: strikes, args: localVolMatrix, args: dc);
2284	}
2285	}
2286
2287	void HestonSLVModelTest::testMoustacheGraph() {
2288	BOOST_TEST_MESSAGE(
2289	"Testing double no touch pricing with SLV and mixing...");
2290
2291	/*
2292	A more detailed description of this test case can found on
2293	https://hpcquantlib.wordpress.com/2016/01/10/monte-carlo-calibration-of-the-heston-stochastic-local-volatiltiy-model/
2294
2295	The comparison of Black-Scholes and SLV prices is derived
2296	from figure 8.8 in Iain J. Clark's book,
2297	Foreign Exchange Option Pricing: A Practitioner’s Guide
2298	*/
2299
2300	const DayCounter dc = ActualActual(ActualActual::ISDA);
2301	const Date todaysDate(5, Jan, 2016);
2302	const Date maturityDate = todaysDate + Period(1, Years);
2303	Settings::instance().evaluationDate() = todaysDate;
2304
2305	const Real s0 = 100;
2306	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
2307	const Rate r = 0.02;
2308	const Rate q = 0.01;
2309
2310	// parameter of the "calibrated" Heston model
2311	const Real kappa = 1.0;
2312	const Real theta = 0.06;
2313	const Real rho = -0.8;
2314	const Real sigma = 0.8;
2315	const Real v0 = 0.09;
2316
2317	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
2318	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
2319
2320	const ext::shared_ptr<HestonModel> hestonModel(
2321	ext::make_shared<HestonModel>(
2322	args: ext::make_shared<HestonProcess>(
2323	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
2324
2325	const ext::shared_ptr<Exercise> europeanExercise(
2326	ext::make_shared<EuropeanExercise>(args: maturityDate));
2327
2328	VanillaOption vanillaOption(
2329	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0),
2330	europeanExercise);
2331
2332	vanillaOption.setPricingEngine(
2333	ext::make_shared<AnalyticHestonEngine>(args: hestonModel));
2334
2335	const Volatility implVol = vanillaOption.impliedVolatility(
2336	price: vanillaOption.NPV(),
2337	process: ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS,
2338	args: Handle<BlackVolTermStructure>(flatVol(volatility: std::sqrt(x: theta), dc))));
2339
2340	const ext::shared_ptr<PricingEngine> analyticEngine(
2341	ext::make_shared<AnalyticDoubleBarrierBinaryEngine>(
2342	args: ext::make_shared<GeneralizedBlackScholesProcess>(
2343	args: spot, args: qTS, args: rTS,
2344	args: Handle<BlackVolTermStructure>(flatVol(volatility: implVol, dc)))));
2345
2346	std::vector<Time> expiries;
2347	const Period timeStepPeriod = Period(1, Weeks);
2348	for (Date expiry = todaysDate + timeStepPeriod;
2349	expiry <= maturityDate; expiry+=timeStepPeriod) {
2350	expiries.push_back(x: dc.yearFraction(d1: todaysDate, d2: expiry));
2351	}
2352
2353	const ext::shared_ptr<TimeGrid> timeGrid(
2354	ext::make_shared<TimeGrid>(args: expiries.begin(), args: expiries.end()));
2355
2356	// first build the true local vol surface from another Heston model
2357	const Handle<LocalVolTermStructure> localVol(
2358	getFixedLocalVolFromHeston(hestonModel, timeGrid));
2359
2360	const ext::shared_ptr<BrownianGeneratorFactory> sobolGeneratorFactory(
2361	ext::make_shared<SobolBrownianGeneratorFactory>(args: SobolBrownianGenerator::Diagonal, args: 1234UL,
2362	args: SobolRsg::JoeKuoD7));
2363
2364	const Size xGrid = 100;
2365	const Size nSim = 20000;
2366
2367	const Real eta = 0.90;
2368
2369	const Handle<HestonModel> modHestonModel(
2370	ext::make_shared<HestonModel>(
2371	args: ext::make_shared<HestonProcess>(
2372	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: eta*sigma, args: rho)));
2373
2374	const ext::shared_ptr<LocalVolTermStructure> leverageFct =
2375	HestonSLVMCModel(localVol, modHestonModel, sobolGeneratorFactory,
2376	maturityDate, 182, xGrid, nSim)
2377	.leverageFunction();
2378
2379	const ext::shared_ptr<PricingEngine> fdEngine(
2380	ext::make_shared<FdHestonDoubleBarrierEngine>(
2381	args: modHestonModel.currentLink(), args: 51, args: 101, args: 31, args: 0,
2382	args: FdmSchemeDesc::Hundsdorfer(), args: leverageFct));
2383
2384	const Real expected[] = {
2385	0.0334, 0.1141, 0.1319, 0.0957, 0.0464, 0.0058,-0.0192,
2386	-0.0293,-0.0297,-0.0251,-0.0192,-0.0134,-0.0084,-0.0045,
2387	-0.0015, 0.0005, 0.0017, 0.0020
2388	};
2389	const Real tol = 1e-2;
2390
2391	for (Size i=0; i < 18; ++i) {
2392	const Real dist = 10.0+5.0*i;
2393
2394	const Real barrier_lo = std::max(a: s0 - dist, b: 1e-2);
2395	const Real barrier_hi = s0 + dist;
2396	DoubleBarrierOption doubleBarrier(
2397	DoubleBarrier::KnockOut, barrier_lo, barrier_hi, 0.0,
2398	ext::make_shared<CashOrNothingPayoff>(
2399	args: Option::Call, args: 0.0, args: 1.0),
2400	europeanExercise);
2401
2402	doubleBarrier.setPricingEngine(analyticEngine);
2403	const Real bsNPV = doubleBarrier.NPV();
2404
2405	doubleBarrier.setPricingEngine(fdEngine);
2406	const Real slvNPV = doubleBarrier.NPV();
2407
2408	const Real diff = slvNPV - bsNPV;
2409	if (std::fabs(x: diff - expected[i]) > tol) {
2410	BOOST_ERROR(
2411	"Failed to reproduce price difference for a Double-No-Touch "
2412	<< "option between Black-Scholes and "
2413	<< "Heston Stochastic Local Volatility model"
2414	<< "\n Barrier Low : " << barrier_lo
2415	<< "\n Barrier High : " << barrier_hi
2416	<< "\n Black-Scholes Price: " << bsNPV
2417	<< "\n Heston SLV Price : " << slvNPV
2418	<< "\n diff : " << diff
2419	<< "\n expected diff : " << expected[i]
2420	<< "\n tolerance : " << tol);
2421	}
2422	}
2423	}
2424
2425	void HestonSLVModelTest::testDiffusionAndDriftSlvProcess() {
2426	BOOST_TEST_MESSAGE(
2427	"Testing diffusion and drift of the SLV process...");
2428
2429	const Date todaysDate(6, June, 2020);
2430	Settings::instance().evaluationDate() = todaysDate;
2431
2432	const DayCounter dc = Actual365Fixed();
2433	const Date maturityDate = todaysDate + Period(6, Months);
2434	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
2435
2436	const Real s0 = 100;
2437	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
2438	const Rate r = -0.005;
2439	const Rate q = 0.04;
2440
2441	const Handle<YieldTermStructure> rTS(flatRate(today: todaysDate, forward: r, dc));
2442	const Handle<YieldTermStructure> qTS(flatRate(today: todaysDate, forward: q, dc));
2443
2444	const ext::shared_ptr<LocalVolTermStructure> localVol =
2445	getFixedLocalVolFromHeston(
2446	hestonModel: ext::make_shared<HestonModel>(
2447	args: ext::make_shared<HestonProcess>(
2448	args: rTS, args: qTS, args: spot, args: 0.1, args: 1.0, args: 0.13, args: 0.8, args: 0.4)),
2449	timeGrid: ext::make_shared<TimeGrid>(args: maturity, args: 20));
2450
2451	const Real kappa = 2.5;
2452	const Real theta = 1.0;
2453	const Real rho = -0.75;
2454	const Real sigma = 2.4;
2455	const Real v0 = 1.0;
2456
2457	const ext::shared_ptr<HestonProcess> hestonProcess =
2458	ext::make_shared<HestonProcess>(
2459	args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho);
2460
2461	const Handle<HestonModel> hestonModel(
2462	ext::make_shared<HestonModel>(args: hestonProcess));
2463
2464	const ext::shared_ptr<HestonSLVProcess> slvProcess =
2465	ext::make_shared<HestonSLVProcess>(args: hestonProcess, args: localVol);
2466
2467	VanillaOption option(
2468	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0),
2469	ext::make_shared<EuropeanExercise>(args: maturityDate));
2470
2471	option.setPricingEngine(
2472	ext::make_shared<FdHestonVanillaEngine>(
2473	args: hestonModel.currentLink(),
2474	args: 26, args: 201, args: 101, args: 0,
2475	args: FdmSchemeDesc::ModifiedCraigSneyd(),
2476	args: localVol));
2477
2478	const Real expected = option.NPV();
2479
2480	const Size nSims = 16733;
2481	const Size nTimeSteps = 40;
2482	const DiscountFactor df = rTS->discount(t: maturity);
2483
2484	SobolBrownianBridgeRsg rsg(2, nTimeSteps, SobolBrownianGenerator::Diagonal, 12345U,
2485	SobolRsg::JoeKuoD7);
2486
2487	Array x(2), xt(2), dw(2);
2488	GeneralStatistics stats;
2489
2490	const Time dt = maturity/nTimeSteps;
2491	const Real sqrtDt = std::sqrt(x: dt);
2492
2493	for (Size i=0; i < nSims; ++i) {
2494	Time t = 0.0;
2495	x[0] = s0; x[1] = v0;
2496
2497	const std::vector<Real> n = rsg.nextSequence().value;
2498
2499	for (Size j=0; j < nTimeSteps; ++j, t+=dt) {
2500
2501	dw[0] = n[j];
2502	dw[1] = n[j+nTimeSteps];
2503
2504	// full truncation scheme
2505	xt[0] = x[0];
2506	xt[1] = (x[1] > 0)? x[1] : 0.0;
2507
2508	x = slvProcess->apply(x0: x,
2509	dx: slvProcess->diffusion(t, x: xt)*sqrtDt*dw
2510	+ slvProcess->drift(t, x: xt)*dt);
2511	}
2512
2513	stats.add(value: df*option.payoff()->operator()(price: x[0]));
2514	}
2515
2516	const Real calculated = stats.mean();
2517	const Real errorEstimate = stats.errorEstimate();
2518
2519	const Real diff = std::fabs(x: expected - calculated);
2520
2521	if (diff > 2.35*errorEstimate) {
2522	BOOST_ERROR(
2523	"Failed to reproduce call option price with HestonSLVProcess "
2524	"diffusion and drift discretization scheme"
2525	<< "\n expected : " << expected
2526	<< "\n calculated : " << calculated
2527	<< "\n error est. : " << errorEstimate
2528	<< "\n diff : " << diff);
2529	}
2530	}
2531
2532	void HestonSLVModelTest::testBarrierPricingMixedModelsMonteCarloVsFdmPricing() {
2533	BOOST_TEST_MESSAGE(
2534	"Testing European and Barrier Pricing for Monte-Carlo and FDM "
2535	"Pricing in Heston SLV models with a mixing factor...");
2536
2537	const Real epsilon = 0.015;
2538
2539	const DayCounter dc = ActualActual(ActualActual::ISDA);
2540	const Date todaysDate(1, Jul, 2021);
2541	const Date maturityDate = todaysDate + Period(2, Years);
2542	const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
2543	Settings::instance().evaluationDate() = todaysDate;
2544
2545	const Real s0 = 100;
2546	const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
2547	const Rate r = 0.02;
2548	const Rate q = 0.01;
2549	const Real mixingFactors[] = {1.0, 0.64, 0.3};
2550	const std::vector<Date>& requiredDates = std::vector<Date>();
2551
2552	// Create two slightly different Heston models. The first will be our stochastic
2553	// vol model, the second is used to create a similar implied vol surface which
2554	// we fit a local vol model to
2555	const Real kappa1 = 2.0;
2556	const Real theta1 = 0.12;
2557	const Real rho1 = -0.25;
2558	const Real sigma1 = 0.8;
2559	const Real v01 = 0.09;
2560
2561	const Real kappa2 = 1.5;
2562	const Real theta2 = 0.11;
2563	const Real rho2 = -0.2;
2564	const Real sigma2 = 0.9;
2565	const Real v02 = 0.1;
2566
2567	const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
2568	const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
2569
2570	const ext::shared_ptr<HestonProcess> hestonProcess
2571	= ext::make_shared<HestonProcess>(
2572	args: rTS, args: qTS, args: spot, args: v01, args: kappa1, args: theta1, args: sigma1, args: rho1);
2573
2574	const ext::shared_ptr<HestonModel> hestonModelPtr
2575	= ext::make_shared<HestonModel>(args: hestonProcess);
2576
2577	const ext::shared_ptr<HestonProcess> hestonProcess2
2578	= ext::make_shared<HestonProcess>(
2579	args: rTS, args: qTS, args: spot, args: v02, args: kappa2, args: theta2, args: sigma2, args: rho2);
2580
2581	const ext::shared_ptr<HestonModel> hestonModelPtr2
2582	= ext::make_shared<HestonModel>(args: hestonProcess2);
2583
2584	const ext::shared_ptr<LocalVolTermStructure> localVolPtr =
2585	getFixedLocalVolFromHeston(hestonModel: hestonModelPtr2,
2586	timeGrid: ext::make_shared<TimeGrid>(args: maturity, args: 20));
2587
2588	const Handle<LocalVolTermStructure> localVol = Handle<LocalVolTermStructure>(localVolPtr);
2589	localVol->enableExtrapolation();
2590	const Handle<HestonModel> hestonModel = Handle<HestonModel>(hestonModelPtr);
2591	const Handle<HestonModel> hestonModel2 = Handle<HestonModel>(hestonModelPtr2);
2592
2593	// Create the options we will price - a vanilla and a barrier
2594	const ext::shared_ptr<Exercise> exercise
2595	= ext::make_shared<EuropeanExercise>(args: maturityDate);
2596
2597	const Real strike = 100;
2598	const ext::shared_ptr<StrikedTypePayoff> payoff =
2599	ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike);
2600
2601	VanillaOption vanillaOption(payoff, exercise);
2602
2603	const Real rebate = 0.0;
2604	const Real barrier = 110.0;
2605	BarrierOption barrierOption(Barrier::UpOut, barrier, rebate, payoff, exercise);
2606
2607	// hestonModel2 is our simulated local vol model, so its vanilla prices
2608	// should match the calibrated SLV model pricers
2609	const ext::shared_ptr<PricingEngine> hestonVanillaEngine
2610	= ext::make_shared<AnalyticHestonEngine>(args: hestonModelPtr2);
2611	vanillaOption.setPricingEngine(hestonVanillaEngine);
2612	const Real localVolPrice = vanillaOption.NPV();
2613
2614	const ext::shared_ptr<BrownianGeneratorFactory> sobolGeneratorFactory(
2615	ext::make_shared<SobolBrownianGeneratorFactory>(args: SobolBrownianGenerator::Diagonal, args: 1234UL,
2616	args: SobolRsg::JoeKuoD7));
2617
2618	for (Real mixingFactor : mixingFactors) {
2619
2620	// Finite Difference calibration
2621	const HestonSLVFokkerPlanckFdmParams logParams = {
2622	.xGrid: 201, .vGrid: 401, .tMaxStepsPerYear: 1000, .tMinStepsPerYear: 30, .tStepNumberDecay: 2.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
2623	.x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
2624	.vLowerEps: 1e-5, .vUpperEps: 1e-5, .vMin: 0.0000025, .v0Density: 1.0, .vLowerBoundDensity: 0.1, .vUpperBoundDensity: 0.9, .leverageFctPropEps: 1e-5,
2625	.greensAlgorithm: FdmHestonGreensFct::Gaussian,
2626	.trafoType: FdmSquareRootFwdOp::Log,
2627	.schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
2628	};
2629
2630	const ext::shared_ptr<LocalVolTermStructure> leverageFctFDM =
2631	HestonSLVFDMModel(
2632	localVol, hestonModel, maturityDate, logParams, false, requiredDates,
2633	mixingFactor).leverageFunction();
2634
2635	// Monte-Carlo calibration
2636	const Size timeStepsPerYear = 365;
2637	const Size nBins = 201;
2638	const Size calibrationPaths = 65536;
2639
2640	const ext::shared_ptr<LocalVolTermStructure> leverageFctMC =
2641	HestonSLVMCModel(
2642	localVol, hestonModel,
2643	sobolGeneratorFactory,
2644	maturityDate, timeStepsPerYear, nBins, calibrationPaths, requiredDates,
2645	mixingFactor).leverageFunction();
2646
2647	// Create SLV pricing engines with both leverage functions
2648	const ext::shared_ptr<PricingEngine> fdEngineWithMixingFactor
2649	= ext::make_shared<FdHestonVanillaEngine>(
2650	args: hestonModelPtr, args: 100, args: 100, args: 50, args: 0,
2651	args: FdmSchemeDesc::Hundsdorfer(), args: leverageFctFDM, args&: mixingFactor);
2652
2653	const ext::shared_ptr<PricingEngine> mcEngineWithMixingFactor
2654	= ext::make_shared<FdHestonVanillaEngine>(
2655	args: hestonModelPtr, args: 100, args: 100, args: 50, args: 0,
2656	args: FdmSchemeDesc::Hundsdorfer(), args: leverageFctMC, args&: mixingFactor);
2657
2658	const ext::shared_ptr<PricingEngine> fdBarrierEngineWithMixingFactor
2659	= ext::make_shared<FdHestonBarrierEngine>(
2660	args: hestonModelPtr, args: 100, args: 100, args: 50, args: 0,
2661	args: FdmSchemeDesc::Hundsdorfer(), args: leverageFctFDM, args&: mixingFactor);
2662
2663	const ext::shared_ptr<PricingEngine> mcBarrierEngineWithMixingFactor
2664	= ext::make_shared<FdHestonBarrierEngine>(
2665	args: hestonModelPtr, args: 100, args: 100, args: 50, args: 0,
2666	args: FdmSchemeDesc::Hundsdorfer(), args: leverageFctMC, args&: mixingFactor);
2667
2668	// Price the vanilla and barrier with both engines
2669	vanillaOption.setPricingEngine(fdEngineWithMixingFactor);
2670	const Real priceFDM = vanillaOption.NPV();
2671
2672	vanillaOption.setPricingEngine(mcEngineWithMixingFactor);
2673	const Real priceMC = vanillaOption.NPV();
2674
2675	barrierOption.setPricingEngine(fdBarrierEngineWithMixingFactor);
2676	const Real barrierPriceFDM = barrierOption.NPV();
2677
2678	barrierOption.setPricingEngine(mcBarrierEngineWithMixingFactor);
2679	const Real barrierPriceMC = barrierOption.NPV();
2680
2681	// Check MC and FDM vanilla prices against local vol, and ensure that the barrier
2682	// prices from MC and FDM are also consistent
2683	if (relativeError(x1: priceFDM, x2: localVolPrice, reference: localVolPrice) > epsilon) {
2684	BOOST_ERROR("FDM price does not match with Local Vol"
2685	<< "\n Local Vol Price: " << localVolPrice
2686	<< "\n FDM Price: " << priceFDM
2687	<< "\n Relative Error: " << relativeError(priceFDM, localVolPrice, localVolPrice)
2688	<< "\n Allowed Error: " << epsilon
2689	<< "\n Mixing Factor: " << mixingFactor);
2690	}
2691
2692	if (relativeError(x1: priceMC, x2: localVolPrice, reference: localVolPrice) > epsilon) {
2693	BOOST_ERROR("MC price does not match with Local Vol"
2694	<< "\n Local Vol Price: " << localVolPrice
2695	<< "\n MC Price: " << priceMC
2696	<< "\n Relative Error: " << relativeError(priceMC, localVolPrice, localVolPrice)
2697	<< "\n Allowed Error: " << epsilon
2698	<< "\n Mixing Factor: " << mixingFactor);
2699	}
2700
2701	if (relativeError(x1: barrierPriceFDM, x2: barrierPriceMC, reference: barrierPriceMC) > epsilon) {
2702	BOOST_ERROR("FDM Barrier Price does not match MC Barrier Price"
2703	<< "\n FDM Barrier Price: " << barrierPriceFDM
2704	<< "\n MC Barrier Price: " << barrierPriceMC
2705	<< "\n Relative Error: " << relativeError(barrierPriceFDM, barrierPriceMC, barrierPriceMC)
2706	<< "\n Allowed Error: " << epsilon
2707	<< "\n Mixing Factor: " << mixingFactor);
2708	}
2709	}
2710	}
2711
2712	test_suite* HestonSLVModelTest::suite(SpeedLevel speed) {
2713	auto* suite = BOOST_TEST_SUITE("Heston Stochastic Local Volatility tests");
2714
2715	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testBlackScholesFokkerPlanckFwdEquation));
2716	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testSquareRootZeroFlowBC));
2717	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testTransformedZeroFlowBC));
2718	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testSquareRootEvolveWithStationaryDensity));
2719	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testSquareRootLogEvolveWithStationaryDensity));
2720	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testSquareRootFokkerPlanckFwdEquation));
2721	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testBarrierPricingViaHestonLocalVol));
2722	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testLocalVolsvSLVPropDensity));
2723	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testDiffusionAndDriftSlvProcess));
2724
2725	if (speed <= Fast) {
2726	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testHestonFokkerPlanckFwdEquationLogLVLeverage));
2727	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testMonteCarloVsFdmPricing));
2728	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testBlackScholesFokkerPlanckFwdEquationLocalVol));
2729	}
2730
2731	if (speed == Slow) {
2732	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testHestonFokkerPlanckFwdEquation));
2733	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testMonteCarloCalibration));
2734	suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testMoustacheGraph));
2735	}
2736
2737	// these tests take very long
2738	// suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testForwardSkewSLV));
2739	// suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testFDMCalibration));
2740	// suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testBarrierPricingMixedModels));
2741	// suite->add(QUANTLIB_TEST_CASE(&HestonSLVModelTest::testBarrierPricingMixedModelsMonteCarloVsFdmPricing)); // ~250s
2742
2743	return suite;
2744	}
2745