[go: up one dir, main page]
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2005, 2008 Klaus Spanderen
5 Copyright (C) 2007 StatPro Italia srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include "batesmodel.hpp"
22#include "utilities.hpp"
23#include <ql/time/calendars/target.hpp>
24#include <ql/processes/batesprocess.hpp>
25#include <ql/processes/merton76process.hpp>
26#include <ql/instruments/europeanoption.hpp>
27#include <ql/time/daycounters/actualactual.hpp>
28#include <ql/termstructures/yield/flatforward.hpp>
29#include <ql/termstructures/yield/zerocurve.hpp>
30#include <ql/pricingengines/blackformula.hpp>
31#include <ql/math/optimization/levenbergmarquardt.hpp>
32#include <ql/pricingengines/vanilla/batesengine.hpp>
33#include <ql/pricingengines/vanilla/jumpdiffusionengine.hpp>
34#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
35#include <ql/pricingengines/vanilla/mceuropeanhestonengine.hpp>
36#include <ql/pricingengines/vanilla/fdbatesvanillaengine.hpp>
37#include <ql/models/equity/batesmodel.hpp>
38#include <ql/models/equity/hestonmodelhelper.hpp>
39#include <ql/time/period.hpp>
40#include <ql/quotes/simplequote.hpp>
41
42using namespace QuantLib;
43using namespace boost::unit_test_framework;
44
45namespace bates_model_test {
46
47 Real getCalibrationError(
48 std::vector<ext::shared_ptr<BlackCalibrationHelper> > & options) {
49 Real sse = 0;
50 for (auto& option : options) {
51 const Real diff = option->calibrationError() * 100.0;
52 sse += diff*diff;
53 }
54 return sse;
55 }
56
57}
58
59
60void BatesModelTest::testAnalyticVsBlack() {
61
62 BOOST_TEST_MESSAGE("Testing analytic Bates engine against Black formula...");
63
64 Date settlementDate = Date::todaysDate();
65 Settings::instance().evaluationDate() = settlementDate;
66
67 DayCounter dayCounter = ActualActual(ActualActual::ISDA);
68 Date exerciseDate = settlementDate + 6*Months;
69
70 ext::shared_ptr<StrikedTypePayoff> payoff(
71 new PlainVanillaPayoff(Option::Put, 30));
72 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));
73
74 Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.1, dc: dayCounter));
75 Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.04, dc: dayCounter));
76 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(32.0)));
77
78 Real yearFraction = dayCounter.yearFraction(d1: settlementDate, d2: exerciseDate);
79 Real forwardPrice = s0->value()*std::exp(x: (0.1-0.04)*yearFraction);
80 Real expected = blackFormula(optionType: payoff->optionType(), strike: payoff->strike(),
81 forward: forwardPrice, stdDev: std::sqrt(x: 0.05*yearFraction)) *
82 std::exp(x: -0.1*yearFraction);
83 const Real v0 = 0.05;
84 const Real kappa = 5.0;
85 const Real theta = 0.05;
86 const Real sigma = 1.0e-4;
87 const Real rho = 0.0;
88 const Real lambda = 0.0001;
89 const Real nu = 0.0;
90 const Real delta = 0.0001;
91
92 VanillaOption option(payoff, exercise);
93
94 ext::shared_ptr<BatesProcess> process(
95 new BatesProcess(riskFreeTS, dividendTS, s0, v0,
96 kappa, theta, sigma, rho, lambda, nu, delta));
97
98 ext::shared_ptr<PricingEngine> engine(new BatesEngine(
99 ext::make_shared<BatesModel>(args&: process), 64));
100
101 option.setPricingEngine(engine);
102 Real calculated = option.NPV();
103
104 Real tolerance = 2.0e-7;
105 Real error = std::fabs(x: calculated - expected);
106 if (error > tolerance) {
107 BOOST_ERROR("failed to reproduce Black price with BatesEngine"
108 << std::fixed
109 << "\n calculated: " << calculated
110 << "\n expected: " << expected
111 << std::scientific
112 << "\n error: " << error);
113 }
114
115 engine = ext::shared_ptr<PricingEngine>(new BatesDetJumpEngine(
116 ext::make_shared<BatesDetJumpModel>(
117 args&: process, args: 1.0, args: 0.0001), 64));
118
119 option.setPricingEngine(engine);
120 calculated = option.NPV();
121
122 error = std::fabs(x: calculated - expected);
123 if (error > tolerance) {
124 BOOST_ERROR("failed to reproduce Black price with " \
125 "BatesDetJumpEngine"
126 << std::fixed
127 << "\n calculated: " << calculated
128 << "\n expected: " << expected
129 << std::scientific
130 << "\n error: " << error);
131 }
132
133 engine = ext::shared_ptr<PricingEngine>(new BatesDoubleExpEngine(
134 ext::make_shared<BatesDoubleExpModel>(
135 args&: process, args: 0.0001, args: 0.0001, args: 0.0001), 64));
136
137 option.setPricingEngine(engine);
138 calculated = option.NPV();
139
140 error = std::fabs(x: calculated - expected);
141 if (error > tolerance) {
142 BOOST_ERROR("failed to reproduce Black price with BatesDoubleExpEngine"
143 << std::fixed
144 << "\n calculated: " << calculated
145 << "\n expected: " << expected
146 << std::scientific
147 << "\n error: " << error);
148 }
149
150 engine = ext::shared_ptr<PricingEngine>(new BatesDoubleExpDetJumpEngine(
151 ext::make_shared<BatesDoubleExpDetJumpModel>(
152
153 args&: process, args: 0.0001, args: 0.0001, args: 0.0001, args: 0.5, args: 1.0, args: 0.0001), 64));
154
155 option.setPricingEngine(engine);
156 calculated = option.NPV();
157
158 error = std::fabs(x: calculated - expected);
159 if (error > tolerance) {
160 BOOST_ERROR("failed to reproduce Black price with " \
161 "BatesDoubleExpDetJumpEngine"
162 << std::fixed
163 << "\n calculated: " << calculated
164 << "\n expected: " << expected
165 << std::scientific
166 << "\n error: " << error);
167 }
168}
169
170
171void BatesModelTest::testAnalyticAndMcVsJumpDiffusion() {
172
173 BOOST_TEST_MESSAGE("Testing analytic Bates engine against Merton-76 engine...");
174
175 Date settlementDate = Date::todaysDate();
176 Settings::instance().evaluationDate() = settlementDate;
177
178 DayCounter dayCounter = ActualActual(ActualActual::ISDA);
179
180 ext::shared_ptr<StrikedTypePayoff> payoff(
181 new PlainVanillaPayoff(Option::Put, 95));
182
183 Handle<YieldTermStructure> riskFreeTS(flatRate(forward: 0.1, dc: dayCounter));
184 Handle<YieldTermStructure> dividendTS(flatRate(forward: 0.04, dc: dayCounter));
185 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
186
187 Real v0 = 0.0433;
188 ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(std::sqrt(x: v0)));
189 ext::shared_ptr<BlackVolTermStructure> volTS =
190 flatVol(today: settlementDate, volatility: vol, dc: dayCounter);
191
192 const Real kappa = 0.5;
193 const Real theta = v0;
194 const Real sigma = 1.0e-4;
195 const Real rho = 0.0;
196
197 ext::shared_ptr<SimpleQuote> jumpIntensity(new SimpleQuote(2));
198 ext::shared_ptr<SimpleQuote> meanLogJump(new SimpleQuote(-0.2));
199 ext::shared_ptr<SimpleQuote> jumpVol(new SimpleQuote(0.2));
200
201 ext::shared_ptr<BatesProcess> batesProcess(new BatesProcess(
202 riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho,
203 jumpIntensity->value(), meanLogJump->value(), jumpVol->value()));
204
205 ext::shared_ptr<Merton76Process> mertonProcess(
206 new Merton76Process(s0, dividendTS, riskFreeTS,
207 Handle<BlackVolTermStructure>(volTS),
208 Handle<Quote>(jumpIntensity),
209 Handle<Quote>(meanLogJump),
210 Handle<Quote>(jumpVol)));
211
212 ext::shared_ptr<PricingEngine> batesEngine(new BatesEngine(
213 ext::make_shared<BatesModel>(args&: batesProcess), 160));
214
215 const Real mcTol = 0.1;
216 ext::shared_ptr<PricingEngine> mcBatesEngine =
217 MakeMCEuropeanHestonEngine<PseudoRandom>(batesProcess)
218 .withStepsPerYear(steps: 2)
219 .withAntitheticVariate()
220 .withAbsoluteTolerance(tolerance: mcTol)
221 .withSeed(seed: 1234);
222
223 ext::shared_ptr<PricingEngine> mertonEngine(
224 new JumpDiffusionEngine(mertonProcess, 1e-10, 1000));
225
226 for (Integer i=1; i<=5; i+=2) {
227 Date exerciseDate = settlementDate + i*Years;
228 ext::shared_ptr<Exercise> exercise(
229 new EuropeanExercise(exerciseDate));
230
231 VanillaOption batesOption(payoff, exercise);
232
233 batesOption.setPricingEngine(batesEngine);
234 Real calculated = batesOption.NPV();
235
236 batesOption.setPricingEngine(mcBatesEngine);
237 Real mcCalculated = batesOption.NPV();
238
239 EuropeanOption mertonOption(payoff, exercise);
240 mertonOption.setPricingEngine(mertonEngine);
241 Real expected = mertonOption.NPV();
242
243 Real tolerance = 2e-8;
244 Real relError = std::fabs(x: calculated - expected)/expected;
245 if (relError > tolerance) {
246 BOOST_FAIL("failed to reproduce Merton76 price with semi "
247 "analytic BatesEngine"
248 << std::fixed << std::setprecision(8)
249 << "\n calculated: " << calculated
250 << "\n expected: " << expected
251 << "\n rel. error: " << relError
252 << "\n tolerance: " << tolerance);
253 }
254
255 Real mcError = std::fabs(x: expected - mcCalculated);
256 if (mcError > 3*mcTol) {
257 BOOST_FAIL("failed to reproduce Merton76 price with Monte-Carlo "
258 "BatesEngine"
259 << std::fixed << std::setprecision(8)
260 << "\n calculated: " << mcCalculated
261 << "\n expected: " << expected
262 << "\n error: " << mcError
263 << "\n tolerance: " << mcTol);
264 }
265 }
266}
267
268namespace bates_model_test {
269 struct HestonModelData {
270 const char* const name;
271 Real v0;
272 Real kappa;
273 Real theta;
274 Real sigma;
275 Real rho;
276 Real r;
277 Real q;
278 };
279
280 HestonModelData hestonModels[] = {
281 // ADI finite difference schemes for option pricing in the
282 // Heston model with correlation, K.J. in t'Hout and S. Foulon,
283 {.name: "'t Hout case 1", .v0: 0.04, .kappa: 1.5, .theta: 0.04, .sigma: 0.3, .rho: -0.9, .r: 0.025, .q: 0.0},
284 // Efficient numerical methods for pricing American options under
285 // stochastic volatility, Samuli Ikonen and Jari Toivanen,
286 {.name: "Ikonen-Toivanen", .v0: 0.0625, .kappa: 5, .theta: 0.16, .sigma: 0.9, .rho: 0.1, .r: 0.1, .q: 0.0},
287 // Not-so-complex logarithms in the Heston model,
288 // Christian Kahl and Peter Jäckel
289 {.name: "Kahl-Jaeckel", .v0: 0.16, .kappa: 1.0, .theta: 0.16, .sigma: 2.0, .rho: -0.8, .r: 0.0, .q: 0.0},
290 // self defined test cases
291 {.name: "Equity case", .v0: 0.07, .kappa: 2.0, .theta: 0.04, .sigma: 0.55, .rho: -0.8, .r: 0.03, .q: 0.035 },
292 };
293}
294
295void BatesModelTest::testAnalyticVsMCPricing() {
296 BOOST_TEST_MESSAGE("Testing analytic Bates engine against Monte-Carlo "
297 "engine...");
298
299 using namespace bates_model_test;
300
301 Date settlementDate(30, March, 2007);
302 Settings::instance().evaluationDate() = settlementDate;
303
304 DayCounter dayCounter = ActualActual(ActualActual::ISDA);
305 Date exerciseDate(30, March, 2012);
306
307 ext::shared_ptr<StrikedTypePayoff> payoff(
308 new PlainVanillaPayoff(Option::Put, 100));
309 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));
310
311
312 for (auto& hestonModel : hestonModels) {
313 Handle<YieldTermStructure> riskFreeTS(flatRate(forward: hestonModel.r, dc: dayCounter));
314 Handle<YieldTermStructure> dividendTS(flatRate(forward: hestonModel.q, dc: dayCounter));
315 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
316
317 ext::shared_ptr<BatesProcess> batesProcess(new BatesProcess(
318 riskFreeTS, dividendTS, s0, hestonModel.v0, hestonModel.kappa, hestonModel.theta,
319 hestonModel.sigma, hestonModel.rho, 2.0, -0.2, 0.1));
320
321 const Real mcTolerance = 0.5;
322 ext::shared_ptr<PricingEngine> mcEngine =
323 MakeMCEuropeanHestonEngine<PseudoRandom>(batesProcess)
324 .withStepsPerYear(steps: 20)
325 .withAntitheticVariate()
326 .withAbsoluteTolerance(tolerance: mcTolerance)
327 .withSeed(seed: 1234);
328
329 ext::shared_ptr<BatesModel> batesModel(new BatesModel(batesProcess));
330
331 ext::shared_ptr<PricingEngine> fdEngine(
332 new FdBatesVanillaEngine(batesModel, 50, 100, 30));
333
334 ext::shared_ptr<PricingEngine> analyticEngine(
335 new BatesEngine(batesModel, 160));
336
337 VanillaOption option(payoff, exercise);
338
339 option.setPricingEngine(mcEngine);
340 const Real calculated = option.NPV();
341
342 option.setPricingEngine(analyticEngine);
343 const Real expected = option.NPV();
344
345 option.setPricingEngine(fdEngine);
346 const Real fdCalculated = option.NPV();
347
348 const Real mcError = std::fabs(x: calculated - expected);
349 if (mcError > 3*mcTolerance) {
350 BOOST_FAIL("failed to reproduce Monte-Carlo price for BatesEngine"
351 << "\n parameter: " << hestonModel.name << std::fixed
352 << std::setprecision(8) << "\n calculated: " << calculated
353 << "\n expected: " << expected << "\n error: " << mcError
354 << "\n tolerance: " << mcTolerance);
355 }
356 const Real fdTolerance = 0.2;
357 const Real fdError = std::fabs(x: fdCalculated - expected);
358 if (fdError > fdTolerance) {
359 BOOST_FAIL("failed to reproduce PIDE price for BatesEngine"
360 << "\n parameter: " << hestonModel.name << std::fixed
361 << std::setprecision(8) << "\n calculated: " << fdCalculated
362 << "\n expected: " << expected << "\n error: " << fdError
363 << "\n tolerance: " << fdTolerance);
364 }
365 }
366}
367
368void BatesModelTest::testDAXCalibration() {
369 /* this example is taken from A. Sepp
370 Pricing European-Style Options under Jump Diffusion Processes
371 with Stochstic Volatility: Applications of Fourier Transform
372 http://math.ut.ee/~spartak/papers/stochjumpvols.pdf
373 */
374
375 BOOST_TEST_MESSAGE(
376 "Testing Bates model calibration using DAX volatility data...");
377
378 using namespace bates_model_test;
379
380 Date settlementDate(5, July, 2002);
381 Settings::instance().evaluationDate() = settlementDate;
382
383 DayCounter dayCounter = Actual365Fixed();
384 Calendar calendar = TARGET();
385
386 Integer t[] = { 13, 41, 75, 165, 256, 345, 524, 703 };
387 Rate r[] = { 0.0357,0.0349,0.0341,0.0355,0.0359,0.0368,0.0386,0.0401 };
388
389 std::vector<Date> dates;
390 std::vector<Rate> rates;
391 dates.push_back(x: settlementDate);
392 rates.push_back(x: 0.0357);
393 for (Size i = 0; i < 8; ++i) {
394 dates.push_back(x: settlementDate + t[i]);
395 rates.push_back(x: r[i]);
396 }
397 Handle<YieldTermStructure> riskFreeTS(
398 ext::shared_ptr<YieldTermStructure>(
399 new ZeroCurve(dates, rates, dayCounter)));
400
401 Handle<YieldTermStructure> dividendTS(
402 flatRate(today: settlementDate, forward: 0.0, dc: dayCounter));
403
404 Volatility v[] =
405 { 0.6625,0.4875,0.4204,0.3667,0.3431,0.3267,0.3121,0.3121,
406 0.6007,0.4543,0.3967,0.3511,0.3279,0.3154,0.2984,0.2921,
407 0.5084,0.4221,0.3718,0.3327,0.3155,0.3027,0.2919,0.2889,
408 0.4541,0.3869,0.3492,0.3149,0.2963,0.2926,0.2819,0.2800,
409 0.4060,0.3607,0.3330,0.2999,0.2887,0.2811,0.2751,0.2775,
410 0.3726,0.3396,0.3108,0.2781,0.2788,0.2722,0.2661,0.2686,
411 0.3550,0.3277,0.3012,0.2781,0.2781,0.2661,0.2661,0.2681,
412 0.3428,0.3209,0.2958,0.2740,0.2688,0.2627,0.2580,0.2620,
413 0.3302,0.3062,0.2799,0.2631,0.2573,0.2533,0.2504,0.2544,
414 0.3343,0.2959,0.2705,0.2540,0.2504,0.2464,0.2448,0.2462,
415 0.3460,0.2845,0.2624,0.2463,0.2425,0.2385,0.2373,0.2422,
416 0.3857,0.2860,0.2578,0.2399,0.2357,0.2327,0.2312,0.2351,
417 0.3976,0.2860,0.2607,0.2356,0.2297,0.2268,0.2241,0.2320 };
418
419 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(4468.17)));
420 Real strike[] = { 3400,3600,3800,4000,4200,4400,
421 4500,4600,4800,5000,5200,5400,5600 };
422
423
424 Real v0 = 0.0433;
425 ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(std::sqrt(x: v0)));
426
427 const Real kappa = 1.0;
428 const Real theta = v0;
429 const Real sigma = 1.0;
430 const Real rho = 0.0;
431 const Real lambda = 1.1098;
432 const Real nu = -0.1285;
433 const Real delta = 0.1702;
434
435 ext::shared_ptr<BatesProcess> process(
436 new BatesProcess(riskFreeTS, dividendTS, s0, v0,
437 kappa, theta, sigma, rho, lambda, nu, delta));
438
439 ext::shared_ptr<BatesModel> batesModel(new BatesModel(process));
440
441 ext::shared_ptr<PricingEngine> batesEngine(
442 new BatesEngine(batesModel, 64));
443
444 std::vector<ext::shared_ptr<BlackCalibrationHelper> > options;
445
446 for (Size s = 0; s < 13; ++s) {
447 for (Size m = 0; m < 8; ++m) {
448 Handle<Quote> vol(ext::shared_ptr<Quote>(
449 new SimpleQuote(v[s*8+m])));
450
451 Period maturity((int)((t[m]+3)/7.), Weeks); // round to weeks
452
453 // this is the calibration helper for the bates models
454 options.push_back(x: ext::shared_ptr<BlackCalibrationHelper>(
455 new HestonModelHelper(maturity, calendar,
456 s0->value(), strike[s], vol,
457 riskFreeTS, dividendTS,
458 BlackCalibrationHelper::ImpliedVolError)));
459 options.back()->setPricingEngine(batesEngine);
460 }
461 }
462
463 // check calibration engine
464 LevenbergMarquardt om;
465 batesModel->calibrate(std::vector<ext::shared_ptr<CalibrationHelper> >(options.begin(), options.end()),
466 method&: om, endCriteria: EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8));
467
468 Real expected = 36.6;
469 Real calculated = getCalibrationError(options);
470
471 if (std::fabs(x: calculated - expected) > 2.5)
472 BOOST_ERROR("failed to calibrate the bates model"
473 << "\n calculated: " << calculated
474 << "\n expected: " << expected);
475
476 //check pricing of derived engines
477 std::vector<ext::shared_ptr<PricingEngine> > pricingEngines;
478
479 process = ext::make_shared<BatesProcess>(
480 args&: riskFreeTS, args&: dividendTS, args&: s0, args&: v0,
481 args: kappa, args: theta, args: sigma, args: rho, args: 1.0, args: -0.1, args: 0.1);
482
483 pricingEngines.push_back(x: ext::shared_ptr<PricingEngine>(
484 new BatesDetJumpEngine(
485 ext::make_shared<BatesDetJumpModel>(
486 args&: process), 64)) );
487
488 ext::shared_ptr<HestonProcess> hestonProcess(new HestonProcess(
489 riskFreeTS, dividendTS, s0, v0, kappa, theta, sigma, rho));
490
491 pricingEngines.push_back(x: ext::shared_ptr<PricingEngine>(
492 new BatesDoubleExpEngine(
493 ext::make_shared<BatesDoubleExpModel>(
494 args&: hestonProcess, args: 1.0), 64)) );
495
496 pricingEngines.push_back(x: ext::shared_ptr<PricingEngine>(
497 new BatesDoubleExpDetJumpEngine(
498 ext::make_shared<BatesDoubleExpDetJumpModel>(
499 args&: hestonProcess, args: 1.0), 64)) );
500
501 Real expectedValues[] = { 5896.37,
502 5499.29,
503 6497.89};
504
505 Real tolerance=0.1;
506 for (Size i = 0; i < pricingEngines.size(); ++i) {
507 for (auto& option : options) {
508 option->setPricingEngine(pricingEngines[i]);
509 }
510
511 Real calculated = std::fabs(x: getCalibrationError(options));
512 if (std::fabs(x: calculated - expectedValues[i]) > tolerance)
513 BOOST_ERROR("failed to calculated prices for derived Bates models"
514 << "\n calculated: " << calculated
515 << "\n expected: " << expectedValues[i]);
516 }
517}
518
519test_suite* BatesModelTest::suite() {
520 auto* suite = BOOST_TEST_SUITE("Bates model tests");
521 suite->add(QUANTLIB_TEST_CASE(&BatesModelTest::testAnalyticVsBlack));
522 suite->add(QUANTLIB_TEST_CASE(&BatesModelTest::testAnalyticAndMcVsJumpDiffusion));
523 suite->add(QUANTLIB_TEST_CASE(&BatesModelTest::testAnalyticVsMCPricing));
524 suite->add(QUANTLIB_TEST_CASE(&BatesModelTest::testDAXCalibration));
525 return suite;
526}
527

source code of quantlib/test-suite/batesmodel.cpp