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