1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2007, 2008, 2009, 2010 Klaus Spanderen
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 "utilities.hpp"
21#include "hybridhestonhullwhiteprocess.hpp"
22
23#include <ql/time/schedule.hpp>
24#include <ql/time/calendars/target.hpp>
25#include <ql/quotes/simplequote.hpp>
26#include <ql/instruments/europeanoption.hpp>
27#include <ql/instruments/impliedvolatility.hpp>
28#include <ql/processes/blackscholesprocess.hpp>
29#include <ql/processes/hybridhestonhullwhiteprocess.hpp>
30#include <ql/math/functional.hpp>
31#include <ql/math/randomnumbers/rngtraits.hpp>
32#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
33#include <ql/math/optimization/simplex.hpp>
34#include <ql/math/optimization/levenbergmarquardt.hpp>
35#include <ql/math/statistics/generalstatistics.hpp>
36#include <ql/math/statistics/sequencestatistics.hpp>
37#include <ql/math/statistics/incrementalstatistics.hpp>
38#include <ql/math/matrixutilities/svd.hpp>
39#include <ql/time/daycounters/actual360.hpp>
40#include <ql/time/daycounters/actual365fixed.hpp>
41#include <ql/methods/montecarlo/multipathgenerator.hpp>
42#include <ql/termstructures/yield/zerocurve.hpp>
43#include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
44#include <ql/models/equity/hestonmodelhelper.hpp>
45#include <ql/models/shortrate/onefactormodels/hullwhite.hpp>
46#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
47#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
48#include <ql/pricingengines/vanilla/analytich1hwengine.hpp>
49#include <ql/pricingengines/vanilla/mchestonhullwhiteengine.hpp>
50#include <ql/pricingengines/vanilla/analyticbsmhullwhiteengine.hpp>
51#include <ql/pricingengines/vanilla/analytichestonhullwhiteengine.hpp>
52#include <ql/pricingengines/vanilla/fdhestonvanillaengine.hpp>
53#include <ql/pricingengines/vanilla/fdhestonhullwhitevanillaengine.hpp>
54#include <cmath>
55
56using namespace QuantLib;
57using namespace boost::unit_test_framework;
58
59void HybridHestonHullWhiteProcessTest::testBsmHullWhiteEngine() {
60 BOOST_TEST_MESSAGE("Testing European option pricing for a BSM process"
61 " with one-factor Hull-White model...");
62
63 DayCounter dc = Actual365Fixed();
64
65 const Date today = Date::todaysDate();
66 const Date maturity = today + Period(20, Years);
67
68 Settings::instance().evaluationDate() = today;
69
70 const Handle<Quote> spot(
71 ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
72 ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.04));
73 const Handle<YieldTermStructure> qTS(flatRate(today, forward: qRate, dc));
74 ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0525));
75 const Handle<YieldTermStructure> rTS(flatRate(today, forward: rRate, dc));
76 ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.25));
77 const Handle<BlackVolTermStructure> volTS(flatVol(today, volatility: vol, dc));
78
79 ext::shared_ptr<HullWhite> hullWhiteModel(
80 new HullWhite(Handle<YieldTermStructure>(rTS), 0.00883, 0.00526));
81
82 ext::shared_ptr<BlackScholesMertonProcess> stochProcess(
83 new BlackScholesMertonProcess(spot, qTS, rTS, volTS));
84
85 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(maturity));
86
87 Real fwd = spot->value()*qTS->discount(d: maturity)/rTS->discount(d: maturity);
88 ext::shared_ptr<StrikedTypePayoff> payoff(new
89 PlainVanillaPayoff(Option::Call, fwd));
90
91 EuropeanOption option(payoff, exercise);
92
93 const Real tol = 1e-8;
94 const Real corr[] = {-0.75, -0.25, 0.0, 0.25, 0.75 };
95 const Volatility expectedVol[] = { 0.217064577, 0.243995801,
96 0.256402830, 0.268236596, 0.290461343 };
97
98 for (Size i=0; i < LENGTH(corr); ++i) {
99 ext::shared_ptr<PricingEngine> bsmhwEngine(
100 new AnalyticBSMHullWhiteEngine(corr[i], stochProcess,
101 hullWhiteModel));
102
103 option.setPricingEngine(bsmhwEngine);
104 const Real npv = option.NPV();
105
106 const Handle<BlackVolTermStructure> compVolTS(
107 flatVol(today, volatility: expectedVol[i], dc));
108
109 ext::shared_ptr<BlackScholesMertonProcess> bsProcess(
110 new BlackScholesMertonProcess(spot, qTS, rTS, compVolTS));
111 ext::shared_ptr<PricingEngine> bsEngine(
112 new AnalyticEuropeanEngine(bsProcess));
113
114 EuropeanOption comp(payoff, exercise);
115 comp.setPricingEngine(bsEngine);
116
117 Volatility impliedVol =
118 comp.impliedVolatility(price: npv, process: bsProcess, accuracy: 1e-10, maxEvaluations: 100);
119
120 if (std::fabs(x: impliedVol - expectedVol[i]) > tol) {
121 BOOST_FAIL("Failed to reproduce implied volatility"
122 << "\n calculated: " << impliedVol
123 << "\n expected : " << expectedVol[i]);
124 }
125 if (std::fabs(x: (comp.NPV() - npv)/npv) > tol) {
126 BOOST_FAIL("Failed to reproduce NPV"
127 << "\n calculated: " << npv
128 << "\n expected : " << comp.NPV());
129 }
130 if (std::fabs(x: comp.delta() - option.delta()) > tol) {
131 BOOST_FAIL("Failed to reproduce NPV"
132 << "\n calculated: " << npv
133 << "\n expected : " << comp.NPV());
134 }
135 if (std::fabs(x: (comp.gamma() - option.gamma())/npv) > tol) {
136 BOOST_FAIL("Failed to reproduce NPV"
137 << "\n calculated: " << npv
138 << "\n expected : " << comp.NPV());
139 }
140 if (std::fabs(x: (comp.theta() - option.theta())/npv) > tol) {
141 BOOST_FAIL("Failed to reproduce NPV"
142 << "\n calculated: " << npv
143 << "\n expected : " << comp.NPV());
144 }
145 if (std::fabs(x: (comp.vega() - option.vega())/npv) > tol) {
146 BOOST_FAIL("Failed to reproduce NPV"
147 << "\n calculated: " << npv
148 << "\n expected : " << comp.NPV());
149 }
150 }
151}
152
153void HybridHestonHullWhiteProcessTest::testCompareBsmHWandHestonHW() {
154 BOOST_TEST_MESSAGE("Comparing European option pricing for a BSM process"
155 " with one-factor Hull-White model...");
156
157 DayCounter dc = Actual365Fixed();
158
159 const Date today = Date::todaysDate();
160
161 Settings::instance().evaluationDate() = today;
162
163 const Handle<Quote> spot(
164 ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
165 std::vector<Date> dates;
166 std::vector<Rate> rates, divRates;
167
168 for (Size i=0; i <= 40; ++i) {
169 dates.push_back(x: today+Period(i, Years));
170 rates.push_back(x: 0.01 + 0.0002*std::exp(x: std::sin(x: i/4.0)));
171 divRates.push_back(x: 0.02 + 0.0001*std::exp(x: std::sin(x: i/5.0)));
172 }
173
174 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
175 const Handle<YieldTermStructure> rTS(
176 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
177 const Handle<YieldTermStructure> qTS(
178 ext::shared_ptr<YieldTermStructure>(
179 new ZeroCurve(dates, divRates, dc)));
180
181 ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.25));
182 const Handle<BlackVolTermStructure> volTS(flatVol(today, volatility: vol, dc));
183
184 ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
185 new BlackScholesMertonProcess(spot, qTS, rTS, volTS));
186
187 ext::shared_ptr<HestonProcess> hestonProcess(
188 new HestonProcess(rTS, qTS, spot,
189 vol->value()*vol->value(), 1.0,
190 vol->value()*vol->value(), 1e-4, 0.0));
191
192 ext::shared_ptr<HestonModel> hestonModel(new HestonModel(hestonProcess));
193
194 ext::shared_ptr<HullWhite> hullWhiteModel(
195 new HullWhite(Handle<YieldTermStructure>(rTS), 0.01, 0.01));
196
197 ext::shared_ptr<PricingEngine> bsmhwEngine(
198 new AnalyticBSMHullWhiteEngine(0.0, bsmProcess, hullWhiteModel));
199
200 ext::shared_ptr<PricingEngine> hestonHwEngine(
201 new AnalyticHestonHullWhiteEngine(hestonModel, hullWhiteModel, 128));
202
203
204 const Real tol = 1e-5;
205 const Real strike[] = { 0.25, 0.5, 0.75, 0.8, 0.9,
206 1.0, 1.1, 1.2, 1.5, 2.0, 4.0 };
207 const Size maturity[] = { 1, 2, 3, 5, 10, 15, 20, 25, 30 };
208 const Option::Type types[] = { Option::Put, Option::Call };
209
210 for (auto type : types) {
211 for (Real j : strike) {
212 for (unsigned long l : maturity) {
213 const Date maturityDate = today + Period(l, Years);
214
215 ext::shared_ptr<Exercise> exercise(
216 new EuropeanExercise(maturityDate));
217
218 Real fwd =
219 j * spot->value() * qTS->discount(d: maturityDate) / rTS->discount(d: maturityDate);
220
221 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(type, fwd));
222
223 EuropeanOption option(payoff, exercise);
224
225 option.setPricingEngine(bsmhwEngine);
226 const Real calculated = option.NPV();
227
228 option.setPricingEngine(hestonHwEngine);
229 const Real expected = option.NPV();
230
231 if (std::fabs(x: calculated-expected) > calculated*tol &&
232 std::fabs(x: calculated-expected) > tol) {
233 BOOST_ERROR("Failed to reproduce npvs"
234 << "\n calculated: " << calculated
235 << "\n expected : " << expected << "\n strike : " << j
236 << "\n maturity : " << l << "\n type : "
237 << ((type == Option::Put) ? "Put" : "Call"));
238 }
239 }
240 }
241 }
242}
243
244void HybridHestonHullWhiteProcessTest::testZeroBondPricing() {
245 BOOST_TEST_MESSAGE("Testing Monte-Carlo zero bond pricing...");
246
247 DayCounter dc = Actual360();
248 const Date today = Date::todaysDate();
249
250 Settings::instance().evaluationDate() = today;
251
252 // construct a strange yield curve to check drifts and discounting
253 // of the joint stochastic process
254
255 std::vector<Date> dates;
256 std::vector<Time> times;
257 std::vector<Rate> rates;
258
259 dates.push_back(x: today);
260 rates.push_back(x: 0.02);
261 times.push_back(x: 0.0);
262 for (Size i=120; i < 240; ++i) {
263 dates.push_back(x: today+Period(i, Months));
264 rates.push_back(x: 0.02 + 0.0002*std::exp(x: std::sin(x: i/8.0)));
265 times.push_back(x: dc.yearFraction(d1: today, d2: dates.back()));
266 }
267
268 const Date maturity = dates.back() + Period(10, Years);
269 dates.push_back(x: maturity);
270 rates.push_back(x: 0.04);
271 times.push_back(x: dc.yearFraction(d1: today, d2: dates.back()));
272
273 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
274
275 const Handle<YieldTermStructure> ts(
276 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
277 const Handle<YieldTermStructure> ds(flatRate(today, forward: 0.0, dc));
278
279 const ext::shared_ptr<HestonProcess> hestonProcess(
280 new HestonProcess(ts, ds, s0, 0.02, 1.0, 0.2, 0.5, -0.8));
281 const ext::shared_ptr<HullWhiteForwardProcess> hwProcess(
282 new HullWhiteForwardProcess(ts, 0.05, 0.05));
283 hwProcess->setForwardMeasureTime(dc.yearFraction(d1: today, d2: maturity));
284 const ext::shared_ptr<HullWhite> hwModel(new HullWhite(ts, 0.05, 0.05));
285
286 const ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess(
287 new HybridHestonHullWhiteProcess(hestonProcess, hwProcess, -0.4));
288
289 TimeGrid grid(times.begin(), times.end()-1);
290
291 typedef SobolBrownianBridgeRsg rsg_type;
292 typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
293
294 const Size factors = jointProcess->factors();
295 const Size steps = grid.size()-1;
296 rsg_type rsg = rsg_type(factors, steps);
297 MultiPathGenerator<rsg_type> generator(jointProcess, grid, rsg, false);
298
299 const Size m = 90;
300 std::vector<GeneralStatistics> zeroStat(m);
301 std::vector<GeneralStatistics> optionStat(m);
302
303 const Size nrTrails = 8191;
304 const Size optionTenor = 24;
305 const DiscountFactor strike = 0.5;
306
307 for (Size i=0; i < nrTrails; ++i) {
308 sample_type path = generator.next();
309
310 for (Size j=1; j < m; ++j) {
311 const Time t = grid[j]; // zero end and option maturity
312 const Time T = grid[j+optionTenor];// maturity of zero bond
313 // of option
314
315 Array states(3);
316 Array optionStates(3);
317 for (Size k=0; k < jointProcess->size(); ++k) {
318 states[k] = path.value[k][j];
319 optionStates[k] = path.value[k][j+optionTenor];
320 }
321
322 const DiscountFactor zeroBond
323 = 1.0/jointProcess->numeraire(t, x: states);
324 const DiscountFactor zeroOption = zeroBond
325 * std::max(a: 0.0, b: hwModel->discountBond(now: t, maturity: T, rate: states[2])-strike);
326
327 zeroStat[j].add(value: zeroBond);
328 optionStat[j].add(value: zeroOption);
329 }
330 }
331
332 for (Size j=1; j < m; ++j) {
333 const Time t = grid[j];
334 Real calculated = zeroStat[j].mean();
335 Real expected = ts->discount(t);
336
337 if (std::fabs(x: calculated - expected) > 0.03) {
338 BOOST_ERROR("Failed to reproduce expected zero bond prices"
339 << "\n t: " << t
340 << "\n calculated: " << calculated
341 << "\n expected: " << expected);
342 }
343
344 const Time T = grid[j+optionTenor];
345
346 calculated = optionStat[j].mean();
347 expected = hwModel->discountBondOption(type: Option::Call, strike, maturity: t, bondMaturity: T);
348
349 if (std::fabs(x: calculated - expected) > 0.0035) {
350 BOOST_ERROR("Failed to reproduce expected zero bond option prices"
351 << "\n t: " << t
352 << "\n T: " << T
353 << "\n calculated: " << calculated
354 << "\n expected: " << expected);
355 }
356 }
357}
358
359void HybridHestonHullWhiteProcessTest::testMcVanillaPricing() {
360 BOOST_TEST_MESSAGE("Testing Monte-Carlo vanilla option pricing...");
361
362 DayCounter dc = Actual360();
363 const Date today = Date::todaysDate();
364
365 Settings::instance().evaluationDate() = today;
366
367 // construct a strange yield curve to check drifts and discounting
368 // of the joint stochastic process
369
370 std::vector<Date> dates;
371 std::vector<Rate> rates, divRates;
372
373 for (Size i=0; i <= 40; ++i) {
374 dates.push_back(x: today+Period(i, Years));
375 rates.push_back(x: 0.03 + 0.0003*std::exp(x: std::sin(x: i/4.0)));
376 divRates.push_back(x: 0.02 + 0.0001*std::exp(x: std::sin(x: i/5.0)));
377 }
378
379 const Date maturity = today + Period(20, Years);
380
381 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
382 const Handle<YieldTermStructure> rTS(
383 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
384 const Handle<YieldTermStructure> qTS(
385 ext::shared_ptr<YieldTermStructure>(
386 new ZeroCurve(dates, divRates, dc)));
387 ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.25));
388 const Handle<BlackVolTermStructure> volTS(flatVol(today, volatility: vol, dc));
389
390 const ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
391 new BlackScholesMertonProcess(s0, qTS, rTS, volTS));
392 const ext::shared_ptr<HestonProcess> hestonProcess(
393 new HestonProcess(rTS, qTS, s0, 0.0625, 0.5, 0.0625, 1e-5, 0.3));
394 const ext::shared_ptr<HullWhiteForwardProcess> hwProcess(
395 new HullWhiteForwardProcess(rTS, 0.01, 0.01));
396 hwProcess->setForwardMeasureTime(dc.yearFraction(d1: today, d2: maturity));
397
398 const Real tol = 0.05;
399 const Real corr[] = {-0.9, -0.5, 0.0, 0.5, 0.9 };
400 const Real strike[] = { 100 };
401
402 for (Real i : corr) {
403 for (Real j : strike) {
404 ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess(
405 new HybridHestonHullWhiteProcess(hestonProcess, hwProcess, i));
406
407 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(Option::Put, j));
408 ext::shared_ptr<Exercise> exercise(
409 new EuropeanExercise(maturity));
410
411 VanillaOption optionHestonHW(payoff, exercise);
412 ext::shared_ptr<PricingEngine> engine =
413 MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
414 .withSteps(steps: 1)
415 .withAntitheticVariate()
416 .withControlVariate()
417 .withAbsoluteTolerance(tolerance: tol)
418 .withSeed(seed: 42);
419
420 optionHestonHW.setPricingEngine(engine);
421
422 const ext::shared_ptr<HullWhite> hwModel(
423 new HullWhite(Handle<YieldTermStructure>(rTS),
424 hwProcess->a(), hwProcess->sigma()));
425
426 VanillaOption optionBsmHW(payoff, exercise);
427 optionBsmHW.setPricingEngine(ext::shared_ptr<PricingEngine>(
428 new AnalyticBSMHullWhiteEngine(i, bsmProcess, hwModel)));
429
430 const Real calculated = optionHestonHW.NPV();
431 const Real error = optionHestonHW.errorEstimate();
432 const Real expected = optionBsmHW.NPV();
433
434 if ((i != 0.0 && std::fabs(x: calculated - expected) > 3 * error) ||
435 (i == 0.0 && std::fabs(x: calculated - expected) > 1e-4)) {
436 BOOST_ERROR("Failed to reproduce BSM-HW vanilla prices"
437 << "\n corr: " << i << "\n strike: " << j
438 << "\n calculated: " << calculated << "\n error: " << error
439 << "\n expected: " << expected);
440 }
441 }
442 }
443}
444
445
446void HybridHestonHullWhiteProcessTest::testMcPureHestonPricing() {
447 BOOST_TEST_MESSAGE("Testing Monte-Carlo Heston option pricing...");
448
449 DayCounter dc = Actual360();
450 const Date today = Date::todaysDate();
451
452 Settings::instance().evaluationDate() = today;
453
454 // construct a strange yield curve to check drifts and discounting
455 // of the joint stochastic process
456
457 std::vector<Date> dates;
458 std::vector<Rate> rates, divRates;
459
460 for (Size i=0; i <= 100; ++i) {
461 dates.push_back(x: today+Period(i, Months));
462 rates.push_back(x: 0.02 + 0.0002*std::exp(x: std::sin(x: i/10.0)));
463 divRates.push_back(x: 0.02 + 0.0001*std::exp(x: std::sin(x: i/20.0)));
464 }
465
466 const Date maturity = today + Period(2, Years);
467
468 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
469 const Handle<YieldTermStructure> rTS(
470 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
471 const Handle<YieldTermStructure> qTS(
472 ext::shared_ptr<YieldTermStructure>(
473 new ZeroCurve(dates, divRates, dc)));
474
475 const ext::shared_ptr<HestonProcess> hestonProcess(
476 new HestonProcess(rTS, qTS, s0, 0.08, 1.5, 0.0625, 0.5, -0.8));
477 const ext::shared_ptr<HullWhiteForwardProcess> hwProcess(
478 new HullWhiteForwardProcess(rTS, 0.1, 1e-8));
479 hwProcess->setForwardMeasureTime(dc.yearFraction(
480 d1: today, d2: maturity+Period(1, Years)));
481
482 const Real tol = 0.001;
483 const Real corr[] = { -0.45, 0.45, 0.25 };
484 const Real strike[] = { 100, 75, 50, 150 };
485
486 for (Real i : corr) {
487 for (Real j : strike) {
488 ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess(
489 new HybridHestonHullWhiteProcess(hestonProcess, hwProcess, i,
490 HybridHestonHullWhiteProcess::Euler));
491
492 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(Option::Put, j));
493 ext::shared_ptr<Exercise> exercise(
494 new EuropeanExercise(maturity));
495
496 VanillaOption optionHestonHW(payoff, exercise);
497 VanillaOption optionPureHeston(payoff, exercise);
498 optionPureHeston.setPricingEngine(
499 ext::shared_ptr<PricingEngine>(
500 new AnalyticHestonEngine(
501 ext::make_shared<HestonModel>(
502 args: hestonProcess))));
503
504 Real expected = optionPureHeston.NPV();
505
506 optionHestonHW.setPricingEngine(
507 MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
508 .withSteps(steps: 2)
509 .withAntitheticVariate()
510 .withControlVariate()
511 .withAbsoluteTolerance(tolerance: tol)
512 .withSeed(seed: 42));
513
514 Real calculated = optionHestonHW.NPV();
515 Real error = optionHestonHW.errorEstimate();
516
517 if ( std::fabs(x: calculated - expected) > 3*error
518 && std::fabs(x: calculated - expected) > tol) {
519 BOOST_ERROR("Failed to reproduce pure heston vanilla prices"
520 << "\n corr: " << i << "\n strike: " << j
521 << "\n calculated: " << calculated << "\n error: " << error
522 << "\n expected: " << expected);
523 }
524 }
525 }
526}
527
528
529void HybridHestonHullWhiteProcessTest::testAnalyticHestonHullWhitePricing() {
530 BOOST_TEST_MESSAGE("Testing analytic Heston Hull-White option pricing...");
531
532 DayCounter dc = Actual360();
533 const Date today = Date::todaysDate();
534
535 Settings::instance().evaluationDate() = today;
536
537 // construct a strange yield curve to check drifts and discounting
538 // of the joint stochastic process
539
540 std::vector<Date> dates;
541 std::vector<Rate> rates, divRates;
542
543 for (Size i=0; i <= 40; ++i) {
544 dates.push_back(x: today+Period(i, Years));
545 rates.push_back(x: 0.03 + 0.0001*std::exp(x: std::sin(x: i/4.0)));
546 divRates.push_back(x: 0.02 + 0.0002*std::exp(x: std::sin(x: i/3.0)));
547 }
548
549 const Date maturity = today + Period(5, Years);
550 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
551 const Handle<YieldTermStructure> rTS(
552 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
553 const Handle<YieldTermStructure> qTS(
554 ext::shared_ptr<YieldTermStructure>(
555 new ZeroCurve(dates, divRates, dc)));
556
557 const ext::shared_ptr<HestonProcess> hestonProcess(
558 new HestonProcess(rTS, qTS, s0, 0.08, 1.5, 0.0625, 0.5, -0.8));
559 const ext::shared_ptr<HestonModel> hestonModel(
560 new HestonModel(hestonProcess));
561
562 const ext::shared_ptr<HullWhiteForwardProcess> hwFwdProcess(
563 new HullWhiteForwardProcess(rTS, 0.01, 0.01));
564 hwFwdProcess->setForwardMeasureTime(dc.yearFraction(d1: today, d2: maturity));
565 const ext::shared_ptr<HullWhite> hullWhiteModel(new HullWhite(
566 rTS, hwFwdProcess->a(), hwFwdProcess->sigma()));
567
568 const Real tol = 0.002;
569 const Real strike[] = { 80, 120 };
570 const Option::Type types[] = { Option::Put, Option::Call };
571
572 for (auto type : types) {
573 for (Real j : strike) {
574 ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess(
575 new HybridHestonHullWhiteProcess(
576 hestonProcess, hwFwdProcess, 0.0,
577 HybridHestonHullWhiteProcess::Euler));
578
579 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(type, j));
580 ext::shared_ptr<Exercise> exercise(
581 new EuropeanExercise(maturity));
582
583 VanillaOption optionHestonHW(payoff, exercise);
584 optionHestonHW.setPricingEngine(
585 MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
586 .withSteps(steps: 1)
587 .withAntitheticVariate()
588 .withControlVariate()
589 .withAbsoluteTolerance(tolerance: tol)
590 .withSeed(seed: 42));
591
592 VanillaOption optionPureHeston(payoff, exercise);
593 optionPureHeston.setPricingEngine(
594 ext::shared_ptr<PricingEngine>(
595 new AnalyticHestonHullWhiteEngine(hestonModel,
596 hullWhiteModel, 128)));
597
598 Real calculated = optionHestonHW.NPV();
599 Real error = optionHestonHW.errorEstimate();
600 Real expected = optionPureHeston.NPV();
601
602 if ( std::fabs(x: calculated - expected) > 3*error
603 && std::fabs(x: calculated - expected) > tol) {
604 BOOST_ERROR("Failed to reproduce hw heston vanilla prices"
605 << "\n strike: " << j << "\n calculated: " << calculated
606 << "\n error: " << error << "\n expected: " << expected);
607 }
608 }
609 }
610}
611
612void HybridHestonHullWhiteProcessTest::testCallableEquityPricing() {
613 BOOST_TEST_MESSAGE("Testing the pricing of a callable equity product...");
614
615 /*
616 For the definition of the example product see
617 Alexander Giese, On the Pricing of Auto-Callable Equity
618 Structures in the Presence of Stochastic Volatility and
619 Stochastic Interest Rates .
620 http://workshop.mathfinance.de/2006/papers/giese/slides.pdf
621 */
622
623 const Size maturity = 7;
624 DayCounter dc = Actual365Fixed();
625 const Date today = Date::todaysDate();
626
627 Settings::instance().evaluationDate() = today;
628
629 Handle<Quote> spot(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
630 ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.04));
631 Handle<YieldTermStructure> qTS(flatRate(today, forward: qRate, dc));
632 ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.04));
633 Handle<YieldTermStructure> rTS(flatRate(today, forward: rRate, dc));
634
635 const ext::shared_ptr<HestonProcess> hestonProcess(
636 new HestonProcess(rTS, qTS, spot, 0.0625, 1.0,
637 0.24*0.24, 1e-4, 0.0));
638 const ext::shared_ptr<HullWhiteForwardProcess> hwProcess(
639 new HullWhiteForwardProcess(rTS, 0.00883, 0.00526));
640 hwProcess->setForwardMeasureTime(
641 dc.yearFraction(d1: today, d2: today+Period(maturity+1, Years)));
642
643 const ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess(
644 new HybridHestonHullWhiteProcess(hestonProcess, hwProcess, -0.4));
645
646 Schedule schedule(today, today + Period(maturity, Years),
647 Period(1, Years), TARGET(),
648 Following, Following,
649 DateGeneration::Forward, false);
650
651 std::vector<Time> times(maturity+1);
652 std::transform(first: schedule.begin(), last: schedule.end(), result: times.begin(),
653 unary_op: [&](const Date& d) { return dc.yearFraction(d1: today, d2: d); });
654
655 for (Size i=0; i<=maturity; ++i)
656 times[i] = static_cast<Time>(i);
657
658 TimeGrid grid(times.begin(), times.end());
659
660 std::vector<Real> redemption(maturity);
661 for (Size i=0; i < maturity; ++i) {
662 redemption[i] = 1.07 + 0.03*i;
663 }
664
665 typedef PseudoRandom::rsg_type rsg_type;
666 typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
667
668 BigNatural seed = 42;
669 rsg_type rsg = PseudoRandom::make_sequence_generator(
670 dimension: jointProcess->factors()*(grid.size()-1), seed);
671
672 MultiPathGenerator<rsg_type> generator(jointProcess, grid, rsg, false);
673 GeneralStatistics stat;
674
675 Real antitheticPayoff=0;
676 const Size nrTrails = 40000;
677 for (Size i=0; i < nrTrails; ++i) {
678 const bool antithetic = (i % 2) != 0;
679
680 sample_type path = antithetic ? generator.antithetic()
681 : generator.next();
682
683 Real payoff=0;
684 for (Size j=1; j <= maturity; ++j) {
685 if (path.value[0][j] > spot->value()) {
686 Array states(3);
687 for (Size k=0; k < 3; ++k) {
688 states[k] = path.value[k][j];
689 }
690 payoff = redemption[j-1]
691 / jointProcess->numeraire(t: grid[j], x: states);
692 break;
693 }
694 else if (j == maturity) {
695 Array states(3);
696 for (Size k=0; k < 3; ++k) {
697 states[k] = path.value[k][j];
698 }
699 payoff = 1.0 / jointProcess->numeraire(t: grid[j], x: states);
700 }
701 }
702
703 if (antithetic){
704 stat.add(value: 0.5*(antitheticPayoff + payoff));
705 }
706 else {
707 antitheticPayoff = payoff;
708 }
709 }
710
711 const Real expected = 0.938;
712 const Real calculated = stat.mean();
713 const Real error = stat.errorEstimate();
714
715 if (std::fabs(x: expected - calculated) > 3*error) {
716 BOOST_ERROR("Failed to reproduce auto-callable equity structure price"
717 << "\n calculated: " << calculated
718 << "\n error: " << error
719 << "\n expected: " << expected);
720 }
721}
722
723void HybridHestonHullWhiteProcessTest::testDiscretizationError() {
724 BOOST_TEST_MESSAGE("Testing the discretization error of the "
725 "Heston Hull-White process...");
726
727 DayCounter dc = Actual360();
728 const Date today = Date::todaysDate();
729
730 Settings::instance().evaluationDate() = today;
731
732 // construct a strange yield curve to check drifts and discounting
733 // of the joint stochastic process
734
735 std::vector<Date> dates;
736 std::vector<Rate> rates, divRates;
737
738 for (Size i=0; i <= 31; ++i) {
739 dates.push_back(x: today+Period(i, Years));
740 rates.push_back(x: 0.04 + 0.0001*std::exp(x: std::sin(x: double(i))));
741 divRates.push_back(x: 0.04 + 0.0001*std::exp(x: std::sin(x: double(i))));
742 }
743
744 const Date maturity = today + Period(10, Years);
745 const Volatility v = 0.25;
746
747 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100)));
748 const ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(v));
749 const Handle<BlackVolTermStructure> volTS(flatVol(today, volatility: vol, dc));
750 const Handle<YieldTermStructure> rTS(
751 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
752 const Handle<YieldTermStructure> qTS(
753 ext::shared_ptr<YieldTermStructure>(
754 new ZeroCurve(dates, divRates, dc)));
755
756 const ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
757 new BlackScholesMertonProcess(s0, qTS, rTS, volTS));
758
759 const ext::shared_ptr<HestonProcess> hestonProcess(
760 new HestonProcess(rTS, qTS, s0, v*v, 1, v*v, 1e-6, -0.4));
761
762 const ext::shared_ptr<HullWhiteForwardProcess> hwProcess(
763 new HullWhiteForwardProcess(rTS, 0.01, 0.01));
764 hwProcess->setForwardMeasureTime(20.1472222222222222);
765
766 const Real tol = 0.05;
767 const Real corr[] = {-0.85, 0.5 };
768 const Real strike[] = { 50, 100, 125 };
769
770 for (Real i : corr) {
771 for (Real j : strike) {
772 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(Option::Put, j));
773 ext::shared_ptr<Exercise> exercise(
774 new EuropeanExercise(maturity));
775
776 VanillaOption optionBsmHW(payoff, exercise);
777 const ext::shared_ptr<HullWhite> hwModel(new HullWhite(
778 rTS, hwProcess->a(), hwProcess->sigma()));
779 optionBsmHW.setPricingEngine(ext::shared_ptr<PricingEngine>(
780 new AnalyticBSMHullWhiteEngine(i, bsmProcess, hwModel)));
781
782 Real expected = optionBsmHW.NPV();
783
784 VanillaOption optionHestonHW(payoff, exercise);
785 ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess(
786 new HybridHestonHullWhiteProcess(hestonProcess, hwProcess, i));
787 optionHestonHW.setPricingEngine(
788 MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
789 .withSteps(steps: 1)
790 .withAntitheticVariate()
791 .withAbsoluteTolerance(tolerance: tol)
792 .withSeed(seed: 42));
793
794 Real calculated = optionHestonHW.NPV();
795 Real error = optionHestonHW.errorEstimate();
796
797 if (( std::fabs(x: calculated - expected) > 3*error
798 && std::fabs(x: calculated - expected) > 1e-5)) {
799 BOOST_ERROR("Failed to reproduce discretization error"
800 << "\n corr: " << i << "\n strike: " << j
801 << "\n calculated: " << calculated << "\n error: " << error
802 << "\n expected: " << expected);
803 }
804 }
805 }
806}
807
808void HybridHestonHullWhiteProcessTest::testFdmHestonHullWhiteEngine() {
809 BOOST_TEST_MESSAGE("Testing the FDM Heston Hull-White engine...");
810
811 const Date today = Date(28, March, 2004);
812 Settings::instance().evaluationDate() = today;
813 const Date exerciseDate = Date(28, March, 2012);
814 DayCounter dc = Actual365Fixed();
815
816 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
817
818 const Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc));
819 const Handle<YieldTermStructure> qTS(flatRate(forward: 0.02, dc));
820
821 const Volatility vol = 0.30;
822 const Handle<BlackVolTermStructure> volTS(flatVol(volatility: vol, dc));
823
824 const Real v0 = vol*vol;
825 ext::shared_ptr<HestonProcess> hestonProcess(
826 new HestonProcess(rTS, qTS, s0, v0, 1.0, v0, 0.000001, 0.0));
827
828 ext::shared_ptr<BlackScholesMertonProcess> stochProcess(
829 new BlackScholesMertonProcess(s0, qTS, rTS, volTS));
830
831 ext::shared_ptr<HullWhiteProcess> hwProcess(
832 new HullWhiteProcess(rTS, 0.00883, 0.01));
833 ext::shared_ptr<HullWhite> hwModel(
834 new HullWhite(rTS, hwProcess->a(), hwProcess->sigma()));
835
836 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exerciseDate));
837 const Real corr[] = {-0.85, 0.5 };
838 const Real strike[] = { 75, 120, 160 };
839
840 for (Real i : corr) {
841 for (Real j : strike) {
842 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(Option::Call, j));
843 VanillaOption option(payoff, exercise);
844
845 option.setPricingEngine(
846 ext::shared_ptr<PricingEngine>(new FdHestonHullWhiteVanillaEngine(
847 ext::make_shared<HestonModel>(args&: hestonProcess), hwProcess, i, 50, 200, 10, 15)));
848 const Real calculated = option.NPV();
849 const Real calculatedDelta = option.delta();
850 const Real calculatedGamma = option.gamma();
851
852 option.setPricingEngine(ext::shared_ptr<PricingEngine>(
853 new AnalyticBSMHullWhiteEngine(i, stochProcess, hwModel)));
854 const Real expected = option.NPV();
855 const Real expectedDelta = option.delta();
856 const Real expectedGamma = option.gamma();
857
858 const Real npvTol = 0.01;
859 if (std::fabs(x: calculated - expected) > npvTol) {
860 BOOST_ERROR("Failed to reproduce analytic npv values"
861 << "\n corr: " << i << "\n strike: " << j
862 << "\n calculated: " << calculated
863 << "\n expected: " << expected);
864 }
865 const Real deltaTol = 0.001;
866 if (std::fabs(x: calculatedDelta - expectedDelta) > deltaTol) {
867 BOOST_ERROR("Failed to reproduce analytic delta values"
868 << "\n corr: " << i << "\n strike: " << j
869 << "\n calculated: " << calculated
870 << "\n expected: " << expected);
871 }
872 const Real gammaTol = 0.001;
873 if (std::fabs(x: calculatedGamma - expectedGamma) > gammaTol) {
874 BOOST_ERROR("Failed to reproduce analytic gamma values"
875 << "\n corr: " << i << "\n strike: " << j
876 << "\n calculated: " << calculated
877 << "\n expected: " << expected);
878 }
879 }
880 }
881}
882
883
884namespace hybrid_heston_hullwhite_process_test {
885
886 struct HestonModelData {
887 const char* const name;
888 Real v0;
889 Real kappa;
890 Real theta;
891 Real sigma;
892 Real rho;
893 Real r;
894 Real q;
895 };
896
897 HestonModelData hestonModels[] = {
898 // ADI finite difference schemes for option pricing in the
899 // Heston model with correlation, K.J. in t'Hout and S. Foulon,
900 {.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},
901 {.name: "'t Hout case 2", .v0: 0.12, .kappa: 3.0, .theta: 0.12, .sigma: 0.04, .rho: 0.6, .r: 0.01, .q: 0.04},
902 {.name: "'t Hout case 3", .v0: 0.0707,.kappa: 0.6067, .theta: 0.0707, .sigma: 0.2928, .rho: -0.7571, .r: 0.03, .q: 0.0},
903 {.name: "'t Hout case 4", .v0: 0.06, .kappa: 2.5, .theta: 0.06, .sigma: 0.5, .rho: -0.1, .r: 0.0507, .q: 0.0469},
904 // Efficient numerical methods for pricing American options under
905 // stochastic volatility, Samuli Ikonen and Jari Toivanen,
906 {.name: "Ikonen-Toivanen", .v0: 0.0625, .kappa: 5, .theta: 0.16, .sigma: 0.9, .rho: 0.1, .r: 0.1, .q: 0.0},
907 // Not-so-complex logarithms in the Heston model,
908 // Christian Kahl and Peter Jäckel
909 {.name: "Kahl-Jaeckel", .v0: 0.16, .kappa: 1.0, .theta: 0.16, .sigma: 2.0, .rho: -0.8, .r: 0.0, .q: 0.0},
910 // self defined test cases
911 {.name: "Equity case", .v0: 0.07, .kappa: 2.0, .theta: 0.04, .sigma: 0.55, .rho: -0.8, .r: 0.03, .q: 0.035 },
912 {.name: "high correlation", .v0: 0.07, .kappa: 1.0, .theta: 0.04, .sigma: 0.55, .rho: 0.995, .r: 0.02, .q: 0.04 },
913 {.name: "low Vol-Of-Vol", .v0: 0.07, .kappa: 1.0, .theta: 0.04, .sigma: 0.001, .rho: -0.75, .r: 0.04, .q: 0.03 },
914 {.name: "kappaEqSigRho", .v0: 0.07, .kappa: 0.4, .theta: 0.04, .sigma: 0.5, .rho: 0.8, .r: 0.03, .q: 0.03 }
915 };
916
917 struct HullWhiteModelData {
918 const char* const name;
919 Real a;
920 Real sigma;
921 };
922
923 HullWhiteModelData hullWhiteModels[] = {
924 {.name: "EUR-2003", .a: 0.00883, .sigma: 0.00631 }
925 };
926
927
928 struct SchemeData {
929 const char* const name;
930 FdmSchemeDesc schemeDesc;
931 };
932
933 SchemeData schemes[] = {
934 { .name: "HV2", .schemeDesc: FdmSchemeDesc::Hundsdorfer() },
935 { .name: "HV1", .schemeDesc: FdmSchemeDesc::ModifiedHundsdorfer() },
936 { .name: "CS" , .schemeDesc: FdmSchemeDesc::CraigSneyd() },
937 { .name: "MCS", .schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd() },
938 { .name: "DS" , .schemeDesc: FdmSchemeDesc::Douglas() }
939 };
940
941 struct VanillaOptionData {
942 Real strike;
943 Time maturity;
944 Option::Type optionType;
945 };
946
947 ext::shared_ptr<HestonProcess> makeHestonProcess(
948 const HestonModelData& params) {
949
950 Handle<Quote> spot(
951 ext::make_shared<SimpleQuote>(args: 100));
952
953 DayCounter dayCounter = Actual365Fixed();
954 Handle<YieldTermStructure> rTS(flatRate(forward: params.r, dc: dayCounter));
955 Handle<YieldTermStructure> qTS(flatRate(forward: params.q, dc: dayCounter));
956
957 return ext::make_shared<HestonProcess>(
958 args&: rTS, args&: qTS, args&: spot, args: params.v0, args: params.kappa,
959 args: params.theta, args: params.sigma, args: params.rho);
960 }
961
962 ext::shared_ptr<VanillaOption> makeVanillaOption(
963 const VanillaOptionData& params) {
964
965 Date maturity = Date(Settings::instance().evaluationDate())
966 + Period(Size(params.maturity*365), Days);
967 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(maturity));
968 ext::shared_ptr<StrikedTypePayoff> payoff(
969 new PlainVanillaPayoff(params.optionType, params.strike));
970
971 return ext::make_shared<VanillaOption>(
972 args&: payoff, args&: exercise);
973 }
974}
975
976void HybridHestonHullWhiteProcessTest::testBsmHullWhitePricing() {
977 BOOST_TEST_MESSAGE("Testing convergence speed of Heston-Hull-White engine...");
978
979 using namespace hybrid_heston_hullwhite_process_test;
980
981 Date today(27, December, 2004);
982 Settings::instance().evaluationDate() = today;
983
984
985 Real maturity = 5.0;
986 Real equityIrCorr = -0.4;
987 std::vector<Real> strikes = {75,85,90,95,100,105,110,115,120,125,130,140,150};
988 Size listOfTimeStepsPerYear[] = { 20 };
989
990 HestonModelData hestonModelData
991 = { .name: "BSM-HW Model", .v0: 0.09, .kappa: 1.0, .theta: 0.09, QL_EPSILON, .rho: 0.0, .r: 0.04, .q: 0.03 };
992 HullWhiteModelData hwModelData = hullWhiteModels[0];
993 bool controlVariate[] = { true, false };
994
995 ext::shared_ptr<HestonProcess> hp(makeHestonProcess(params: hestonModelData));
996 ext::shared_ptr<HestonModel> hestonModel(new HestonModel(hp));
997
998 ext::shared_ptr<HullWhiteProcess> hwProcess(
999 new HullWhiteProcess(hp->riskFreeRate(),
1000 hwModelData.a, hwModelData.sigma));
1001 ext::shared_ptr<HullWhite> hullWhiteModel(
1002 new HullWhite(hp->riskFreeRate(),
1003 hwProcess->a(), hwProcess->sigma()));
1004
1005
1006 ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
1007 new BlackScholesMertonProcess(
1008 hp->s0(), hp->dividendYield(), hp->riskFreeRate(),
1009 Handle<BlackVolTermStructure>(
1010 flatVol(today, volatility: std::sqrt(x: hestonModelData.theta),
1011 dc: hp->riskFreeRate()->dayCounter()))));
1012
1013 ext::shared_ptr<PricingEngine> bsmhwEngine(
1014 new AnalyticBSMHullWhiteEngine(equityIrCorr, bsmProcess,
1015 hullWhiteModel));
1016
1017 Real tolWithCV[] = { 2e-4, 2e-4, 2e-4, 2e-4, 0.01 };
1018 Real tolWithOutCV[] = { 5e-3, 5e-3, 5e-3, 5e-3, 0.02 };
1019 for (Size l=0; l < LENGTH(schemes); ++l) {
1020 SchemeData scheme = schemes[l];
1021 for (bool i : controlVariate) {
1022 for (unsigned long u : listOfTimeStepsPerYear) {
1023 Size tSteps = Size(maturity * u);
1024
1025 ext::shared_ptr<FdHestonHullWhiteVanillaEngine> fdEngine(
1026 new FdHestonHullWhiteVanillaEngine(hestonModel, hwProcess, equityIrCorr, tSteps,
1027 400, 2, 10, 0, i, scheme.schemeDesc));
1028 fdEngine->enableMultipleStrikesCaching(strikes);
1029
1030 Real avgPriceDiff = 0.0;
1031 for (Real& strike : strikes) {
1032 VanillaOptionData optionData = {.strike: strike, .maturity: maturity, .optionType: Option::Call};
1033 ext::shared_ptr<VanillaOption> option
1034 = makeVanillaOption(params: optionData);
1035 option->setPricingEngine(bsmhwEngine);
1036 Real expected = option->NPV();
1037
1038 option->setPricingEngine(fdEngine);
1039 Real calculated = option->NPV();
1040 avgPriceDiff
1041 += std::fabs(x: expected-calculated)/strikes.size(); // NOLINT(bugprone-integer-division)
1042 }
1043
1044 if (i && tolWithCV[l] < avgPriceDiff) {
1045 BOOST_ERROR("Failed to reproduce BSM-Hull-White prices"
1046 << "\n scheme : " << scheme.name << "\n model : "
1047 << hestonModelData.name << "\n CV : on");
1048 }
1049
1050
1051 if (!i && tolWithOutCV[l] < avgPriceDiff) {
1052 BOOST_ERROR("Failed to reproduce BSM-Hull-White prices"
1053 << "\n scheme : " << scheme.name
1054 << "\n model : " << hestonModelData.name
1055 << "\n CV : off");
1056 }
1057 }
1058 }
1059 }
1060}
1061
1062void HybridHestonHullWhiteProcessTest::testSpatialDiscretizatinError() {
1063 BOOST_TEST_MESSAGE("Testing spatial convergence speed of Heston engine...");
1064
1065 using namespace hybrid_heston_hullwhite_process_test;
1066
1067 Date today(27, December, 2004);
1068 Settings::instance().evaluationDate() = today;
1069
1070 Real maturity=1.0;
1071 std::vector<Real> strikes = {75,85,90,95,100,105,110,115,120,125,130,140,150};
1072 Size listOfTimeStepsPerYear[] = { 40 };
1073
1074 const Real tol[] = { 0.02, 0.02, 0.02, 0.02, 0.05 };
1075 for (unsigned long u : listOfTimeStepsPerYear) {
1076 for (Size i=0; i < LENGTH(schemes); ++i) {
1077 for (auto& j : hestonModels) {
1078 Real avgPriceDiff = 0;
1079 ext::shared_ptr<HestonProcess> hestonProcess(makeHestonProcess(params: j));
1080 ext::shared_ptr<HestonModel> hestonModel(
1081 new HestonModel(hestonProcess));
1082
1083 ext::shared_ptr<PricingEngine> analyticEngine(
1084 new AnalyticHestonEngine(hestonModel, 172));
1085
1086 Size tSteps = Size(maturity * u);
1087
1088 ext::shared_ptr<FdHestonVanillaEngine> fdEngine(
1089 new FdHestonVanillaEngine(
1090 hestonModel, tSteps, 200, 40, 0,
1091 schemes[i].schemeDesc));
1092 fdEngine->enableMultipleStrikesCaching(strikes);
1093
1094 for (Real& strike : strikes) {
1095 VanillaOptionData optionData = {.strike: strike, .maturity: maturity, .optionType: Option::Call};
1096 ext::shared_ptr<VanillaOption> option
1097 = makeVanillaOption(params: optionData);
1098 option->setPricingEngine(analyticEngine);
1099 Real expected = option->NPV();
1100
1101 option->setPricingEngine(fdEngine);
1102 Real calculated = option->NPV();
1103
1104 avgPriceDiff
1105 += std::fabs(x: expected-calculated)/strikes.size(); // NOLINT(bugprone-integer-division)
1106 }
1107
1108 if (avgPriceDiff > tol[i]) {
1109 BOOST_ERROR("\nFailed to reproduce Heston prices"
1110 << "\n scheme : " << schemes[i].name
1111 << "\n model : " << j.name << "\n error : " << avgPriceDiff
1112 << "\n tolerance : " << tol[i]);
1113 }
1114 }
1115 }
1116 }
1117}
1118
1119
1120
1121
1122namespace hybrid_heston_hullwhite_process_test {
1123 class HestonHullWhiteCorrelationConstraint : public Constraint {
1124 private:
1125 class Impl : public Constraint::Impl {
1126 public:
1127 explicit Impl(Real equityShortRateCorr)
1128 : equityShortRateCorr_(equityShortRateCorr) {}
1129
1130 bool test(const Array& params) const override {
1131 const Real rho = params[3];
1132
1133 return (squared(x: rho) + squared(x: equityShortRateCorr_) <= 1.0);
1134 }
1135
1136 private:
1137 const Real equityShortRateCorr_;
1138 };
1139 public:
1140 explicit HestonHullWhiteCorrelationConstraint(
1141 Real equityShortRateCorr)
1142 : Constraint(ext::shared_ptr<Constraint::Impl>(
1143 new HestonHullWhiteCorrelationConstraint::Impl(
1144 equityShortRateCorr))) {}
1145 };
1146}
1147
1148
1149void HybridHestonHullWhiteProcessTest::testHestonHullWhiteCalibration() {
1150 BOOST_TEST_MESSAGE("Testing the Heston Hull-White calibration...");
1151
1152 using namespace hybrid_heston_hullwhite_process_test;
1153
1154 // Calibration of a hybrid Heston-Hull-White model using
1155 // the finite difference HestonHullWhite pricing engine
1156 //
1157 // Input surface is based on a Heston-Hull-White model with
1158 // Hull-White: a = 0.00883, \sigma = 0.00631
1159 // Heston : \nu = 0.12, \kappa = 2.0,
1160 // \theta = 0.09, \sigma = 0.5, \rho=-0.75
1161 // Equity Short rate correlation: -0.5
1162
1163 const DayCounter dc = Actual365Fixed();
1164 const Calendar calendar = TARGET();
1165 const Date today = Date(28, March, 2004);
1166 Settings::instance().evaluationDate() = today;
1167
1168 const Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc));
1169
1170 // assuming, that the Hull-White process is already calibrated
1171 // on a given set of pure interest rate calibration instruments.
1172 ext::shared_ptr<HullWhiteProcess> hwProcess(
1173 new HullWhiteProcess(rTS, 0.00883, 0.00631));
1174 ext::shared_ptr<HullWhite> hullWhiteModel(
1175 new HullWhite(rTS, hwProcess->a(), hwProcess->sigma()));
1176
1177 const Handle<YieldTermStructure> qTS(flatRate(forward: 0.02, dc));
1178 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
1179
1180 // starting point of the pure Heston calibration
1181 const Real start_v0 = 0.2*0.2;
1182 const Real start_theta = start_v0;
1183 const Real start_kappa = 0.5;
1184 const Real start_sigma = 0.25;
1185 const Real start_rho = -0.5;
1186
1187 const ext::shared_ptr<HestonProcess> hestonProcess(
1188 new HestonProcess(rTS, qTS, s0, start_v0, start_kappa,
1189 start_theta, start_sigma, start_rho));
1190 const ext::shared_ptr<HestonModel> analyticHestonModel
1191 (new HestonModel(hestonProcess));
1192 const ext::shared_ptr<PricingEngine> analyticHestonEngine(
1193 new AnalyticHestonEngine(analyticHestonModel, 164));
1194 const ext::shared_ptr<HestonModel> fdmHestonModel
1195 (new HestonModel(hestonProcess));
1196
1197
1198 const Real equityShortRateCorr = -0.5;
1199
1200 std::vector<Real> strikes = { 50, 75, 90, 100, 110, 125, 150, 200 };
1201 std::vector<Time> maturities = { 1/12., 3/12., 0.5, 1.0, 2.0, 3.0, 5.0, 7.5, 10};
1202
1203 std::vector<Volatility> vol = {
1204 0.482627,0.407617,0.366682,0.340110,0.314266,0.280241,0.252471,0.325552,
1205 0.464811,0.393336,0.354664,0.329758,0.305668,0.273563,0.244024,0.244886,
1206 0.441864,0.375618,0.340464,0.318249,0.297127,0.268839,0.237972,0.225553,
1207 0.407506,0.351125,0.322571,0.305173,0.289034,0.267361,0.239315,0.213761,
1208 0.366761,0.326166,0.306764,0.295279,0.284765,0.270592,0.250702,0.222928,
1209 0.345671,0.314748,0.300259,0.291744,0.283971,0.273475,0.258503,0.235683,
1210 0.324512,0.303631,0.293981,0.288338,0.283193,0.276248,0.266271,0.250506,
1211 0.311278,0.296340,0.289481,0.285482,0.281840,0.276924,0.269856,0.258609,
1212 0.303219,0.291534,0.286187,0.283073,0.280239,0.276414,0.270926,0.262173
1213 };
1214
1215 std::vector<ext::shared_ptr<CalibrationHelper> > options;
1216
1217 for (Size i=0; i < maturities.size(); ++i) {
1218 const Period maturity((int)std::lround(x: maturities[i]*12.0), Months);
1219 ext::shared_ptr<Exercise> exercise(
1220 new EuropeanExercise(today + maturity));
1221
1222 for (Size j=0; j < strikes.size(); ++j) {
1223 ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(
1224 strikes[j] * rTS->discount(t: maturities[i]) >=
1225 s0->value() * qTS->discount(t: maturities[i])
1226 ? Option::Call
1227 : Option::Put,
1228 strikes[j]));
1229 RelinkableHandle<Quote> v(ext::shared_ptr<Quote>(new SimpleQuote(vol[i*strikes.size()+j])));
1230
1231 ext::shared_ptr<BlackCalibrationHelper> helper(
1232 new HestonModelHelper(maturity, calendar, s0,
1233 strikes[j], v, rTS, qTS,
1234 BlackCalibrationHelper::PriceError));
1235 options.push_back(x: helper);
1236 const Real marketValue = helper->marketValue();
1237
1238 // Improve the quality of the starting point
1239 // for the full Heston-Hull-White calibration
1240 ext::shared_ptr<SimpleQuote> volQuote(new SimpleQuote);
1241 ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess =
1242 QuantLib::detail::ImpliedVolatilityHelper::clone(
1243 ext::make_shared<GeneralizedBlackScholesProcess>(
1244 args&: s0, args: qTS, args: rTS, args: Handle<BlackVolTermStructure>(
1245 flatVol(volatility: v->value(), dc))),
1246 volQuote);
1247
1248 VanillaOption dummyOption(payoff, exercise);
1249
1250 ext::shared_ptr<PricingEngine> bshwEngine(
1251 new AnalyticBSMHullWhiteEngine(equityShortRateCorr,
1252 bsProcess, hullWhiteModel));
1253
1254 Volatility vt = QuantLib::detail::ImpliedVolatilityHelper::calculate(
1255 instrument: dummyOption, engine: *bshwEngine, volQuote&: *volQuote,
1256 targetValue: marketValue, accuracy: 1e-8, maxEvaluations: 100, minVol: 0.0001, maxVol: 10);
1257
1258 v.linkTo(h: ext::shared_ptr<Quote>(new SimpleQuote(vt)));
1259
1260 helper->setPricingEngine(
1261 ext::shared_ptr<PricingEngine>(analyticHestonEngine));
1262 }
1263 }
1264
1265 HestonHullWhiteCorrelationConstraint corrConstraint(equityShortRateCorr);
1266 LevenbergMarquardt om(1e-6, 1e-8, 1e-8);
1267 analyticHestonModel->calibrate(options, method&: om,
1268 endCriteria: EndCriteria(400, 40, 1.0e-8, 1.0e-4, 1.0e-8),
1269 constraint: corrConstraint);
1270
1271 options.clear();
1272 fdmHestonModel->setParams(analyticHestonModel->params());
1273
1274 for (Size i=0; i < maturities.size(); ++i) {
1275 const Size tGrid = static_cast<Size>(std::max(a: 5.0, b: maturities[i]*5.0));
1276 ext::shared_ptr<FdHestonHullWhiteVanillaEngine> engine(
1277 new FdHestonHullWhiteVanillaEngine(fdmHestonModel, hwProcess,
1278 equityShortRateCorr,
1279 tGrid, 45, 11, 5, 0, true));
1280
1281 engine->enableMultipleStrikesCaching(strikes);
1282
1283 const Period maturity((int)std::lround(x: maturities[i]*12.0), Months);
1284
1285 for (Size j=0; j < strikes.size(); ++j) {
1286 // multiple strikes engine works best if the first option
1287 // per maturity has the average strike (because the first option
1288 // is priced first during the calibration and the first pricing
1289 // is used to calculate the prices for all strikes
1290 const Size js = (j + (strikes.size()-1)/2) % strikes.size();
1291
1292 ext::shared_ptr<StrikedTypePayoff> payoff(
1293 new PlainVanillaPayoff(Option::Call, strikes[js]));
1294 Handle<Quote> v(ext::shared_ptr<Quote>(new SimpleQuote(vol[i*strikes.size()+js])));
1295 ext::shared_ptr<BlackCalibrationHelper> helper(
1296 new HestonModelHelper(maturity, calendar, s0,
1297 strikes[js], v, rTS, qTS,
1298 BlackCalibrationHelper::PriceError));
1299 options.push_back(x: helper);
1300
1301 helper->setPricingEngine(engine);
1302 }
1303 }
1304
1305 LevenbergMarquardt vm(1e-6, 1e-2, 1e-2);
1306 fdmHestonModel->calibrate(options, method&: vm,
1307 endCriteria: EndCriteria(400, 40, 1.0e-8, 1.0e-4, 1.0e-8),
1308 constraint: corrConstraint);
1309
1310 const Real relTol = 0.01;
1311 const Real expected_v0 = 0.12;
1312 const Real expected_kappa = 2.0;
1313 const Real expected_theta = 0.09;
1314 const Real expected_sigma = 0.5;
1315 const Real expected_rho = -0.75;
1316
1317 if (std::fabs(x: fdmHestonModel->v0() - expected_v0)/expected_v0 > relTol) {
1318 BOOST_ERROR("Failed to reproduce Heston-Hull-White model"
1319 << "\n v0 calculated: " << fdmHestonModel->v0()
1320 << "\n v0 expected : " << expected_v0
1321 << "\n relatove tol : " << relTol);
1322 }
1323 if (std::fabs(x: fdmHestonModel->theta() - expected_theta)/expected_theta > relTol) {
1324 BOOST_ERROR("Failed to reproduce Heston-Hull-White model"
1325 << "\n theta calculated: " << fdmHestonModel->theta()
1326 << "\n theta expected : " << expected_theta
1327 << "\n relatove tol : " << relTol);
1328 }
1329 if (std::fabs(x: fdmHestonModel->kappa() - expected_kappa)/expected_kappa > relTol) {
1330 BOOST_ERROR("Failed to reproduce Heston-Hull-White model"
1331 << "\n kappa calculated: " << fdmHestonModel->kappa()
1332 << "\n kappa expected : " << expected_kappa
1333 << "\n relatove tol : " << relTol);
1334 }
1335 if (std::fabs(x: fdmHestonModel->sigma() - expected_sigma)/expected_sigma > relTol) {
1336 BOOST_ERROR("Failed to reproduce Heston-Hull-White model"
1337 << "\n sigma calculated: " << fdmHestonModel->sigma()
1338 << "\n sigma expected : " << expected_sigma
1339 << "\n relatove tol : " << relTol);
1340 }
1341 if (std::fabs(x: fdmHestonModel->rho() - expected_rho)/expected_rho > relTol) {
1342 BOOST_ERROR("Failed to reproduce Heston-Hull-White model"
1343 << "\n rho calculated: " << fdmHestonModel->rho()
1344 << "\n rho expected : " << expected_rho
1345 << "\n relatove tol : " << relTol);
1346 }
1347}
1348
1349void HybridHestonHullWhiteProcessTest::testH1HWPricingEngine() {
1350 BOOST_TEST_MESSAGE("Testing the H1-HW approximation engine...");
1351
1352 /*
1353 * Example taken from Lech Aleksander Grzelak,
1354 * Equity and Foreign Exchange Hybrid Models for Pricing Long-Maturity
1355 * Financial Derivatives,
1356 * http://repository.tudelft.nl/assets/uuid:a8e1a007-bd89-481a-aee3-0e22f15ade6b/PhDThesis_main.pdf
1357 */
1358
1359 const Date today = Date(15, July, 2012);
1360 Settings::instance().evaluationDate() = today;
1361 const Date exerciseDate = Date(13, July, 2022);
1362 const DayCounter dc = Actual365Fixed();
1363
1364 const ext::shared_ptr<Exercise> exercise(
1365 new EuropeanExercise(exerciseDate));
1366
1367 const Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
1368
1369 const Real r = 0.02;
1370 const Real q = 0.00;
1371 const Real v0 = 0.05;
1372 const Real theta = 0.05;
1373 const Real kappa_v = 0.3;
1374 const Real sigma_v[] = { 0.3, 0.6 };
1375 const Real rho_sv =-0.30;
1376 const Real rho_sr = 0.6;
1377 const Real kappa_r = 0.01;
1378 const Real sigma_r = 0.01;
1379
1380 const Handle<YieldTermStructure> rTS(flatRate(today, forward: r, dc));
1381 const Handle<YieldTermStructure> qTS(flatRate(today, forward: q, dc));
1382
1383 const Handle<BlackVolTermStructure> flatVolTS(flatVol(today, volatility: 0.20, dc));
1384 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
1385 new GeneralizedBlackScholesProcess(s0, qTS, rTS, flatVolTS));
1386
1387 const ext::shared_ptr<HullWhiteProcess> hwProcess(
1388 new HullWhiteProcess(rTS, kappa_r, sigma_r));
1389 const ext::shared_ptr<HullWhite> hullWhiteModel(
1390 new HullWhite(Handle<YieldTermStructure>(rTS), kappa_r, sigma_r));
1391
1392 const Real tol = 0.0001;
1393 const Real strikes[] = {40, 80, 100, 120, 180 };
1394 const Real expected[LENGTH(sigma_v)][LENGTH(strikes)]
1395 = { { 0.267503, 0.235742, 0.228223, 0.223461, 0.217855 },
1396 { 0.263626, 0.211625, 0.199907, 0.193502, 0.190025 } };
1397
1398 for (Size j=0; j < LENGTH(sigma_v); ++j) {
1399 const ext::shared_ptr<HestonProcess> hestonProcess(
1400 new HestonProcess(rTS, qTS, s0, v0, kappa_v, theta,
1401 sigma_v[j], rho_sv));
1402 const ext::shared_ptr<HestonModel> hestonModel(
1403 new HestonModel(hestonProcess));
1404
1405 for (Size i=0; i < LENGTH(strikes); ++i) {
1406 const ext::shared_ptr<StrikedTypePayoff> payoff(
1407 new PlainVanillaPayoff(Option::Call, strikes[i]));
1408
1409 VanillaOption option(payoff, exercise);
1410
1411 const ext::shared_ptr<PricingEngine> analyticH1HWEngine(
1412 new AnalyticH1HWEngine(hestonModel, hullWhiteModel,
1413 rho_sr, 144));
1414 option.setPricingEngine(analyticH1HWEngine);
1415 const Real impliedH1HW
1416 = option.impliedVolatility(price: option.NPV(), process: bsProcess);
1417
1418 if (std::fabs(x: expected[j][i] - impliedH1HW) > tol) {
1419 BOOST_ERROR("Failed to reproduce H1HW implied volatility"
1420 << "\n expected : " << expected[j][i]
1421 << "\n calculated : " << impliedH1HW
1422 << "\n tol : " << tol
1423 << "\n strike : " << strikes[i]
1424 << "\n sigma : " << sigma_v[j]);
1425 }
1426 }
1427 }
1428}
1429
1430test_suite* HybridHestonHullWhiteProcessTest::suite(SpeedLevel speed) {
1431 auto* suite = BOOST_TEST_SUITE("Hybrid Heston-HullWhite tests");
1432
1433 suite->add(QUANTLIB_TEST_CASE(
1434 &HybridHestonHullWhiteProcessTest::testBsmHullWhiteEngine));
1435 suite->add(QUANTLIB_TEST_CASE(
1436 &HybridHestonHullWhiteProcessTest::testCompareBsmHWandHestonHW));
1437 suite->add(QUANTLIB_TEST_CASE(
1438 &HybridHestonHullWhiteProcessTest::testZeroBondPricing));
1439 suite->add(QUANTLIB_TEST_CASE(
1440 &HybridHestonHullWhiteProcessTest::testMcVanillaPricing));
1441 suite->add(QUANTLIB_TEST_CASE(
1442 &HybridHestonHullWhiteProcessTest::testMcPureHestonPricing));
1443 suite->add(QUANTLIB_TEST_CASE(
1444 &HybridHestonHullWhiteProcessTest::testAnalyticHestonHullWhitePricing));
1445 suite->add(QUANTLIB_TEST_CASE(
1446 &HybridHestonHullWhiteProcessTest::testCallableEquityPricing));
1447 suite->add(QUANTLIB_TEST_CASE(
1448 &HybridHestonHullWhiteProcessTest::testDiscretizationError));
1449 suite->add(QUANTLIB_TEST_CASE(
1450 &HybridHestonHullWhiteProcessTest::testBsmHullWhitePricing));
1451 suite->add(QUANTLIB_TEST_CASE(
1452 &HybridHestonHullWhiteProcessTest::testH1HWPricingEngine));
1453
1454 if (speed <= Fast) {
1455 suite->add(QUANTLIB_TEST_CASE(
1456 &HybridHestonHullWhiteProcessTest::testSpatialDiscretizatinError));
1457 suite->add(QUANTLIB_TEST_CASE(
1458 &HybridHestonHullWhiteProcessTest::testFdmHestonHullWhiteEngine));
1459 }
1460
1461 if (speed == Slow) {
1462 suite->add(QUANTLIB_TEST_CASE(
1463 &HybridHestonHullWhiteProcessTest::testHestonHullWhiteCalibration));
1464 }
1465
1466 return suite;
1467}
1468

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