1	/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3	/*
4	Copyright (C) 2012 Peter Caspers
5
6	This file is part of QuantLib, a free-software/open-source library
7	for financial quantitative analysts and developers - http://quantlib.org/
8
9	QuantLib is free software: you can redistribute it and/or modify it
10	under the terms of the QuantLib license. You should have received a
11	copy of the license along with this program; if not, please email
12	<quantlib-dev@lists.sf.net>. The license is also available online at
13	<http://quantlib.org/license.shtml>.
14
15	This program is distributed in the hope that it will be useful, but WITHOUT
16	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17	FOR A PARTICULAR PURPOSE. See the license for more details.
18	*/
19
20	#include "markovfunctional.hpp"
21	#include "utilities.hpp"
22	#include <ql/processes/mfstateprocess.hpp>
23	#include <ql/models/shortrate/onefactormodels/markovfunctional.hpp>
24	#include <ql/pricingengines/swaption/gaussian1dswaptionengine.hpp>
25	#include <ql/pricingengines/capfloor/gaussian1dcapfloorengine.hpp>
26	#include <ql/termstructures/yield/flatforward.hpp>
27	#include <ql/termstructures/yield/piecewiseyieldcurve.hpp>
28	#include <ql/termstructures/volatility/swaption/swaptionconstantvol.hpp>
29	#include <ql/termstructures/volatility/optionlet/constantoptionletvol.hpp>
30	#include <ql/termstructures/volatility/swaption/swaptionvolmatrix.hpp>
31	#include <ql/termstructures/volatility/swaption/sabrswaptionvolatilitycube.hpp>
32	#include <ql/termstructures/volatility/swaption/interpolatedswaptionvolatilitycube.hpp>
33	#include <ql/termstructures/volatility/capfloor/capfloortermvolsurface.hpp>
34	#include <ql/termstructures/volatility/optionlet/optionletstripper1.hpp>
35	#include <ql/termstructures/volatility/optionlet/strippedoptionletadapter.hpp>
36	#include <ql/termstructures/volatility/interpolatedsmilesection.hpp>
37	#include <ql/termstructures/volatility/kahalesmilesection.hpp>
38	#include <ql/time/calendars/target.hpp>
39	#include <ql/indexes/swap/euriborswap.hpp>
40	#include <ql/indexes/ibor/euribor.hpp>
41	#include <ql/termstructures/yield/ratehelpers.hpp>
42	#include <ql/time/daycounters/actual360.hpp>
43	#include <ql/time/daycounters/thirty360.hpp>
44	#include <ql/time/daycounters/actualactual.hpp>
45	#include <ql/instruments/makeswaption.hpp>
46	#include <ql/instruments/makevanillaswap.hpp>
47	#include <ql/instruments/makecapfloor.hpp>
48	#include <ql/cashflows/cashflowvectors.hpp>
49	#include <ql/pricingengines/swaption/blackswaptionengine.hpp>
50	#include <ql/pricingengines/capfloor/blackcapfloorengine.hpp>
51	#include <ql/models/shortrate/calibrationhelpers/swaptionhelper.hpp>
52	#include <ql/models/shortrate/calibrationhelpers/caphelper.hpp>
53
54	using namespace QuantLib;
55	using namespace boost::unit_test_framework;
56
57	using std::fabs;
58
59	void MarkovFunctionalTest::testMfStateProcess() {
60
61	const Real tolerance = 1E-10;
62	BOOST_TEST_MESSAGE("Testing Markov functional state process...");
63
64	Array times1(0), vols1(1, 1.0);
65	MfStateProcess sp1(0.00, times1, vols1);
66	Real var11 = sp1.variance(t0: 0.0, x0: 0.0, dt: 1.0);
67	Real var12 = sp1.variance(t0: 0.0, x0: 0.0, dt: 2.0);
68	if (std::fabs(x: var11 - 1.0) > tolerance)
69	BOOST_ERROR("process 1 has not variance 1.0 for dt = 1.0 but "
70	<< var11);
71	if (std::fabs(x: var12 - 2.0) > tolerance)
72	BOOST_ERROR("process 1 has not variance 1.0 for dt = 1.0 but "
73	<< var12);
74
75	Array times2(2), vols2(3);
76	times2[0] = 1.0;
77	times2[1] = 2.0;
78	vols2[0] = 1.0;
79	vols2[1] = 2.0;
80	vols2[2] = 3.0;
81	MfStateProcess sp2(0.00, times2, vols2);
82	Real dif21 = sp2.diffusion(t: 0.0, x: 0.0);
83	Real dif22 = sp2.diffusion(t: 0.99, x: 0.0);
84	Real dif23 = sp2.diffusion(t: 1.0, x: 0.0);
85	Real dif24 = sp2.diffusion(t: 1.9, x: 0.0);
86	Real dif25 = sp2.diffusion(t: 2.0, x: 0.0);
87	Real dif26 = sp2.diffusion(t: 3.0, x: 0.0);
88	Real dif27 = sp2.diffusion(t: 5.0, x: 0.0);
89	if (std::fabs(x: dif21 - 1.0) > tolerance)
90	BOOST_ERROR("process 2 has wrong drift at 0.0, should be 1.0 but is "
91	<< dif21);
92	if (std::fabs(x: dif22 - 1.0) > tolerance)
93	BOOST_ERROR("process 2 has wrong drift at 0.99, should be 1.0 but is "
94	<< dif22);
95	if (std::fabs(x: dif23 - 2.0) > tolerance)
96	BOOST_ERROR("process 2 has wrong drift at 1.0, should be 2.0 but is "
97	<< dif23);
98	if (std::fabs(x: dif24 - 2.0) > tolerance)
99	BOOST_ERROR("process 2 has wrong drift at 1.9, should be 2.0 but is "
100	<< dif24);
101	if (std::fabs(x: dif25 - 3.0) > tolerance)
102	BOOST_ERROR("process 2 has wrong drift at 2.0, should be 3.0 but is "
103	<< dif25);
104	if (std::fabs(x: dif26 - 3.0) > tolerance)
105	BOOST_ERROR("process 2 has wrong drift at 3.0, should be 3.0 but is "
106	<< dif26);
107	if (std::fabs(x: dif27 - 3.0) > tolerance)
108	BOOST_ERROR("process 2 has wrong drift at 5.0, should be 3.0 but is "
109	<< dif27);
110	Real var21 = sp2.variance(t0: 0.0, x0: 0.0, dt: 0.0);
111	Real var22 = sp2.variance(t0: 0.0, x0: 0.0, dt: 0.5);
112	Real var23 = sp2.variance(t0: 0.0, x0: 0.0, dt: 1.0);
113	Real var24 = sp2.variance(t0: 0.0, x0: 0.0, dt: 1.5);
114	Real var25 = sp2.variance(t0: 0.0, x0: 0.0, dt: 3.0);
115	Real var26 = sp2.variance(t0: 0.0, x0: 0.0, dt: 5.0);
116	Real var27 = sp2.variance(t0: 1.2, x0: 0.0, dt: 1.0);
117	if (std::fabs(x: var21 - 0.0) > tolerance)
118	BOOST_ERROR("process 2 has wrong variance at 0.0, should be 0.0 but is "
119	<< var21);
120	if (std::fabs(x: var22 - 0.5) > tolerance)
121	BOOST_ERROR("process 2 has wrong variance at 0.5, should be 0.5 but is "
122	<< var22);
123	if (std::fabs(x: var23 - 1.0) > tolerance)
124	BOOST_ERROR("process 2 has wrong variance at 1.0, should be 1.0 but is "
125	<< var23);
126	if (std::fabs(x: var24 - 3.0) > tolerance)
127	BOOST_ERROR("process 2 has wrong variance at 1.5, should be 3.0 but is "
128	<< var24);
129	if (std::fabs(x: var25 - 14.0) > tolerance)
130	BOOST_ERROR(
131	"process 2 has wrong variance at 3.0, should be 14.0 but is "
132	<< var25);
133	if (std::fabs(x: var26 - 32.0) > tolerance)
134	BOOST_ERROR(
135	"process 2 has wrong variance at 5.0, should be 32.0 but is "
136	<< var26);
137	if (std::fabs(x: var27 - 5.0) > tolerance)
138	BOOST_ERROR("process 2 has wrong variance between 1.2 and 2.2, should "
139	"be 5.0 but is "
140	<< var27);
141
142	MfStateProcess sp3(0.01, times2, vols2);
143	Real var31 = sp3.variance(t0: 0.0, x0: 0.0, dt: 0.0);
144	Real var32 = sp3.variance(t0: 0.0, x0: 0.0, dt: 0.5);
145	Real var33 = sp3.variance(t0: 0.0, x0: 0.0, dt: 1.0);
146	Real var34 = sp3.variance(t0: 0.0, x0: 0.0, dt: 1.5);
147	Real var35 = sp3.variance(t0: 0.0, x0: 0.0, dt: 3.0);
148	Real var36 = sp3.variance(t0: 0.0, x0: 0.0, dt: 5.0);
149	Real var37 = sp3.variance(t0: 1.2, x0: 0.0, dt: 1.0);
150	if (std::fabs(x: var31 - 0.0) > tolerance)
151	BOOST_ERROR("process 3 has wrong variance at 0.0, should be 0.0 but is "
152	<< std::setprecision(12) << var31);
153	if (std::fabs(x: var32 - 0.502508354208) > tolerance)
154	BOOST_ERROR("process 3 has wrong variance at 0.5, should be 0.5 but it "
155	<< std::setprecision(12) << var32);
156	if (std::fabs(x: var33 - 1.01006700134) > tolerance)
157	BOOST_ERROR("process 3 has wrong variance at 1.0, should be 1.0 but it "
158	<< std::setprecision(12) << var33);
159	if (std::fabs(x: var34 - 3.06070578669) > tolerance)
160	BOOST_ERROR("process 3 has wrong variance at 1.5, should be 3.0 but it "
161	<< std::setprecision(12) << var34);
162	if (std::fabs(x: var35 - 14.5935513933) > tolerance)
163	BOOST_ERROR(
164	"process 3 has wrong variance at 3.0, should be 14.0 but it "
165	<< std::setprecision(12) << var35);
166	if (std::fabs(x: var36 - 34.0940185819) > tolerance)
167	BOOST_ERROR(
168	"process 3 has wrong variance at 5.0, should be 32.0 but it "
169	<< std::setprecision(12) << var36);
170	if (std::fabs(x: var37 - 5.18130257358) > tolerance)
171	BOOST_ERROR("process 3 has wrong variance between 1.2 and 2.2, should "
172	"be 5.0 but it "
173	<< std::setprecision(12) << var37);
174	}
175
176	namespace {
177
178	// Flat yield term structure at 3%
179
180	Handle<YieldTermStructure> flatYts() {
181
182	return Handle<YieldTermStructure>(ext::shared_ptr<YieldTermStructure>(
183	new FlatForward(0, TARGET(), 0.03, Actual365Fixed())));
184	}
185
186	// Flat swaption volatility structure at 20%
187
188	Handle<SwaptionVolatilityStructure> flatSwaptionVts() {
189
190	return Handle<SwaptionVolatilityStructure>(
191	ext::shared_ptr<SwaptionVolatilityStructure>(
192	new ConstantSwaptionVolatility(0, TARGET(), ModifiedFollowing,
193	0.20, Actual365Fixed())));
194	}
195
196	// Flat cap volatility structure at 20%
197
198	Handle<OptionletVolatilityStructure> flatOptionletVts() {
199
200	return Handle<OptionletVolatilityStructure>(
201	ext::shared_ptr<OptionletVolatilityStructure>(
202	new ConstantOptionletVolatility(0, TARGET(), ModifiedFollowing,
203	0.20, Actual365Fixed())));
204	}
205
206	// Yield term structure as of 14.11.2012 (6m discounting)
207
208	Handle<YieldTermStructure> md0Yts() {
209
210	ext::shared_ptr<IborIndex> euribor6mEmpty(new Euribor(6 * Months));
211
212	std::vector<ext::shared_ptr<Quote> > q6m;
213	std::vector<ext::shared_ptr<RateHelper> > r6m;
214
215	double q6mh[] = { 0.0001, 0.0001, 0.0001, 0.0003, 0.00055, 0.0009,
216	0.0014, 0.0019, 0.0025, 0.0031, 0.00325, 0.00313,
217	0.0031, 0.00307, 0.00309, 0.00339, 0.00316, 0.00326,
218	0.00335, 0.00343, 0.00358, 0.00351, 0.00388, 0.00404,
219	0.00425, 0.00442, 0.00462, 0.00386, 0.00491, 0.00647,
220	0.00837, 0.01033, 0.01218, 0.01382, 0.01527, 0.01654,
221	0.0177, 0.01872, 0.01959, 0.0203, 0.02088, 0.02132,
222	0.02164, 0.02186, 0.02202, 0.02213, 0.02222, 0.02229,
223	0.02234, 0.02238, 0.02241, 0.02243, 0.02244, 0.02245,
224	0.02247, 0.0225, 0.02284, 0.02336, 0.02407, 0.0245 };
225
226	Period q6mh1[] = { 1 * Days, 1 * Days, 1 * Days, 1 * Weeks,
227	1 * Months, 2 * Months, 3 * Months, 4 * Months,
228	5 * Months, 6 * Months };
229
230	Period q6mh2[] = {
231	7 * Months, 8 * Months, 9 * Months, 10 * Months, 11 * Months,
232	1 * Years, 13 * Months, 14 * Months, 15 * Months, 16 * Months,
233	17 * Months, 18 * Months, 19 * Months, 20 * Months, 21 * Months,
234	22 * Months, 23 * Months, 2 * Years, 3 * Years, 4 * Years,
235	5 * Years, 6 * Years, 7 * Years, 8 * Years, 9 * Years,
236	10 * Years, 11 * Years, 12 * Years, 13 * Years, 14 * Years,
237	15 * Years, 16 * Years, 17 * Years, 18 * Years, 19 * Years,
238	20 * Years, 21 * Years, 22 * Years, 23 * Years, 24 * Years,
239	25 * Years, 26 * Years, 27 * Years, 28 * Years, 29 * Years,
240	30 * Years, 35 * Years, 40 * Years, 50 * Years, 60 * Years
241	};
242
243	q6m.reserve(n: 10 + 15 + 35);
244	for (Real i : q6mh) {
245	q6m.push_back(x: ext::shared_ptr<Quote>(new SimpleQuote(i)));
246	}
247
248	r6m.reserve(n: 10);
249	for (int i = 0; i < 10; i++) {
250	r6m.push_back(x: ext::make_shared<DepositRateHelper>(
251	args: Handle<Quote>(q6m[i]), args&: q6mh1[i],
252	args: i < 2 ? i : 2, args: TARGET(),
253	args: ModifiedFollowing, args: false, args: Actual360()));
254	}
255
256	for (int i = 0; i < 18; i++) {
257	if (i + 1 != 6 && i + 1 != 12 && i + 1 != 18) {
258	r6m.push_back(x: ext::make_shared<FraRateHelper>(
259	args: Handle<Quote>(q6m[10 + i]), args: i + 1,
260	args: i + 7, args: 2, args: TARGET(), args: ModifiedFollowing,
261	args: false, args: Actual360()));
262	}
263	}
264
265	for (int i = 0; i < 15 + 35; i++) {
266	if (i + 7 == 12 || i + 7 == 18 || i + 7 >= 24) {
267	r6m.push_back(x: ext::make_shared<SwapRateHelper>(
268	args: Handle<Quote>(q6m[10 + i]), args&: q6mh2[i],
269	args: TARGET(), args: Annual, args: ModifiedFollowing,
270	args: Actual360(), args&: euribor6mEmpty));
271	}
272	}
273
274	Handle<YieldTermStructure> res(ext::shared_ptr<YieldTermStructure>(
275	new PiecewiseYieldCurve<Discount, LogLinear>(0, TARGET(), r6m,
276	Actual365Fixed())));
277	res->enableExtrapolation();
278
279	return res;
280	}
281
282	// Swaption volatility cube as of 14.11.2012, 1y underlying vols are not
283	// converted here from 3m to 6m
284
285	Handle<SwaptionVolatilityStructure> md0SwaptionVts() {
286
287	std::vector<Period> optionTenors = {
288	1 * Months, 2 * Months, 3 * Months, 6 * Months,
289	9 * Months, 1 * Years, 18 * Months, 2 * Years,
290	3 * Years, 4 * Years, 5 * Years, 6 * Years,
291	7 * Years, 8 * Years, 9 * Years, 10 * Years,
292	15 * Years, 20 * Years, 25 * Years, 30 * Years };
293
294	std::vector<Period> swapTenors = {
295	1 * Years, 2 * Years, 3 * Years, 4 * Years,
296	5 * Years, 6 * Years, 7 * Years, 8 * Years,
297	9 * Years, 10 * Years, 15 * Years, 20 * Years,
298	25 * Years, 30 * Years };
299
300	double qSwAtmh[] = {
301	1.81, 0.897, 0.819, 0.692, 0.551, 0.47, 0.416, 0.379, 0.357,
302	0.335, 0.283, 0.279, 0.283, 0.287, 1.717, 0.886, 0.79, 0.69,
303	0.562, 0.481, 0.425, 0.386, 0.359, 0.339, 0.29, 0.287, 0.292,
304	0.296, 1.762, 0.903, 0.804, 0.693, 0.582, 0.5, 0.448, 0.411,
305	0.387, 0.365, 0.31, 0.307, 0.312, 0.317, 1.662, 0.882, 0.764,
306	0.67, 0.586, 0.513, 0.468, 0.432, 0.408, 0.388, 0.331, 0.325,
307	0.33, 0.334, 1.53, 0.854, 0.728, 0.643, 0.565, 0.503, 0.464,
308	0.435, 0.415, 0.393, 0.337, 0.33, 0.333, 0.338, 1.344, 0.786,
309	0.683, 0.609, 0.54, 0.488, 0.453, 0.429, 0.411, 0.39, 0.335,
310	0.329, 0.332, 0.336, 1.1, 0.711, 0.617, 0.548, 0.497, 0.456,
311	0.43, 0.408, 0.392, 0.374, 0.328, 0.323, 0.326, 0.33, 0.956,
312	0.638, 0.553, 0.496, 0.459, 0.427, 0.403, 0.385, 0.371, 0.359,
313	0.321, 0.318, 0.323, 0.327, 0.671, 0.505, 0.45, 0.42, 0.397,
314	0.375, 0.36, 0.347, 0.337, 0.329, 0.305, 0.303, 0.309, 0.313,
315	0.497, 0.406, 0.378, 0.36, 0.348, 0.334, 0.323, 0.315, 0.309,
316	0.304, 0.289, 0.289, 0.294, 0.297, 0.404, 0.352, 0.334, 0.322,
317	0.313, 0.304, 0.296, 0.291, 0.288, 0.286, 0.278, 0.277, 0.281,
318	0.282, 0.345, 0.312, 0.302, 0.294, 0.286, 0.28, 0.276, 0.274,
319	0.273, 0.273, 0.267, 0.265, 0.268, 0.269, 0.305, 0.285, 0.277,
320	0.271, 0.266, 0.262, 0.26, 0.259, 0.26, 0.262, 0.259, 0.256,
321	0.257, 0.256, 0.282, 0.265, 0.259, 0.254, 0.251, 0.25, 0.25,
322	0.251, 0.253, 0.256, 0.253, 0.25, 0.249, 0.246, 0.263, 0.248,
323	0.244, 0.241, 0.24, 0.24, 0.242, 0.245, 0.249, 0.252, 0.249,
324	0.245, 0.243, 0.238, 0.242, 0.234, 0.232, 0.232, 0.233, 0.235,
325	0.239, 0.243, 0.247, 0.249, 0.246, 0.242, 0.237, 0.231, 0.233,
326	0.234, 0.241, 0.246, 0.249, 0.253, 0.257, 0.261, 0.263, 0.264,
327	0.251, 0.236, 0.222, 0.214, 0.262, 0.26, 0.262, 0.263, 0.263,
328	0.266, 0.268, 0.269, 0.269, 0.265, 0.237, 0.214, 0.202, 0.196,
329	0.26, 0.26, 0.261, 0.261, 0.258, 0.255, 0.252, 0.248, 0.245,
330	0.24, 0.207, 0.187, 0.182, 0.176, 0.236, 0.223, 0.221, 0.218,
331	0.214, 0.21, 0.207, 0.204, 0.202, 0.2, 0.175, 0.167, 0.163,
332	0.158
333	};
334
335	std::vector<std::vector<Handle<Quote> > > qSwAtm;
336	for (int i = 0; i < 20; i++) {
337	std::vector<Handle<Quote> > qSwAtmTmp;
338	qSwAtmTmp.reserve(n: 14);
339	for (int j = 0; j < 14; j++) {
340	qSwAtmTmp.emplace_back(
341	args: ext::shared_ptr<Quote>(new SimpleQuote(qSwAtmh[i * 14 + j])));
342	}
343	qSwAtm.push_back(x: qSwAtmTmp);
344	}
345
346	Handle<SwaptionVolatilityStructure> swaptionVolAtm(
347	ext::shared_ptr<SwaptionVolatilityStructure>(
348	new SwaptionVolatilityMatrix(TARGET(), ModifiedFollowing,
349	optionTenors, swapTenors, qSwAtm,
350	Actual365Fixed())));
351
352	std::vector<Period> optionTenorsSmile = { 3 * Months, 1 * Years, 5 * Years,
353	10 * Years, 20 * Years, 30 * Years };
354	std::vector<Period> swapTenorsSmile = { 2 * Years, 5 * Years, 10 * Years,
355	20 * Years, 30 * Years };
356	std::vector<Real> strikeSpreads = { -0.02, -0.01, -0.0050, -0.0025, 0.0,
357	0.0025, 0.0050, 0.01, 0.02 };
358
359	std::vector<std::vector<Handle<Quote> > > qSwSmile;
360
361	double qSwSmileh[] = {
362	2.2562, 2.2562, 2.2562, 0.1851, 0.0, -0.0389, -0.0507,
363	-0.0571, -0.06, 14.9619, 14.9619, 0.1249, 0.0328, 0.0,
364	-0.0075, -0.005, 0.0078, 0.0328, 0.2296, 0.2296, 0.0717,
365	0.0267, 0.0, -0.0115, -0.0126, -0.0002, 0.0345, 0.6665,
366	0.1607, 0.0593, 0.0245, 0.0, -0.0145, -0.0204, -0.0164,
367	0.0102, 0.6509, 0.1649, 0.0632, 0.027, 0.0, -0.018,
368	-0.0278, -0.0303, -0.0105, 0.6303, 0.6303, 0.6303, 0.1169,
369	0.0, -0.0469, -0.0699, -0.091, -0.1065, 0.4437, 0.4437,
370	0.0944, 0.0348, 0.0, -0.0206, -0.0327, -0.0439, -0.0472,
371	2.1557, 0.1501, 0.0531, 0.0225, 0.0, -0.0161, -0.0272,
372	-0.0391, -0.0429, 0.4365, 0.1077, 0.0414, 0.0181, 0.0,
373	-0.0137, -0.0237, -0.0354, -0.0401, 0.4415, 0.1117, 0.0437,
374	0.0193, 0.0, -0.015, -0.0264, -0.0407, -0.0491, 0.4301,
375	0.0776, 0.0283, 0.0122, 0.0, -0.0094, -0.0165, -0.0262,
376	-0.035, 0.2496, 0.0637, 0.0246, 0.0109, 0.0, -0.0086,
377	-0.0153, -0.0247, -0.0334, 0.1912, 0.0569, 0.023, 0.0104,
378	0.0, -0.0085, -0.0155, -0.0256, -0.0361, 0.2095, 0.06,
379	0.0239, 0.0107, 0.0, -0.0087, -0.0156, -0.0254, -0.0348,
380	0.2357, 0.0669, 0.0267, 0.012, 0.0, -0.0097, -0.0174,
381	-0.0282, -0.0383, 0.1291, 0.0397, 0.0158, 0.007, 0.0,
382	-0.0056, -0.01, -0.0158, -0.0203, 0.1281, 0.0397, 0.0159,
383	0.0071, 0.0, -0.0057, -0.0102, -0.0164, -0.0215, 0.1547,
384	0.0468, 0.0189, 0.0085, 0.0, -0.0069, -0.0125, -0.0205,
385	-0.0283, 0.1851, 0.0522, 0.0208, 0.0093, 0.0, -0.0075,
386	-0.0135, -0.0221, -0.0304, 0.1782, 0.0506, 0.02, 0.0089,
387	0.0, -0.0071, -0.0128, -0.0206, -0.0276, 0.2665, 0.0654,
388	0.0255, 0.0113, 0.0, -0.0091, -0.0163, -0.0265, -0.0367,
389	0.2873, 0.0686, 0.0269, 0.0121, 0.0, -0.0098, -0.0179,
390	-0.0298, -0.043, 0.2993, 0.0688, 0.0273, 0.0123, 0.0,
391	-0.0103, -0.0189, -0.0324, -0.0494, 0.1869, 0.0501, 0.0202,
392	0.0091, 0.0, -0.0076, -0.014, -0.0239, -0.0358, 0.1573,
393	0.0441, 0.0178, 0.008, 0.0, -0.0066, -0.0121, -0.0202,
394	-0.0294, 0.196, 0.0525, 0.0204, 0.009, 0.0, -0.0071,
395	-0.0125, -0.0197, -0.0253, 0.1795, 0.0497, 0.0197, 0.0088,
396	0.0, -0.0071, -0.0128, -0.0208, -0.0286, 0.1401, 0.0415,
397	0.0171, 0.0078, 0.0, -0.0066, -0.0122, -0.0209, -0.0318,
398	0.112, 0.0344, 0.0142, 0.0065, 0.0, -0.0055, -0.01,
399	-0.0171, -0.0256, 0.1077, 0.0328, 0.0134, 0.0061, 0.0,
400	-0.005, -0.0091, -0.0152, -0.0216,
401	};
402
403	for (int i = 0; i < 30; i++) {
404	std::vector<Handle<Quote> > qSwSmileTmp;
405	qSwSmileTmp.reserve(n: 9);
406	for (int j = 0; j < 9; j++) {
407	qSwSmileTmp.emplace_back(
408	args: ext::shared_ptr<Quote>(new SimpleQuote(qSwSmileh[i * 9 + j])));
409	}
410	qSwSmile.push_back(x: qSwSmileTmp);
411	}
412
413	double qSwSmileh1[] = {
414	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
415	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
416	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
417	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
418	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
419	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
420	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
421	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
422	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2,
423	0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2, 0.01, 0.2, 0.8, -0.2
424	};
425
426	std::vector<bool> parameterFixed = {false, false, false, false};
427
428	std::vector<std::vector<Handle<Quote> > > parameterGuess;
429	for (int i = 0; i < 30; i++) {
430	std::vector<Handle<Quote> > parameterGuessTmp;
431	parameterGuessTmp.reserve(n: 4);
432	for (int j = 0; j < 4; j++) {
433	parameterGuessTmp.emplace_back(
434	args: ext::shared_ptr<Quote>(new SimpleQuote(qSwSmileh1[i * 4 + j])));
435	}
436	parameterGuess.push_back(x: parameterGuessTmp);
437	}
438
439	ext::shared_ptr<EndCriteria> ec(
440	new EndCriteria(50000, 250, 1E-6, 1E-6, 1E-6));
441
442	ext::shared_ptr<SwapIndex> swapIndex(new EuriborSwapIsdaFixA(
443	30 * Years, Handle<YieldTermStructure>(md0Yts())));
444	ext::shared_ptr<SwapIndex> shortSwapIndex(new EuriborSwapIsdaFixA(
445	1 * Years,
446	Handle<YieldTermStructure>(md0Yts()))); // We assume that we have 6m
447	// vols (which we actually
448	// don't have for 1y
449	// underlying, but this is
450	// just a test...)
451
452	// return Handle<SwaptionVolatilityStructure>(new
453	// InterpolatedSwaptionVolatilityCube(swaptionVolAtm,optionTenorsSmile,swapTenorsSmile,strikeSpreads,qSwSmile,swapIndex,shortSwapIndex,false));
454	// // bilinear interpolation gives nasty digitals
455	Handle<SwaptionVolatilityStructure> res(
456	ext::shared_ptr<SwaptionVolatilityStructure>(new SabrSwaptionVolatilityCube(
457	swaptionVolAtm, optionTenorsSmile, swapTenorsSmile,
458	strikeSpreads, qSwSmile, swapIndex, shortSwapIndex, true,
459	parameterGuess, parameterFixed, true, ec,
460	0.0050))); // put a big error tolerance here ... we just want a
461	// smooth cube for testing
462	res->enableExtrapolation();
463	return res;
464	}
465
466	// Cap volatility surface as of 14.11.2012. Par vols up to 2y are converted
467	// to 6m to get a consistent caplet surface.
468
469	Handle<OptionletVolatilityStructure> md0OptionletVts() {
470
471	// with the thread safe observer it takes very long to destruct
472	// the cap floor instruments created in OptionletStripper1
473	#ifdef QL_ENABLE_THREAD_SAFE_OBSERVER_PATTERN
474	return flatOptionletVts();
475	#endif
476
477	Size nOptTen = 16;
478	Size nStrikes = 12; // leave out last strike 10% because it causes an
479	// exception in bootstrapping
480
481	std::vector<Period> optionTenors = {
482	1 * Years, 18 * Months, 2 * Years, 3 * Years,
483	4 * Years, 5 * Years, 6 * Years, 7 * Years,
484	8 * Years, 9 * Years, 10 * Years, 12 * Years,
485	15 * Years, 20 * Years, 25 * Years, 30 * Years };
486
487	std::vector<Real> strikes = {
488	0.0025, 0.0050, 0.0100, 0.0150, 0.0200, 0.0225,
489	0.0250, 0.0300, 0.0350, 0.0400, 0.0500, 0.0600 };
490
491	Matrix vols(nOptTen, nStrikes);
492	Real volsa[13][16] = { { 1.3378, 1.3032, 1.2514, 1.081, 1.019, 0.961,
493	0.907, 0.862, 0.822, 0.788, 0.758, 0.709,
494	0.66, 0.619, 0.597, 0.579 }, // strike1
495	{ 1.1882, 1.1057, 0.9823, 0.879, 0.828, 0.779,
496	0.736, 0.7, 0.67, 0.644, 0.621, 0.582,
497	0.544, 0.513, 0.496, 0.482 }, // strike2
498	{ 1.1646, 1.0356, 0.857, 0.742, 0.682, 0.626,
499	0.585, 0.553, 0.527, 0.506, 0.488, 0.459,
500	0.43, 0.408, 0.396, 0.386 }, // ...
501	{ 1.1932, 1.0364, 0.8291, 0.691, 0.618, 0.553,
502	0.509, 0.477, 0.452, 0.433, 0.417, 0.391,
503	0.367, 0.35, 0.342, 0.335 },
504	{ 1.2233, 1.0489, 0.8268, 0.666, 0.582, 0.51,
505	0.463, 0.43, 0.405, 0.387, 0.372, 0.348,
506	0.326, 0.312, 0.306, 0.301 },
507	{ 1.2369, 1.0555, 0.8283, 0.659, 0.57, 0.495,
508	0.447, 0.414, 0.388, 0.37, 0.355, 0.331,
509	0.31, 0.298, 0.293, 0.289 },
510	{ 1.2498, 1.0622, 0.8307, 0.653, 0.56, 0.483,
511	0.434, 0.4, 0.374, 0.356, 0.341, 0.318,
512	0.297, 0.286, 0.282, 0.279 },
513	{ 1.2719, 1.0747, 0.8368, 0.646, 0.546, 0.465,
514	0.415, 0.38, 0.353, 0.335, 0.32, 0.296,
515	0.277, 0.268, 0.265, 0.263 },
516	{ 1.2905, 1.0858, 0.8438, 0.643, 0.536, 0.453,
517	0.403, 0.367, 0.339, 0.32, 0.305, 0.281,
518	0.262, 0.255, 0.254, 0.252 },
519	{ 1.3063, 1.0953, 0.8508, 0.642, 0.53, 0.445,
520	0.395, 0.358, 0.329, 0.31, 0.294, 0.271,
521	0.252, 0.246, 0.246, 0.244 },
522	{ 1.332, 1.1108, 0.8631, 0.642, 0.521, 0.436,
523	0.386, 0.348, 0.319, 0.298, 0.282, 0.258,
524	0.24, 0.237, 0.237, 0.236 },
525	{ 1.3513, 1.1226, 0.8732, 0.645, 0.517, 0.43,
526	0.381, 0.344, 0.314, 0.293, 0.277, 0.252,
527	0.235, 0.233, 0.234, 0.233 },
528	{ 1.395, 1.1491, 0.9003, 0.661, 0.511, 0.425,
529	0.38, 0.344, 0.314, 0.292, 0.275, 0.251,
530	0.236, 0.236, 0.238, 0.235 } };
531
532	for (Size i = 0; i < nStrikes; i++) {
533	for (Size j = 0; j < nOptTen; j++) {
534	vols[j][i] = volsa[i][j];
535	}
536	}
537
538	ext::shared_ptr<IborIndex> iborIndex(
539	new Euribor(6 * Months, md0Yts()));
540	ext::shared_ptr<CapFloorTermVolSurface> cf(new CapFloorTermVolSurface(
541	0, TARGET(), ModifiedFollowing, optionTenors, strikes, vols));
542	ext::shared_ptr<OptionletStripper> stripper(
543	new OptionletStripper1(cf, iborIndex));
544
545	return Handle<OptionletVolatilityStructure>(
546	ext::shared_ptr<OptionletVolatilityStructure>(
547	new StrippedOptionletAdapter(stripper)));
548	}
549
550	// Calibration Basket 1: CMS10y Swaptions, 5 yearly fixings
551
552	std::vector<Date> expiriesCalBasket1() {
553
554	std::vector<Date> res;
555	Date referenceDate_ = Settings::instance().evaluationDate();
556
557	for (int i = 1; i <= 5; i++)
558	res.push_back(x: TARGET().advance(date: referenceDate_, period: i * Years));
559
560	return res;
561	}
562
563	std::vector<Period> tenorsCalBasket1() {
564
565	std::vector<Period> res(5, 10 * Years);
566
567	return res;
568	}
569
570	// Calibration Basket 2: 6m caplets, 5 years
571
572	std::vector<Date> expiriesCalBasket2() {
573
574	Date referenceDate_ = Settings::instance().evaluationDate();
575
576	std::vector<Date> res = {
577	TARGET().advance(date: referenceDate_, period: 6 * Months),
578	TARGET().advance(date: referenceDate_, period: 12 * Months),
579	TARGET().advance(date: referenceDate_, period: 18 * Months),
580	TARGET().advance(date: referenceDate_, period: 24 * Months),
581	TARGET().advance(date: referenceDate_, period: 30 * Months),
582	TARGET().advance(date: referenceDate_, period: 36 * Months),
583	TARGET().advance(date: referenceDate_, period: 42 * Months),
584	TARGET().advance(date: referenceDate_, period: 48 * Months),
585	TARGET().advance(date: referenceDate_, period: 54 * Months),
586	TARGET().advance(date: referenceDate_, period: 60 * Months)
587	};
588
589	return res;
590	}
591
592	// Calibration Basket 3: Coterminal Swaptions 10y
593
594	std::vector<Date> expiriesCalBasket3() {
595
596	Date referenceDate_ = Settings::instance().evaluationDate();
597
598	std::vector<Date> res = {
599	TARGET().advance(date: referenceDate_, period: 1 * Years),
600	TARGET().advance(date: referenceDate_, period: 2 * Years),
601	TARGET().advance(date: referenceDate_, period: 3 * Years),
602	TARGET().advance(date: referenceDate_, period: 4 * Years),
603	TARGET().advance(date: referenceDate_, period: 5 * Years),
604	TARGET().advance(date: referenceDate_, period: 6 * Years),
605	TARGET().advance(date: referenceDate_, period: 7 * Years),
606	TARGET().advance(date: referenceDate_, period: 8 * Years),
607	TARGET().advance(date: referenceDate_, period: 9 * Years)
608	};
609
610	return res;
611	}
612
613	std::vector<Period> tenorsCalBasket3() {
614
615	std::vector<Period> res = {9 * Years, 8 * Years, 7 * Years, 6 * Years, 5 * Years,
616	4 * Years, 3 * Years, 2 * Years, 1 * Years};
617	return res;
618	}
619
620	std::vector<Real> impliedStdDevs(const Real atm, const std::vector<Real> &strikes,
621	const std::vector<Real> &prices) {
622
623	std::vector<Real> result;
624
625	for (Size i = 0; i < prices.size(); i++) {
626	result.push_back(x: blackFormulaImpliedStdDev(optionType: Option::Call, strike: strikes[i],
627	forward: atm, blackPrice: prices[i], discount: 1.0, displacement: 0.0,
628	guess: 0.2, accuracy: 1E-8, maxIterations: 1000));
629	}
630
631	return result;
632	}
633	}
634
635	void MarkovFunctionalTest::testKahaleSmileSection() {
636
637	BOOST_TEST_MESSAGE("Testing Kahale smile section...");
638
639	const Real tol = 1E-8;
640
641	// arbitrage free sample smile data
642
643	const Real atm = 0.05;
644	const Real t = 1.0;
645
646	const Real strikes0[] = { 0.01, 0.02, 0.03, 0.04, 0.05,
647	0.06, 0.07, 0.08, 0.09, 0.10 };
648
649	const std::vector<Real> strikes(strikes0, strikes0 + 10);
650	std::vector<Real> money;
651	std::vector<Real> calls0;
652
653	for (Real strike : strikes) {
654	money.push_back(x: strike / atm);
655	calls0.push_back(x: blackFormula(optionType: Option::Call, strike, forward: atm, stdDev: 0.50 * std::sqrt(x: t), discount: 1.0, displacement: 0.0));
656	}
657
658	std::vector<Real> stdDevs0 = impliedStdDevs(atm, strikes, prices: calls0);
659	ext::shared_ptr<SmileSection> sec1(
660	new InterpolatedSmileSection<Linear>(t, strikes, stdDevs0, atm));
661
662	// test arbitrage free smile reproduction
663
664	ext::shared_ptr<KahaleSmileSection> ksec11(
665	new KahaleSmileSection(sec1, atm, false, false, false, money));
666
667	if (std::fabs(x: ksec11->leftCoreStrike() - 0.01) > tol)
668	BOOST_ERROR("smile11 left af strike is " << ksec11->leftCoreStrike()
669	<< " expected 0.01");
670
671	if (std::fabs(x: ksec11->rightCoreStrike() - 0.10) > tol)
672	BOOST_ERROR("smile11 right af strike is " << ksec11->rightCoreStrike()
673	<< " expected 0.10");
674
675	Real k = strikes[0];
676	while (k <= strikes.back() + tol) {
677	Real pric0 = sec1->optionPrice(strike: k);
678	Real pric1 = ksec11->optionPrice(strike: k);
679	if (std::fabs(x: pric0 - pric1) > tol)
680	BOOST_ERROR("smile11 is not reprocduced at strike "
681	<< k << " input smile call price is " << pric0
682	<< " kahale smile call price is " << pric1);
683	k += 0.0001;
684	}
685
686	// test interpolation
687
688	ext::shared_ptr<KahaleSmileSection> ksec12(
689	new KahaleSmileSection(sec1, atm, true, false, false, money));
690
691	// sanity check for left point extrapolation may mark 0.01 as bad as well as
692	// good depending
693	// on platform and compiler due to numerical differences, so we have to
694	// admit two possible results
695	if (std::fabs(x: ksec12->leftCoreStrike() - 0.02) > tol &&
696	std::fabs(x: ksec12->leftCoreStrike() - 0.01) > tol)
697	BOOST_ERROR("smile12 left af strike is " << ksec12->leftCoreStrike()
698	<< "expected 0.01 or 0.02");
699
700	if (std::fabs(x: ksec12->rightCoreStrike() - 0.10) > tol)
701	BOOST_ERROR("smile12 right af strike is " << ksec12->rightCoreStrike()
702	<< "expected 0.10");
703
704	for (Size i = 1; i < strikes.size(); i++) {
705	Real pric0 = sec1->optionPrice(strike: strikes[i]);
706	Real pric1 = ksec12->optionPrice(strike: strikes[i]);
707	if (std::fabs(x: pric0 - pric1) > tol)
708	BOOST_ERROR("smile12 is not reproduced on grid at strike "
709	<< strikes[i] << " input smile call price is " << pric0
710	<< " kahale smile call price is " << pric1);
711	}
712
713	// test global no arbitrageability
714
715	k = 0.0010;
716	Real dig00 = 1.0, dig10 = 1.0;
717	while (k <= 2.0 * strikes.back() + tol) {
718	Real dig0 = ksec11->digitalOptionPrice(strike: k);
719	Real dig1 = ksec12->digitalOptionPrice(strike: k);
720	if (!(dig0 <= dig00 + tol && dig0 >= 0.0))
721	BOOST_ERROR("arbitrage in digitals11 (" << dig00 << "," << dig0
722	<< ") at strike " << k);
723	if (!(dig1 <= dig10 + tol && dig1 >= 0.0))
724	BOOST_ERROR("arbitrage in digitals12 (" << dig10 << "," << dig1
725	<< ") at strike " << k);
726	dig00 = dig0;
727	dig10 = dig1;
728	k += 0.0001;
729	}
730
731	// test exponential extrapolation
732
733	ext::shared_ptr<KahaleSmileSection> ksec13(
734	new KahaleSmileSection(sec1, atm, false, true, false, money));
735
736	k = strikes.back();
737	Real dig0 = ksec13->digitalOptionPrice(strike: k - 0.0010);
738	while (k <= 10.0 * strikes.back() + tol) {
739	Real dig = ksec13->digitalOptionPrice(strike: k);
740	if (!(dig <= dig0 + tol && dig >= 0.0))
741	BOOST_ERROR("arbitrage in digitals13 (" << dig0 << "," << dig
742	<< ") at strike " << k);
743	k += 0.0001;
744	}
745
746	// test arbitrageable smile (leftmost point)
747
748	std::vector<Real> calls1(calls0);
749	calls1[0] = (atm - strikes[0]) +
750	0.0010; // introduce arbitrage by changing call price
751	std::vector<Real> stdDevs1 = impliedStdDevs(atm, strikes, prices: calls1);
752	ext::shared_ptr<SmileSection> sec2(
753	new InterpolatedSmileSection<Linear>(t, strikes, stdDevs1, atm));
754
755	ext::shared_ptr<KahaleSmileSection> ksec21(
756	new KahaleSmileSection(sec2, atm, false, false, false, money));
757	ext::shared_ptr<KahaleSmileSection> ksec22(
758	new KahaleSmileSection(sec2, atm, true, false, true, money));
759
760	if (std::fabs(x: ksec21->leftCoreStrike() - 0.02) > tol)
761	BOOST_ERROR("smile21 left af strike is " << ksec21->leftCoreStrike()
762	<< " expected 0.02");
763	if (std::fabs(x: ksec22->leftCoreStrike() - 0.02) > tol)
764	BOOST_ERROR("smile22 left af strike is " << ksec22->leftCoreStrike()
765	<< " expected 0.02");
766
767	if (std::fabs(x: ksec21->rightCoreStrike() - 0.10) > tol)
768	BOOST_ERROR("smile21 right af strike is " << ksec21->rightCoreStrike()
769	<< " expected 0.10");
770	if (std::fabs(x: ksec22->rightCoreStrike() - 0.10) > tol)
771	BOOST_ERROR("smile22 right af strike is " << ksec22->rightCoreStrike()
772	<< " expected 0.10");
773
774	k = 0.0010;
775	dig00 = dig10 = 1.0;
776	while (k <= 2.0 * strikes.back() + tol) {
777	Real dig0 = ksec21->digitalOptionPrice(strike: k);
778	Real dig1 = ksec22->digitalOptionPrice(strike: k);
779	if (!(dig0 <= dig00 + tol && dig0 >= 0.0))
780	BOOST_ERROR("arbitrage in digitals21 (" << dig00 << "," << dig0
781	<< ") at strike " << k);
782	if (!(dig1 <= dig10 + tol && dig1 >= 0.0))
783	BOOST_ERROR("arbitrage in digitals22 (" << dig10 << "," << dig1
784	<< ") at strike " << k);
785	dig00 = dig0;
786	dig10 = dig1;
787	k += 0.0001;
788	}
789
790	// test arbitrageable smile (second but rightmost point)
791
792	std::vector<Real> calls2(calls0);
793	calls2[8] = 0.9 * calls2[9] +
794	0.1 * calls2[8]; // introduce arbitrage by changing call price
795	std::vector<Real> stdDevs2 = impliedStdDevs(atm, strikes, prices: calls2);
796	ext::shared_ptr<SmileSection> sec3(
797	new InterpolatedSmileSection<Linear>(t, strikes, stdDevs2, atm));
798
799	ext::shared_ptr<KahaleSmileSection> ksec31(
800	new KahaleSmileSection(sec3, atm, false, false, false, money));
801	ext::shared_ptr<KahaleSmileSection> ksec32(
802	new KahaleSmileSection(sec3, atm, true, false, true, money));
803
804	if (std::fabs(x: ksec31->leftCoreStrike() - 0.01) > tol)
805	BOOST_ERROR("smile31 left af strike is " << ksec31->leftCoreStrike()
806	<< " expected 0.01");
807
808	// sanity check for left point extrapolation may mark 0.01 as bad as well as
809	// good depending
810	// on platform and compiler due to numerical differences, so we have to
811	// admit two possible results
812	if (std::fabs(x: ksec32->leftCoreStrike() - 0.02) > tol &&
813	std::fabs(x: ksec32->leftCoreStrike() - 0.01) > tol)
814	BOOST_ERROR("smile32 left af strike is " << ksec32->leftCoreStrike()
815	<< " expected 0.01 or 0.02");
816
817	if (std::fabs(x: ksec31->rightCoreStrike() - 0.08) > tol)
818	BOOST_ERROR("smile31 right af strike is " << ksec31->rightCoreStrike()
819	<< " expected 0.08");
820	if (std::fabs(x: ksec32->rightCoreStrike() - 0.10) > tol)
821	BOOST_ERROR("smile32 right af strike is " << ksec32->rightCoreStrike()
822	<< " expected 0.10");
823	k = 0.0010;
824	dig00 = dig10 = 1.0;
825	while (k <= 2.0 * strikes.back() + tol) {
826	Real dig0 = ksec31->digitalOptionPrice(strike: k);
827	Real dig1 = ksec32->digitalOptionPrice(strike: k);
828	if (!(dig0 <= dig00 + tol && dig0 >= 0.0))
829	BOOST_ERROR("arbitrage in digitals31 (" << dig00 << "," << dig0
830	<< ") at strike " << k);
831	if (!(dig1 <= dig10 + tol && dig1 >= 0.0))
832	BOOST_ERROR("arbitrage in digitals32 (" << dig10 << "," << dig1
833	<< ") at strike " << k);
834	dig00 = dig0;
835	dig10 = dig1;
836	k += 0.0001;
837	}
838	}
839
840	void MarkovFunctionalTest::testCalibrationOneInstrumentSet() {
841
842	const Real tol0 = 0.0001; // 1bp tolerance for model zero rates vs. market
843	// zero rates (note that model zero rates are
844	// implied by the calibration of the numeraire to
845	// the smile)
846	const Real tol1 =
847	0.0001; // 1bp tolerance for model call put premia vs. market premia
848
849	BOOST_TEST_MESSAGE(
850	"Testing Markov functional calibration to one instrument set...");
851
852	Date savedEvalDate = Settings::instance().evaluationDate();
853	Date referenceDate(14, November, 2012);
854	Settings::instance().evaluationDate() = referenceDate;
855
856	Handle<YieldTermStructure> flatYts_ = flatYts();
857	Handle<YieldTermStructure> md0Yts_ = md0Yts();
858	Handle<SwaptionVolatilityStructure> flatSwaptionVts_ = flatSwaptionVts();
859	Handle<SwaptionVolatilityStructure> md0SwaptionVts_ = md0SwaptionVts();
860	Handle<OptionletVolatilityStructure> flatOptionletVts_ = flatOptionletVts();
861	Handle<OptionletVolatilityStructure> md0OptionletVts_ = md0OptionletVts();
862
863	ext::shared_ptr<SwapIndex> swapIndexBase(
864	new EuriborSwapIsdaFixA(1 * Years));
865	ext::shared_ptr<IborIndex> iborIndex(new Euribor(6 * Months));
866
867	std::vector<Date> volStepDates;
868	std::vector<Real> vols = {1.0};
869
870	// use a grid with fewer points for smile arbitrage
871	// testing and model outputs than the default grid
872	// for the testing here
873	std::vector<Real> money = {0.1, 0.25, 0.50, 0.75, 1.0, 1.25, 1.50, 2.0, 5.0};
874
875	// Calibration Basket 1 / flat yts, vts
876
877	ext::shared_ptr<MarkovFunctional> mf1(new MarkovFunctional(
878	flatYts_, 0.01, volStepDates, vols, flatSwaptionVts_,
879	expiriesCalBasket1(), tenorsCalBasket1(), swapIndexBase,
880	MarkovFunctional::ModelSettings()
881	.withYGridPoints(n: 64) // we use the default values more or less, this
882	// is just to demonstrate how to set the model
883	// parameters
884	.withYStdDevs(s: 7.0)
885	.withGaussHermitePoints(n: 32)
886	.withDigitalGap(d: 1e-5)
887	.withMarketRateAccuracy(a: 1e-7)
888	.withLowerRateBound(l: 0.0)
889	.withUpperRateBound(u: 2.0)
890	.withAdjustments(
891	a: MarkovFunctional::ModelSettings::KahaleSmile |
892	MarkovFunctional::ModelSettings::SmileExponentialExtrapolation)
893	.withSmileMoneynessCheckpoints(m: money)));
894
895	MarkovFunctional::ModelOutputs outputs1 =
896	mf1->modelOutputs(); // this costs a lot of time, so only use it if you
897	// want to check the calibration
898	// BOOST_TEST_MESSAGE(outputs1);
899
900	for (Size i = 0; i < outputs1.expiries_.size(); i++) {
901	if (fabs(x: outputs1.marketZerorate_[i] - outputs1.modelZerorate_[i]) >
902	tol0)
903	BOOST_ERROR("Basket 1 / flat termstructures : Market zero rate ("
904	<< outputs1.marketZerorate_[i]
905	<< ") and model zero rate ("
906	<< outputs1.modelZerorate_[i] << ") do not agree.");
907	}
908
909	for (Size i = 0; i < outputs1.expiries_.size(); i++) {
910	for (Size j = 0; j < outputs1.smileStrikes_[i].size(); j++) {
911	if (fabs(x: outputs1.marketCallPremium_[i][j] -
912	outputs1.modelCallPremium_[i][j]) > tol1)
913	BOOST_ERROR(
914	"Basket 1 / flat termstructures : Market call premium ("
915	<< outputs1.marketCallPremium_[i][j]
916	<< ") does not match model premium ("
917	<< outputs1.modelCallPremium_[i][j] << ")");
918	if (fabs(x: outputs1.marketPutPremium_[i][j] -
919	outputs1.modelPutPremium_[i][j]) > tol1)
920	BOOST_ERROR(
921	"Basket 1 / flat termstructures : Market put premium ("
922	<< outputs1.marketPutPremium_[i][j]
923	<< ") does not match model premium ("
924	<< outputs1.modelPutPremium_[i][j] << ")");
925	}
926	}
927
928	// Calibration Basket 2 / flat yts, vts
929
930	ext::shared_ptr<MarkovFunctional> mf2(new MarkovFunctional(
931	flatYts_, 0.01, volStepDates, vols, flatOptionletVts_,
932	expiriesCalBasket2(), iborIndex,
933	MarkovFunctional::ModelSettings()
934	.withYGridPoints(n: 64)
935	.withYStdDevs(s: 7.0)
936	.withGaussHermitePoints(n: 32)
937	.withDigitalGap(d: 1e-5)
938	.withMarketRateAccuracy(a: 1e-7)
939	.withLowerRateBound(l: 0.0)
940	.withUpperRateBound(u: 2.0)
941	.withAdjustments(a: MarkovFunctional::ModelSettings::AdjustNone)
942	.withSmileMoneynessCheckpoints(m: money)));
943
944	MarkovFunctional::ModelOutputs outputs2 = mf2->modelOutputs();
945	// BOOST_TEST_MESSAGE(outputs2);
946
947	for (Size i = 0; i < outputs2.expiries_.size(); i++) {
948	if (fabs(x: outputs2.marketZerorate_[i] - outputs2.modelZerorate_[i]) >
949	tol0)
950	BOOST_ERROR("Basket 2 / flat termstructures : Market zero rate ("
951	<< outputs2.marketZerorate_[i]
952	<< ") and model zero rate ("
953	<< outputs2.modelZerorate_[i] << ") do not agree.");
954	}
955
956	for (Size i = 0; i < outputs2.expiries_.size(); i++) {
957	for (Size j = 0; j < outputs2.smileStrikes_[i].size(); j++) {
958	if (fabs(x: outputs2.marketCallPremium_[i][j] -
959	outputs2.modelCallPremium_[i][j]) > tol1)
960	BOOST_ERROR(
961	"Basket 2 / flat termstructures : Market call premium ("
962	<< outputs2.marketCallPremium_[i][j]
963	<< ") does not match model premium ("
964	<< outputs2.modelCallPremium_[i][j] << ")");
965	if (fabs(x: outputs2.marketPutPremium_[i][j] -
966	outputs2.modelPutPremium_[i][j]) > tol1)
967	BOOST_ERROR(
968	"Basket 2/ flat termstructures : Market put premium ("
969	<< outputs2.marketPutPremium_[i][j]
970	<< ") does not match model premium ("
971	<< outputs2.modelPutPremium_[i][j] << ")");
972	}
973	}
974
975	// Calibration Basket 1 / real yts, vts
976
977	ext::shared_ptr<MarkovFunctional> mf3(new MarkovFunctional(
978	md0Yts_, 0.01, volStepDates, vols, md0SwaptionVts_,
979	expiriesCalBasket1(), tenorsCalBasket1(), swapIndexBase,
980	MarkovFunctional::ModelSettings()
981	.withYGridPoints(n: 128) // use more points to increase accuracy
982	.withYStdDevs(s: 7.0)
983	.withGaussHermitePoints(n: 64)
984	.withDigitalGap(d: 1e-5)
985	.withMarketRateAccuracy(a: 1e-7)
986	.withLowerRateBound(l: 0.0)
987	.withUpperRateBound(u: 2.0)
988	.withSmileMoneynessCheckpoints(m: money)));
989
990	MarkovFunctional::ModelOutputs outputs3 = mf3->modelOutputs();
991	// BOOST_TEST_MESSAGE(outputs3);
992	// outputSurfaces(mf3,md0Yts_);
993
994	for (Size i = 0; i < outputs3.expiries_.size(); i++) {
995	if (fabs(x: outputs3.marketZerorate_[i] - outputs3.modelZerorate_[i]) >
996	tol0)
997	BOOST_ERROR("Basket 1 / real termstructures: Market zero rate ("
998	<< outputs3.marketZerorate_[i]
999	<< ") and model zero rate ("
1000	<< outputs3.modelZerorate_[i] << ") do not agree.");
1001	}
1002
1003	for (Size i = 0; i < outputs3.expiries_.size(); i++) {
1004	for (Size j = 0; j < outputs3.smileStrikes_[i].size(); j++) {
1005	if (fabs(x: outputs3.marketCallPremium_[i][j] -
1006	outputs3.modelCallPremium_[i][j]) > tol1)
1007	BOOST_ERROR(
1008	"Basket 1 / real termstructures: Market call premium ("
1009	<< outputs3.marketCallPremium_[i][j]
1010	<< ") does not match model premium ("
1011	<< outputs3.modelCallPremium_[i][j] << ")");
1012	if (fabs(x: outputs3.marketPutPremium_[i][j] -
1013	outputs3.modelPutPremium_[i][j]) > tol1)
1014	BOOST_ERROR(
1015	"Basket 1 / real termstructures: Market put premium ("
1016	<< outputs3.marketPutPremium_[i][j]
1017	<< ") does not match model premium ("
1018	<< outputs3.modelPutPremium_[i][j] << ")");
1019	}
1020	}
1021
1022	// Calibration Basket 2 / real yts, vts
1023
1024	ext::shared_ptr<MarkovFunctional> mf4(
1025	new MarkovFunctional(md0Yts_, 0.01, volStepDates, vols,
1026	md0OptionletVts_, expiriesCalBasket2(), iborIndex,
1027	MarkovFunctional::ModelSettings()
1028	.withYGridPoints(n: 64)
1029	.withYStdDevs(s: 7.0)
1030	.withGaussHermitePoints(n: 32)
1031	.withDigitalGap(d: 1e-5)
1032	.withMarketRateAccuracy(a: 1e-7)
1033	.withLowerRateBound(l: 0.0)
1034	.withUpperRateBound(u: 2.0)
1035	.withSmileMoneynessCheckpoints(m: money)));
1036
1037	MarkovFunctional::ModelOutputs outputs4 = mf4->modelOutputs();
1038	// BOOST_TEST_MESSAGE(outputs4);
1039
1040	for (Size i = 0; i < outputs4.expiries_.size(); i++) {
1041	if (fabs(x: outputs4.marketZerorate_[i] - outputs4.modelZerorate_[i]) >
1042	tol0)
1043	BOOST_ERROR("Basket 2 / real termstructures : Market zero rate ("
1044	<< outputs4.marketZerorate_[i]
1045	<< ") and model zero rate ("
1046	<< outputs4.modelZerorate_[i] << ") do not agree.");
1047	}
1048
1049	for (Size i = 0; i < outputs4.expiries_.size(); i++) {
1050	for (Size j = 0; j < outputs4.smileStrikes_[i].size(); j++) {
1051	if (fabs(x: outputs4.marketCallPremium_[i][j] -
1052	outputs4.modelCallPremium_[i][j]) > tol1)
1053	BOOST_ERROR(
1054	"Basket 2 / real termstructures : Market call premium ("
1055	<< outputs4.marketCallPremium_[i][j]
1056	<< ") does not match model premium ("
1057	<< outputs4.modelCallPremium_[i][j] << ")");
1058	if (fabs(x: outputs4.marketPutPremium_[i][j] -
1059	outputs4.modelPutPremium_[i][j]) > tol1)
1060	BOOST_ERROR(
1061	"Basket 2/ real termstructures : Market put premium ("
1062	<< outputs4.marketPutPremium_[i][j]
1063	<< ") does not match model premium ("
1064	<< outputs4.modelPutPremium_[i][j] << ")");
1065	}
1066	}
1067
1068	Settings::instance().evaluationDate() = savedEvalDate;
1069	}
1070
1071	void MarkovFunctionalTest::testVanillaEngines() {
1072
1073	const Real tol1 = 0.0001; // 1bp tolerance for model engine call put premia
1074	// vs. black premia
1075	// note that we use the real market conventions here (i.e. 2 fixing days),
1076	// different from the calibration approach where 0 fixing days must be used.
1077	// therefore higher errors compared to the calibration results are expected.
1078
1079	BOOST_TEST_MESSAGE("Testing Markov functional vanilla engines...");
1080
1081	Date savedEvalDate = Settings::instance().evaluationDate();
1082	Date referenceDate(14, November, 2012);
1083	Settings::instance().evaluationDate() = referenceDate;
1084
1085	Handle<YieldTermStructure> flatYts_ = flatYts();
1086	Handle<YieldTermStructure> md0Yts_ = md0Yts();
1087	Handle<SwaptionVolatilityStructure> flatSwaptionVts_ = flatSwaptionVts();
1088	Handle<SwaptionVolatilityStructure> md0SwaptionVts_ = md0SwaptionVts();
1089	Handle<OptionletVolatilityStructure> flatOptionletVts_ = flatOptionletVts();
1090	Handle<OptionletVolatilityStructure> md0OptionletVts_ = md0OptionletVts();
1091
1092	ext::shared_ptr<SwapIndex> swapIndexBase(
1093	new EuriborSwapIsdaFixA(1 * Years));
1094
1095	std::vector<Date> volStepDates;
1096	std::vector<Real> vols = {1.0};
1097
1098	// use a grid with few points for the check here
1099	std::vector<Real> money = {0.1, 0.25, 0.50, 0.75, 1.0, 1.25, 1.50, 2.0, 5.0};
1100
1101	// Calibration Basket 1 / flat yts, vts
1102
1103	ext::shared_ptr<IborIndex> iborIndex1(new Euribor(6 * Months, flatYts_));
1104
1105	ext::shared_ptr<MarkovFunctional> mf1(new MarkovFunctional(
1106	flatYts_, 0.01, volStepDates, vols, flatSwaptionVts_,
1107	expiriesCalBasket1(), tenorsCalBasket1(), swapIndexBase,
1108	MarkovFunctional::ModelSettings()
1109	.withYGridPoints(n: 64)
1110	.withYStdDevs(s: 7.0)
1111	.withGaussHermitePoints(n: 32)
1112	.withDigitalGap(d: 1e-5)
1113	.withMarketRateAccuracy(a: 1e-7)
1114	.withLowerRateBound(l: 0.0)
1115	.withUpperRateBound(u: 2.0)
1116	.withSmileMoneynessCheckpoints(m: money)));
1117
1118	MarkovFunctional::ModelOutputs outputs1 = mf1->modelOutputs();
1119	// BOOST_TEST_MESSAGE(outputs1);
1120
1121	ext::shared_ptr<Gaussian1dSwaptionEngine> mfSwaptionEngine1(
1122	new Gaussian1dSwaptionEngine(mf1, 64, 7.0));
1123	ext::shared_ptr<BlackSwaptionEngine> blackSwaptionEngine1(
1124	new BlackSwaptionEngine(flatYts_, flatSwaptionVts_));
1125
1126	for (Size i = 0; i < outputs1.expiries_.size(); i++) {
1127	for (Size j = 0; j < outputs1.smileStrikes_[0].size(); j++) {
1128	ext::shared_ptr<VanillaSwap> underlyingCall =
1129	MakeVanillaSwap(outputs1.tenors_[i], iborIndex1,
1130	outputs1.smileStrikes_[i][j])
1131	.withEffectiveDate(
1132	TARGET().advance(outputs1.expiries_[i], n: 2, unit: Days))
1133	.receiveFixed(flag: false);
1134	ext::shared_ptr<VanillaSwap> underlyingPut =
1135	MakeVanillaSwap(outputs1.tenors_[i], iborIndex1,
1136	outputs1.smileStrikes_[i][j])
1137	.withEffectiveDate(
1138	TARGET().advance(outputs1.expiries_[i], n: 2, unit: Days))
1139	.receiveFixed(flag: true);
1140	ext::shared_ptr<Exercise> exercise(
1141	new EuropeanExercise(outputs1.expiries_[i]));
1142	Swaption swaptionC(underlyingCall, exercise);
1143	Swaption swaptionP(underlyingPut, exercise);
1144	swaptionC.setPricingEngine(blackSwaptionEngine1);
1145	swaptionP.setPricingEngine(blackSwaptionEngine1);
1146	Real blackPriceCall = swaptionC.NPV();
1147	Real blackPricePut = swaptionP.NPV();
1148	swaptionC.setPricingEngine(mfSwaptionEngine1);
1149	swaptionP.setPricingEngine(mfSwaptionEngine1);
1150	Real mfPriceCall = swaptionC.NPV();
1151	Real mfPricePut = swaptionP.NPV();
1152	if (fabs(x: blackPriceCall - mfPriceCall) > tol1)
1153	BOOST_ERROR(
1154	"Basket 1 / flat termstructures: Call premium market ("
1155	<< blackPriceCall << ") does not match model premium ("
1156	<< mfPriceCall << ")");
1157	if (fabs(x: blackPricePut - mfPricePut) > tol1)
1158	BOOST_ERROR(
1159	"Basket 1 / flat termstructures: Put premium market ("
1160	<< blackPricePut << ") does not match model premium ("
1161	<< mfPricePut << ")");
1162	}
1163	}
1164
1165	// Calibration Basket 2 / flat yts, vts
1166
1167	ext::shared_ptr<IborIndex> iborIndex2(new Euribor(6 * Months, flatYts_));
1168
1169	ext::shared_ptr<MarkovFunctional> mf2(new MarkovFunctional(
1170	flatYts_, 0.01, volStepDates, vols, flatOptionletVts_,
1171	expiriesCalBasket2(), iborIndex2,
1172	MarkovFunctional::ModelSettings()
1173	.withYGridPoints(n: 64)
1174	.withYStdDevs(s: 7.0)
1175	.withGaussHermitePoints(n: 16)
1176	.withDigitalGap(d: 1e-5)
1177	.withMarketRateAccuracy(a: 1e-7)
1178	.withLowerRateBound(l: 0.0)
1179	.withUpperRateBound(u: 2.0)
1180	.withSmileMoneynessCheckpoints(m: money)));
1181
1182	MarkovFunctional::ModelOutputs outputs2 = mf2->modelOutputs();
1183	// BOOST_TEST_MESSAGE(outputs2);
1184
1185	ext::shared_ptr<BlackCapFloorEngine> blackCapFloorEngine2(
1186	new BlackCapFloorEngine(flatYts_, flatOptionletVts_));
1187	ext::shared_ptr<Gaussian1dCapFloorEngine> mfCapFloorEngine2(
1188	new Gaussian1dCapFloorEngine(mf2, 64, 7.0));
1189	std::vector<CapFloor> c2 = {
1190	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.01),
1191	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.02),
1192	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.03),
1193	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.04),
1194	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.05),
1195	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.07),
1196	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex2, 0.10),
1197	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.01),
1198	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.02),
1199	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.03),
1200	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.04),
1201	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.05),
1202	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.07),
1203	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex2, 0.10)
1204	};
1205
1206	for (auto& i : c2) {
1207	i.setPricingEngine(blackCapFloorEngine2);
1208	Real blackPrice = i.NPV();
1209	i.setPricingEngine(mfCapFloorEngine2);
1210	Real mfPrice = i.NPV();
1211	if (fabs(x: blackPrice - mfPrice) > tol1)
1212	BOOST_ERROR(
1213	"Basket 2 / flat termstructures: Cap/Floor premium market ("
1214	<< blackPrice << ") does not match model premium (" << mfPrice
1215	<< ")");
1216	}
1217
1218	// Calibration Basket 1 / real yts, vts
1219
1220	ext::shared_ptr<IborIndex> iborIndex3(new Euribor(6 * Months, md0Yts_));
1221
1222	ext::shared_ptr<MarkovFunctional> mf3(new MarkovFunctional(
1223	md0Yts_, 0.01, volStepDates, vols, md0SwaptionVts_,
1224	expiriesCalBasket1(), tenorsCalBasket1(), swapIndexBase,
1225	MarkovFunctional::ModelSettings()
1226	.withYGridPoints(n: 64)
1227	.withYStdDevs(s: 7.0)
1228	.withGaussHermitePoints(n: 32)
1229	.withDigitalGap(d: 1e-5)
1230	.withMarketRateAccuracy(a: 1e-7)
1231	.withLowerRateBound(l: 0.0)
1232	.withUpperRateBound(u: 2.0)
1233	.withSmileMoneynessCheckpoints(m: money)));
1234
1235	ext::shared_ptr<Gaussian1dSwaptionEngine> mfSwaptionEngine3(
1236	new Gaussian1dSwaptionEngine(mf3, 64, 7.0));
1237	ext::shared_ptr<BlackSwaptionEngine> blackSwaptionEngine3(
1238	new BlackSwaptionEngine(md0Yts_, md0SwaptionVts_));
1239
1240	MarkovFunctional::ModelOutputs outputs3 = mf3->modelOutputs();
1241	// BOOST_TEST_MESSAGE(outputs3);
1242
1243	for (Size i = 0; i < outputs3.expiries_.size(); i++) {
1244	for (Size j = 0; j < outputs3.smileStrikes_[0].size(); j++) {
1245	ext::shared_ptr<VanillaSwap> underlyingCall =
1246	MakeVanillaSwap(outputs3.tenors_[i], iborIndex3,
1247	outputs3.smileStrikes_[i][j])
1248	.withEffectiveDate(
1249	TARGET().advance(outputs3.expiries_[i], n: 2, unit: Days))
1250	.receiveFixed(flag: false);
1251	ext::shared_ptr<VanillaSwap> underlyingPut =
1252	MakeVanillaSwap(outputs3.tenors_[i], iborIndex3,
1253	outputs3.smileStrikes_[i][j])
1254	.withEffectiveDate(
1255	TARGET().advance(outputs3.expiries_[i], n: 2, unit: Days))
1256	.receiveFixed(flag: true);
1257	ext::shared_ptr<Exercise> exercise(
1258	new EuropeanExercise(outputs3.expiries_[i]));
1259	Swaption swaptionC(underlyingCall, exercise);
1260	Swaption swaptionP(underlyingPut, exercise);
1261	swaptionC.setPricingEngine(blackSwaptionEngine3);
1262	swaptionP.setPricingEngine(blackSwaptionEngine3);
1263	Real blackPriceCall = swaptionC.NPV();
1264	Real blackPricePut = swaptionP.NPV();
1265	swaptionC.setPricingEngine(mfSwaptionEngine3);
1266	swaptionP.setPricingEngine(mfSwaptionEngine3);
1267	Real mfPriceCall = swaptionC.NPV();
1268	Real mfPricePut = swaptionP.NPV();
1269	Real smileCorrectionCall =
1270	(outputs3.marketCallPremium_[i][j] -
1271	outputs3.marketRawCallPremium_[i][j]); // we can not expect to
1272	// match the black
1273	// scholes price where
1274	// the smile is adjusted
1275	Real smileCorrectionPut = (outputs3.marketPutPremium_[i][j] -
1276	outputs3.marketRawPutPremium_[i][j]);
1277	if (fabs(x: blackPriceCall - mfPriceCall + smileCorrectionCall) > tol1)
1278	BOOST_ERROR(
1279	"Basket 1 / real termstructures: Call premium market ("
1280	<< blackPriceCall << ") does not match model premium ("
1281	<< mfPriceCall << ")");
1282	if (fabs(x: blackPricePut - mfPricePut + smileCorrectionPut) > tol1)
1283	BOOST_ERROR(
1284	"Basket 1 / real termstructures: Put premium market ("
1285	<< blackPricePut << ") does not match model premium ("
1286	<< mfPricePut << ")");
1287	}
1288	}
1289
1290	// Calibration Basket 2 / real yts, vts
1291
1292	ext::shared_ptr<IborIndex> iborIndex4(new Euribor(6 * Months, md0Yts_));
1293
1294	ext::shared_ptr<MarkovFunctional> mf4(new MarkovFunctional(
1295	md0Yts_, 0.01, volStepDates, vols, md0OptionletVts_,
1296	expiriesCalBasket2(), iborIndex4,
1297	MarkovFunctional::ModelSettings()
1298	.withYGridPoints(n: 64)
1299	.withYStdDevs(s: 7.0)
1300	.withGaussHermitePoints(n: 32)
1301	.withDigitalGap(d: 1e-5)
1302	.withMarketRateAccuracy(a: 1e-7)
1303	.withLowerRateBound(l: 0.0)
1304	.withUpperRateBound(u: 2.0)
1305	.withSmileMoneynessCheckpoints(m: money)
1306	//.withAdjustments(MarkovFunctional::ModelSettings::KahaleSmile
1307	// |
1308	// MarkovFunctional::ModelSettings::KahaleInterpolation
1309	//| MarkovFunctional::ModelSettings::KahaleExponentialExtrapolation
1310	//)
1311	));
1312
1313	MarkovFunctional::ModelOutputs outputs4 = mf4->modelOutputs();
1314	// BOOST_TEST_MESSAGE(outputs4);
1315
1316	ext::shared_ptr<BlackCapFloorEngine> blackCapFloorEngine4(
1317	new BlackCapFloorEngine(md0Yts_, md0OptionletVts_));
1318	ext::shared_ptr<Gaussian1dCapFloorEngine> mfCapFloorEngine4(
1319	new Gaussian1dCapFloorEngine(mf4, 64, 7.0));
1320
1321	std::vector<CapFloor> c4 = {
1322	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex4, 0.01),
1323	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex4, 0.02),
1324	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex4, 0.03),
1325	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex4, 0.04),
1326	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex4, 0.05),
1327	MakeCapFloor(CapFloor::Cap, 5 * Years, iborIndex4, 0.06),
1328	// //exclude because caplet stripper fails for this strike
1329	// MakeCapFloor(CapFloor::Cap,5*Years,iborIndex4,0.10),
1330	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex4, 0.01),
1331	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex4, 0.02),
1332	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex4, 0.03),
1333	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex4, 0.04),
1334	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex4, 0.05),
1335	MakeCapFloor(CapFloor::Floor, 5 * Years, iborIndex4, 0.06)
1336	// //exclude because caplet stripper fails for this strike
1337	// MakeCapFloor(CapFloor::Floor,5*Years,iborIndex4,0.10)
1338	};
1339
1340	for (auto& i : c4) {
1341	i.setPricingEngine(blackCapFloorEngine4);
1342	Real blackPrice = i.NPV();
1343	i.setPricingEngine(mfCapFloorEngine4);
1344	Real mfPrice = i.NPV();
1345	if (fabs(x: blackPrice - mfPrice) > tol1)
1346	BOOST_ERROR(
1347	"Basket 2 / real termstructures: Cap/Floor premium market ("
1348	<< blackPrice << ") does not match model premium (" << mfPrice
1349	<< ")");
1350	}
1351
1352	Settings::instance().evaluationDate() = savedEvalDate;
1353	}
1354
1355	void MarkovFunctionalTest::testCalibrationTwoInstrumentSets() {
1356
1357	const Real tol1 = 0.1; // 0.1 times vega tolerance for model vs. market in
1358	// second instrument set
1359	BOOST_TEST_MESSAGE(
1360	"Testing Markov functional calibration to two instrument sets...");
1361
1362	Date savedEvalDate = Settings::instance().evaluationDate();
1363	Date referenceDate(14, November, 2012);
1364	Settings::instance().evaluationDate() = referenceDate;
1365
1366	Handle<YieldTermStructure> flatYts_ = flatYts();
1367	Handle<YieldTermStructure> md0Yts_ = md0Yts();
1368	Handle<SwaptionVolatilityStructure> flatSwaptionVts_ = flatSwaptionVts();
1369	Handle<SwaptionVolatilityStructure> md0SwaptionVts_ = md0SwaptionVts();
1370	Handle<OptionletVolatilityStructure> flatOptionletVts_ = flatOptionletVts();
1371	Handle<OptionletVolatilityStructure> md0OptionletVts_ = md0OptionletVts();
1372
1373	ext::shared_ptr<SwapIndex> swapIndexBase(
1374	new EuriborSwapIsdaFixA(1 * Years));
1375
1376	std::vector<Date> volStepDates = {
1377	TARGET().advance(date: referenceDate, period: 1 * Years),
1378	TARGET().advance(date: referenceDate, period: 2 * Years),
1379	TARGET().advance(date: referenceDate, period: 3 * Years),
1380	TARGET().advance(date: referenceDate, period: 4 * Years)
1381	};
1382	std::vector<Real> vols = {1.0, 1.0, 1.0, 1.0, 1.0};
1383	std::vector<Real> money = {0.1, 0.25, 0.50, 0.75, 1.0, 1.25, 1.50, 2.0, 5.0};
1384
1385	LevenbergMarquardt om;
1386	// ConjugateGradient om;
1387	EndCriteria ec(1000, 500, 1e-2, 1e-2, 1e-2);
1388
1389	// Calibration Basket 1 / flat yts, vts / Secondary calibration set consists
1390	// of coterminal swaptions
1391
1392	ext::shared_ptr<IborIndex> iborIndex1(new Euribor(6 * Months, flatYts_));
1393
1394	std::vector<ext::shared_ptr<BlackCalibrationHelper> > calibrationHelper1;
1395	std::vector<Real> calibrationHelperVols1 = {0.20, 0.20, 0.20, 0.20};
1396
1397	calibrationHelper1.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1398	new SwaptionHelper(1 * Years, 4 * Years,
1399	Handle<Quote>(ext::shared_ptr<Quote>(
1400	new SimpleQuote(calibrationHelperVols1[0]))),
1401	iborIndex1, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1402	flatYts_)));
1403	calibrationHelper1.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1404	new SwaptionHelper(2 * Years, 3 * Years,
1405	Handle<Quote>(ext::shared_ptr<Quote>(
1406	new SimpleQuote(calibrationHelperVols1[1]))),
1407	iborIndex1, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1408	flatYts_)));
1409	calibrationHelper1.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1410	new SwaptionHelper(3 * Years, 2 * Years,
1411	Handle<Quote>(ext::shared_ptr<Quote>(
1412	new SimpleQuote(calibrationHelperVols1[2]))),
1413	iborIndex1, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1414	flatYts_)));
1415	calibrationHelper1.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1416	new SwaptionHelper(4 * Years, 1 * Years,
1417	Handle<Quote>(ext::shared_ptr<Quote>(
1418	new SimpleQuote(calibrationHelperVols1[3]))),
1419	iborIndex1, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1420	flatYts_)));
1421
1422	ext::shared_ptr<MarkovFunctional> mf1(new MarkovFunctional(
1423	flatYts_, 0.01, volStepDates, vols, flatSwaptionVts_,
1424	expiriesCalBasket1(), tenorsCalBasket1(), swapIndexBase,
1425	MarkovFunctional::ModelSettings()
1426	.withYGridPoints(n: 64)
1427	.withYStdDevs(s: 7.0)
1428	.withGaussHermitePoints(n: 32)
1429	.withDigitalGap(d: 1e-5)
1430	.withMarketRateAccuracy(a: 1e-7)
1431	.withLowerRateBound(l: 0.0)
1432	.withUpperRateBound(u: 2.0)
1433	.withSmileMoneynessCheckpoints(m: money)));
1434
1435	ext::shared_ptr<Gaussian1dSwaptionEngine> mfSwaptionEngine1(
1436	new Gaussian1dSwaptionEngine(mf1, 64, 7.0));
1437	calibrationHelper1[0]->setPricingEngine(mfSwaptionEngine1);
1438	calibrationHelper1[1]->setPricingEngine(mfSwaptionEngine1);
1439	calibrationHelper1[2]->setPricingEngine(mfSwaptionEngine1);
1440	calibrationHelper1[3]->setPricingEngine(mfSwaptionEngine1);
1441
1442	mf1->calibrate(helpers: calibrationHelper1, method&: om, endCriteria: ec);
1443
1444	// std::cout << "Calibrated parameters 1: ";
1445	// Array params1 = mf1->params();
1446	// for(Size i=0;i<params1.size();i++) std::cout << params1[i] << ";";
1447	// std::cout << std::endl;
1448
1449	std::vector<Swaption> ch1;
1450	ch1.push_back(
1451	x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1452	new EuriborSwapIsdaFixA(4 * Years, flatYts_)),
1453	1 * Years));
1454	ch1.push_back(
1455	x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1456	new EuriborSwapIsdaFixA(3 * Years, flatYts_)),
1457	2 * Years));
1458	ch1.push_back(
1459	x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1460	new EuriborSwapIsdaFixA(2 * Years, flatYts_)),
1461	3 * Years));
1462	ch1.push_back(
1463	x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1464	new EuriborSwapIsdaFixA(1 * Years, flatYts_)),
1465	4 * Years));
1466
1467	for (Size i = 0; i < ch1.size(); i++) {
1468	ext::shared_ptr<BlackSwaptionEngine> blackEngine(
1469	new BlackSwaptionEngine(flatYts_, calibrationHelperVols1[i]));
1470	ch1[i].setPricingEngine(blackEngine);
1471	Real blackPrice = ch1[i].NPV();
1472	Real blackVega = ch1[i].result<Real>(tag: "vega");
1473	ch1[i].setPricingEngine(mfSwaptionEngine1);
1474	Real mfPrice = ch1[i].NPV();
1475	if (fabs(x: blackPrice - mfPrice) / blackVega > tol1)
1476	BOOST_TEST_MESSAGE("Basket 1 / flat yts, vts: Secondary instrument set "
1477	"calibration failed for instrument #"
1478	<< i << " black premium is " << blackPrice
1479	<< " while model premium is " << mfPrice
1480	<< " (market vega is " << blackVega << ")");
1481	}
1482
1483	// MarkovFunctional::ModelOutputs outputs1 = mf1->modelOutputs();
1484	// BOOST_TEST_MESSAGE(outputs1);
1485
1486	// Calibration Basket 1 / real yts, vts / Secondary calibration set consists
1487	// of coterminal swaptions
1488
1489	ext::shared_ptr<IborIndex> iborIndex2(new Euribor(6 * Months, md0Yts_));
1490
1491	ext::shared_ptr<MarkovFunctional> mf2(new MarkovFunctional(
1492	md0Yts_, 0.01, volStepDates, vols, md0SwaptionVts_,
1493	expiriesCalBasket1(), tenorsCalBasket1(), swapIndexBase,
1494	MarkovFunctional::ModelSettings()
1495	.withYGridPoints(n: 64)
1496	.withYStdDevs(s: 7.0)
1497	.withGaussHermitePoints(n: 32)
1498	.withDigitalGap(d: 1e-5)
1499	.withMarketRateAccuracy(a: 1e-7)
1500	.withLowerRateBound(l: 0.0)
1501	.withUpperRateBound(u: 2.0)
1502	.withSmileMoneynessCheckpoints(m: money)));
1503
1504	std::vector<ext::shared_ptr<BlackCalibrationHelper> > calibrationHelper2;
1505	std::vector<Real> calibrationHelperVols2;
1506	calibrationHelperVols2.push_back(x: md0SwaptionVts_->volatility(
1507	optionTenor: 1 * Years, swapTenor: 4 * Years,
1508	strike: ext::dynamic_pointer_cast<SwaptionVolatilityCube>(
1509	r: md0SwaptionVts_.currentLink())->atmStrike(optionTenor: 1 * Years, swapTenor: 4 * Years)));
1510	calibrationHelperVols2.push_back(x: md0SwaptionVts_->volatility(
1511	optionTenor: 2 * Years, swapTenor: 3 * Years,
1512	strike: ext::dynamic_pointer_cast<SwaptionVolatilityCube>(
1513	r: md0SwaptionVts_.currentLink())->atmStrike(optionTenor: 2 * Years, swapTenor: 3 * Years)));
1514	calibrationHelperVols2.push_back(x: md0SwaptionVts_->volatility(
1515	optionTenor: 3 * Years, swapTenor: 2 * Years,
1516	strike: ext::dynamic_pointer_cast<SwaptionVolatilityCube>(
1517	r: md0SwaptionVts_.currentLink())->atmStrike(optionTenor: 3 * Years, swapTenor: 2 * Years)));
1518	calibrationHelperVols2.push_back(x: md0SwaptionVts_->volatility(
1519	optionTenor: 4 * Years, swapTenor: 1 * Years,
1520	strike: ext::dynamic_pointer_cast<SwaptionVolatilityCube>(
1521	r: md0SwaptionVts_.currentLink())->atmStrike(optionTenor: 4 * Years, swapTenor: 1 * Years)));
1522
1523	calibrationHelper2.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1524	new SwaptionHelper(1 * Years, 4 * Years,
1525	Handle<Quote>(ext::shared_ptr<Quote>(
1526	new SimpleQuote(calibrationHelperVols2[0]))),
1527	iborIndex2, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1528	md0Yts_)));
1529	calibrationHelper2.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1530	new SwaptionHelper(2 * Years, 3 * Years,
1531	Handle<Quote>(ext::shared_ptr<Quote>(
1532	new SimpleQuote(calibrationHelperVols2[1]))),
1533	iborIndex2, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1534	md0Yts_)));
1535	calibrationHelper2.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1536	new SwaptionHelper(3 * Years, 2 * Years,
1537	Handle<Quote>(ext::shared_ptr<Quote>(
1538	new SimpleQuote(calibrationHelperVols2[2]))),
1539	iborIndex2, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1540	md0Yts_)));
1541	calibrationHelper2.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
1542	new SwaptionHelper(4 * Years, 1 * Years,
1543	Handle<Quote>(ext::shared_ptr<Quote>(
1544	new SimpleQuote(calibrationHelperVols2[3]))),
1545	iborIndex2, 1 * Years, Thirty360(Thirty360::BondBasis), Actual360(),
1546	md0Yts_)));
1547
1548	ext::shared_ptr<Gaussian1dSwaptionEngine> mfSwaptionEngine2(
1549	new Gaussian1dSwaptionEngine(mf2, 64, 7.0));
1550	calibrationHelper2[0]->setPricingEngine(mfSwaptionEngine2);
1551	calibrationHelper2[1]->setPricingEngine(mfSwaptionEngine2);
1552	calibrationHelper2[2]->setPricingEngine(mfSwaptionEngine2);
1553	calibrationHelper2[3]->setPricingEngine(mfSwaptionEngine2);
1554
1555	mf2->calibrate(helpers: calibrationHelper2, method&: om, endCriteria: ec);
1556
1557	// std::cout << "Calibrated parameters 2: ";
1558	// Array params2 = mf2->params();
1559	// for(Size i=0;i<params2.size();i++) std::cout << params2[i] << ";";
1560	// std::cout << std::endl;
1561
1562	std::vector<Swaption> ch2;
1563	ch2.push_back(x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1564	new EuriborSwapIsdaFixA(4 * Years, md0Yts_)),
1565	1 * Years));
1566	ch2.push_back(x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1567	new EuriborSwapIsdaFixA(3 * Years, md0Yts_)),
1568	2 * Years));
1569	ch2.push_back(x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1570	new EuriborSwapIsdaFixA(2 * Years, md0Yts_)),
1571	3 * Years));
1572	ch2.push_back(x: MakeSwaption(ext::shared_ptr<SwapIndex>(
1573	new EuriborSwapIsdaFixA(1 * Years, md0Yts_)),
1574	4 * Years));
1575
1576	for (Size i = 0; i < ch2.size(); i++) {
1577	ext::shared_ptr<BlackSwaptionEngine> blackEngine(
1578	new BlackSwaptionEngine(md0Yts_, calibrationHelperVols2[i]));
1579	ch2[i].setPricingEngine(blackEngine);
1580	Real blackPrice = ch2[i].NPV();
1581	Real blackVega = ch2[i].result<Real>(tag: "vega");
1582	ch2[i].setPricingEngine(mfSwaptionEngine2);
1583	Real mfPrice = ch2[i].NPV();
1584	if (fabs(x: blackPrice - mfPrice) / blackVega > tol1)
1585	BOOST_TEST_MESSAGE("Basket 1 / real yts, vts: Secondary instrument set "
1586	"calibration failed for instrument #"
1587	<< i << " black premium is " << blackPrice
1588	<< " while model premium is " << mfPrice
1589	<< " (market vega is " << blackVega << ")");
1590	}
1591
1592	// MarkovFunctional::ModelOutputs outputs2 = mf2->modelOutputs();
1593	// BOOST_TEST_MESSAGE(outputs2);
1594
1595	Settings::instance().evaluationDate() = savedEvalDate;
1596	}
1597
1598	void MarkovFunctionalTest::testBermudanSwaption() {
1599
1600	Real tol0 = 0.0001; // 1bp tolerance against cached values
1601
1602	BOOST_TEST_MESSAGE("Testing Markov functional Bermudan swaption engine...");
1603
1604	Date savedEvalDate = Settings::instance().evaluationDate();
1605	Date referenceDate(14, November, 2012);
1606	Settings::instance().evaluationDate() = referenceDate;
1607
1608	Handle<YieldTermStructure> flatYts_ = flatYts();
1609	Handle<YieldTermStructure> md0Yts_ = md0Yts();
1610	Handle<SwaptionVolatilityStructure> flatSwaptionVts_ = flatSwaptionVts();
1611	Handle<SwaptionVolatilityStructure> md0SwaptionVts_ = md0SwaptionVts();
1612	Handle<OptionletVolatilityStructure> flatOptionletVts_ = flatOptionletVts();
1613	Handle<OptionletVolatilityStructure> md0OptionletVts_ = md0OptionletVts();
1614
1615	ext::shared_ptr<SwapIndex> swapIndexBase(
1616	new EuriborSwapIsdaFixA(1 * Years));
1617
1618	std::vector<Date> volStepDates;
1619	std::vector<Real> vols = {1.0};
1620
1621	ext::shared_ptr<IborIndex> iborIndex1(new Euribor(6 * Months, md0Yts_));
1622
1623	ext::shared_ptr<MarkovFunctional> mf1(
1624	new MarkovFunctional(md0Yts_, 0.01, volStepDates, vols, md0SwaptionVts_,
1625	expiriesCalBasket3(), tenorsCalBasket3(),
1626	swapIndexBase, MarkovFunctional::ModelSettings()
1627	.withYGridPoints(n: 32)
1628	.withYStdDevs(s: 7.0)
1629	.withGaussHermitePoints(n: 16)
1630	.withMarketRateAccuracy(a: 1e-7)
1631	.withDigitalGap(d: 1e-5)
1632	.withLowerRateBound(l: 0.0)
1633	.withUpperRateBound(u: 2.0)));
1634
1635	ext::shared_ptr<PricingEngine> mfSwaptionEngine1(
1636	new Gaussian1dSwaptionEngine(mf1, 64, 7.0));
1637
1638	ext::shared_ptr<VanillaSwap> underlyingCall =
1639	MakeVanillaSwap(10 * Years, iborIndex1, 0.03)
1640	.withEffectiveDate(TARGET().advance(referenceDate, n: 2, unit: Days))
1641	//.withNominal(100000000.0)
1642	.receiveFixed(flag: false);
1643
1644	std::vector<ext::shared_ptr<Exercise> > europeanExercises;
1645	std::vector<Date> expiries = expiriesCalBasket3();
1646	std::vector<Swaption> europeanSwaptions;
1647	for (Size i = 0; i < expiries.size(); i++) {
1648	europeanExercises.push_back(
1649	x: ext::shared_ptr<Exercise>(new EuropeanExercise(expiries[i])));
1650	europeanSwaptions.emplace_back(args&: underlyingCall, args&: europeanExercises[i]);
1651	europeanSwaptions.back().setPricingEngine(mfSwaptionEngine1);
1652	}
1653
1654	ext::shared_ptr<Exercise> bermudanExercise(
1655	new BermudanExercise(expiries));
1656	Swaption bermudanSwaption(underlyingCall, bermudanExercise);
1657	bermudanSwaption.setPricingEngine(mfSwaptionEngine1);
1658
1659	Real cachedValues[] = { 0.0030757, 0.0107344, 0.0179862,
1660	0.0225881, 0.0243215, 0.0229148,
1661	0.0191415, 0.0139035, 0.0076354 };
1662	Real cachedValue = 0.0327776;
1663
1664	for (Size i = 0; i < expiries.size(); i++) {
1665	Real npv = europeanSwaptions[i].NPV();
1666	if (fabs(x: npv - cachedValues[i]) > tol0)
1667	BOOST_ERROR("European swaption value ("
1668	<< npv << ") deviates from cached value ("
1669	<< cachedValues[i] << ")");
1670	}
1671
1672	Real npv = bermudanSwaption.NPV();
1673	if (fabs(x: npv - cachedValue) > tol0)
1674	BOOST_ERROR("Bermudan swaption value ("
1675	<< npv << ") deviates from cached value (" << cachedValue
1676	<< ")");
1677
1678	Settings::instance().evaluationDate() = savedEvalDate;
1679	}
1680
1681	test_suite *MarkovFunctionalTest::suite(SpeedLevel speed) {
1682	auto* suite = BOOST_TEST_SUITE("Markov functional model tests");
1683
1684	suite->add(QUANTLIB_TEST_CASE(&MarkovFunctionalTest::testMfStateProcess));
1685	suite->add(QUANTLIB_TEST_CASE(&MarkovFunctionalTest::testKahaleSmileSection));
1686	suite->add(QUANTLIB_TEST_CASE(&MarkovFunctionalTest::testBermudanSwaption));
1687
1688	if (speed <= Fast) {
1689	suite->add(QUANTLIB_TEST_CASE(&MarkovFunctionalTest::testCalibrationTwoInstrumentSets));
1690	}
1691
1692	if (speed == Slow) {
1693	suite->add(QUANTLIB_TEST_CASE(&MarkovFunctionalTest::testCalibrationOneInstrumentSet));
1694	suite->add(QUANTLIB_TEST_CASE(&MarkovFunctionalTest::testVanillaEngines));
1695	}
1696
1697	return suite;
1698	}
1699