1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2016 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
21#include "normalclvmodel.hpp"
22#include "utilities.hpp"
23#include <ql/pricingengines/barrier/analyticdoublebarrierbinaryengine.hpp>
24#include <ql/instruments/doublebarrieroption.hpp>
25#include <ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp>
26#include <ql/experimental/models/normalclvmodel.hpp>
27#include <ql/experimental/volatility/sabrvoltermstructure.hpp>
28#include <ql/functional.hpp>
29#include <ql/instruments/forwardvanillaoption.hpp>
30#include <ql/instruments/impliedvolatility.hpp>
31#include <ql/math/integrals/gausslobattointegral.hpp>
32#include <ql/math/randomnumbers/rngtraits.hpp>
33#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
34#include <ql/math/statistics/statistics.hpp>
35#include <ql/methods/finitedifferences/utilities/bsmrndcalculator.hpp>
36#include <ql/methods/finitedifferences/utilities/hestonrndcalculator.hpp>
37#include <ql/methods/montecarlo/pathgenerator.hpp>
38#include <ql/pricingengines/blackcalculator.hpp>
39#include <ql/pricingengines/forward/forwardengine.hpp>
40#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
41#include <ql/processes/blackscholesprocess.hpp>
42#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
43#include <ql/quotes/simplequote.hpp>
44#include <ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp>
45#include <ql/time/calendars/nullcalendar.hpp>
46#include <ql/time/daycounters/actual360.hpp>
47#include <ql/time/daycounters/actual365fixed.hpp>
48#include <ql/time/daycounters/actualactual.hpp>
49#include <utility>
50
51using namespace QuantLib;
52using namespace boost::unit_test_framework;
53
54void NormalCLVModelTest::testBSCumlativeDistributionFunction() {
55 BOOST_TEST_MESSAGE("Testing Black-Scholes cumulative distribution function"
56 " with constant volatility...");
57
58 const DayCounter dc = Actual365Fixed();
59 const Date today = Date(22, June, 2016);
60 const Date maturity = today + Period(6, Months);
61
62 const Real s0 = 100;
63 const Real rRate = 0.1;
64 const Real qRate = 0.05;
65 const Volatility vol = 0.25;
66
67 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
68 const Handle<YieldTermStructure> qTS(flatRate(today, forward: qRate, dc));
69 const Handle<YieldTermStructure> rTS(flatRate(today, forward: rRate, dc));
70 const Handle<BlackVolTermStructure> volTS(flatVol(today, volatility: vol, dc));
71
72 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
73 ext::make_shared<GeneralizedBlackScholesProcess>(
74 args: spot, args: qTS, args: rTS, args: volTS));
75 const ext::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess;
76
77 const NormalCLVModel m(
78 bsProcess, ouProcess, std::vector<Date>(), 5);
79 const BSMRNDCalculator rndCalculator(bsProcess);
80
81
82 constexpr double tol = 1e5 * QL_EPSILON;
83 const Time t = dc.yearFraction(d1: today, d2: maturity);
84 for (Real x=10; x < 400; x+=10) {
85 const Real calculated = m.cdf(d: maturity, x);
86 const Real expected = rndCalculator.cdf(x: std::log(x: x), t);
87
88 if (std::fabs(x: calculated - expected) > tol) {
89 BOOST_FAIL("Failed to reproduce CDF for "
90 << "\n strike: " << x
91 << "\n calculated: " << calculated
92 << "\n expected: " << expected);
93 }
94 }
95}
96
97void NormalCLVModelTest::testHestonCumlativeDistributionFunction() {
98 BOOST_TEST_MESSAGE("Testing Heston cumulative distribution function...");
99
100 const DayCounter dc = Actual365Fixed();
101 const Date today = Date(22, June, 2016);
102 const Date maturity = today + Period(1, Years);
103
104 const Real s0 = 100;
105 const Real v0 = 0.01;
106 const Real rRate = 0.1;
107 const Real qRate = 0.05;
108 const Real kappa = 2.0;
109 const Real theta = 0.09;
110 const Real sigma = 0.4;
111 const Real rho = -0.75;
112
113 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
114 const Handle<YieldTermStructure> qTS(flatRate(today, forward: qRate, dc));
115 const Handle<YieldTermStructure> rTS(flatRate(today, forward: rRate, dc));
116
117 const ext::shared_ptr<HestonProcess> process(
118 ext::make_shared<HestonProcess>(
119 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
120
121 const Handle<BlackVolTermStructure> hestonVolTS(
122 ext::make_shared<HestonBlackVolSurface>(
123 args: Handle<HestonModel>(ext::make_shared<HestonModel>(args: process))));
124
125 const NormalCLVModel m(
126 ext::make_shared<GeneralizedBlackScholesProcess>(
127 args: spot, args: qTS, args: rTS, args: hestonVolTS),
128 ext::shared_ptr<OrnsteinUhlenbeckProcess>(),
129 std::vector<Date>(), 5);
130
131 const HestonRNDCalculator rndCalculator(process);
132
133 const Real tol = 1e-6;
134 const Time t = dc.yearFraction(d1: today, d2: maturity);
135 for (Real x=10; x < 400; x+=25) {
136 const Real calculated = m.cdf(d: maturity, x);
137 const Real expected = rndCalculator.cdf(x: std::log(x: x), t);
138
139 if (std::fabs(x: calculated - expected) > tol) {
140 BOOST_FAIL("Failed to reproduce CDF for "
141 << "\n strike: " << x
142 << "\n calculated: " << calculated
143 << "\n expected: " << expected);
144 }
145 }
146}
147
148
149
150void NormalCLVModelTest::testIllustrative1DExample() {
151 BOOST_TEST_MESSAGE(
152 "Testing illustrative 1D example of normal CLV model...");
153
154 // example taken from:
155 // A. Grzelak, 2015, The CLV Framework -
156 // A Fresh Look at Efficient Pricing with Smile
157 // http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2747541
158
159 const DayCounter dc = Actual360();
160 const Date today = Date(22, June, 2016);
161
162 //SABR
163 const Real beta = 0.5;
164 const Real alpha= 0.2;
165 const Real rho = -0.9;
166 const Real gamma= 0.2;
167
168 // Ornstein-Uhlenbeck
169 const Real speed = 1.3;
170 const Real level = 0.1;
171 const Real vol = 0.25;
172 const Real x0 = 1.0;
173
174 const Real s0 = 1.0;
175 const Real rRate = 0.03;
176 const Real qRate = 0.0;
177
178 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
179 const Handle<YieldTermStructure> qTS(flatRate(today, forward: qRate, dc));
180 const Handle<YieldTermStructure> rTS(flatRate(today, forward: rRate, dc));
181
182 const Handle<BlackVolTermStructure> sabrVol(
183 ext::make_shared<SABRVolTermStructure>(
184 args: alpha, args: beta, args: gamma, args: rho, args: s0, args: rRate, args: today, args: dc));
185
186 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
187 ext::make_shared<GeneralizedBlackScholesProcess>(
188 args: spot, args: qTS, args: rTS, args: sabrVol));
189
190 const ext::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess(
191 ext::make_shared<OrnsteinUhlenbeckProcess>(
192 args: speed, args: vol, args: x0, args: level));
193
194 std::vector<Date> maturityDates = {
195 today + Period(18, Days),
196 today + Period(90, Days),
197 today + Period(180, Days),
198 today + Period(360, Days),
199 today + Period(720, Days)
200 };
201
202 const NormalCLVModel m(bsProcess, ouProcess, maturityDates, 4);
203 const ext::function<Real(Real, Real)> g = m.g();
204
205 // test collocation points in x_ij
206 std::vector<Date> maturities = { maturityDates[0], maturityDates[2], maturityDates[4] };
207
208 std::vector<std::vector<Real> > x = {
209 { 1.070, 0.984, 0.903, 0.817 },
210 { 0.879, 0.668, 0.472, 0.261 },
211 { 0.528, 0.282, 0.052,-0.194 }
212 };
213
214 std::vector<std::vector<Real> > s = {
215 {1.104, 1.035, 0.969, 0.895},
216 {1.328, 1.122, 0.911, 0.668},
217 {1.657, 1.283, 0.854, 0.339}
218 };
219
220 std::vector<Real> c = { 2.3344, 0.7420, -0.7420, -2.3344 };
221
222 const Real tol = 0.001;
223 for (Size i=0; i < maturities.size(); ++i) {
224 const Time t = dc.yearFraction(d1: today, d2: maturities[i]);
225
226 for (Size j=0; j < x.front().size(); ++j) {
227 const Real calculatedX = m.collocationPointsX(d: maturities[i])[j];
228 const Real expectedX = x[i][j];
229
230 if (std::fabs(x: calculatedX - expectedX) > tol) {
231 BOOST_FAIL("Failed to reproduce collocation x points for "
232 << "\n time: " << maturities[i]
233 << "\n j " << j
234 << "\n calculated: " << calculatedX
235 << "\n expected: " << expectedX);
236 }
237
238 const Real calculatedS = m.collocationPointsY(d: maturities[i])[j];
239 const Real expectedS = s[i][j];
240 if (std::fabs(x: calculatedS - expectedS) > tol) {
241 BOOST_FAIL("Failed to reproduce collocation s points for "
242 << "\n time: " << maturities[i]
243 << "\n j " << j
244 << "\n calculated: " << calculatedS
245 << "\n expected: " << expectedS);
246 }
247
248 const Real expectation
249 = ouProcess->expectation(0.0, x0: ouProcess->x0(), dt: t);
250 const Real stdDeviation
251 = ouProcess->stdDeviation(t: 0.0, x0: ouProcess->x0(), dt: t);
252
253 const Real calculatedG = g(t, expectation + stdDeviation*c[j]);
254 if (std::fabs(x: calculatedG - expectedS) > tol) {
255 BOOST_FAIL("Failed to reproduce g values "
256 "at collocation points for "
257 << "\n time: " << maturities[i]
258 << "\n j " << j
259 << "\n calculated: " << calculatedG
260 << "\n expected: " << expectedS);
261 }
262 }
263 }
264}
265
266namespace normal_clv_model_test {
267 class CLVModelPayoff : public PlainVanillaPayoff {
268 public:
269 CLVModelPayoff(Option::Type type, Real strike, ext::function<Real(Real)> g)
270 : PlainVanillaPayoff(type, strike), g_(std::move(g)) {}
271
272 Real operator()(Real x) const override { return PlainVanillaPayoff::operator()(price: g_(x)); }
273
274 private:
275 const ext::function<Real(Real)> g_;
276 };
277}
278
279void NormalCLVModelTest::testMonteCarloBSOptionPricing() {
280 BOOST_TEST_MESSAGE("Testing Monte Carlo BS option pricing...");
281
282 using namespace normal_clv_model_test;
283
284 const DayCounter dc = Actual365Fixed();
285 const Date today = Date(22, June, 2016);
286 const Date maturity = today + Period(1, Years);
287 const Time t = dc.yearFraction(d1: today, d2: maturity);
288
289 const Real strike = 110;
290 const ext::shared_ptr<StrikedTypePayoff> payoff =
291 ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike);
292 const ext::shared_ptr<Exercise> exercise =
293 ext::make_shared<EuropeanExercise>(args: maturity);
294
295 // Ornstein-Uhlenbeck
296 const Real speed = 2.3;
297 const Real level = 100;
298 const Real sigma = 0.35;
299 const Real x0 = 100.0;
300
301 const Real s0 = x0;
302 const Volatility vol = 0.25;
303 const Real rRate = 0.10;
304 const Real qRate = 0.04;
305
306 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
307 const Handle<YieldTermStructure> qTS(flatRate(today, forward: qRate, dc));
308 const Handle<YieldTermStructure> rTS(flatRate(today, forward: rRate, dc));
309 const Handle<BlackVolTermStructure> vTS(flatVol(today, volatility: vol, dc));
310
311 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
312 ext::make_shared<GeneralizedBlackScholesProcess>(
313 args: spot, args: qTS, args: rTS, args: vTS));
314
315 const ext::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess(
316 ext::make_shared<OrnsteinUhlenbeckProcess>(
317 args: speed, args: sigma, args: x0, args: level));
318
319 std::vector<Date> maturities = { today + Period(6, Months), maturity };
320
321 const NormalCLVModel m(bsProcess, ouProcess, maturities, 8);
322 const ext::function<Real(Real, Real)> g = m.g();
323
324 const Size nSims = 32767;
325 const LowDiscrepancy::rsg_type ld
326 = LowDiscrepancy::make_sequence_generator(dimension: 1, seed: 23455);
327
328 Statistics stat;
329 for (Size i=0; i < nSims; ++i) {
330 const Real dw = ld.nextSequence().value.front();
331
332 const Real o_t = ouProcess->evolve(t0: 0, x0, dt: t, dw);
333 const Real s = g(t, o_t);
334
335 stat.add(value: (*payoff)(s));
336
337 }
338
339 Real calculated = stat.mean() * rTS->discount(d: maturity);
340
341 VanillaOption option(payoff, exercise);
342 option.setPricingEngine(
343 ext::make_shared<AnalyticEuropeanEngine>(args: bsProcess));
344 const Real expected = option.NPV();
345
346 const Real tol = 0.01;
347 if (std::fabs(x: calculated - expected) > tol) {
348 BOOST_FAIL("Failed to reproduce Monte-Carlo vanilla option price "
349 << "\n time: " << maturity
350 << "\n strike: " << strike
351 << "\n calculated: " << calculated
352 << "\n expected: " << expected);
353 }
354
355 VanillaOption fdmOption(
356 ext::make_shared<CLVModelPayoff>(args: payoff->optionType(), args: payoff->strike(),
357 args: [&](Real _x) { return g(t, _x); }),
358 exercise);
359
360 fdmOption.setPricingEngine(
361 ext::make_shared<FdOrnsteinUhlenbeckVanillaEngine>(
362 args: ouProcess, args: rTS.currentLink(), args: 50, args: 800));
363
364 calculated = fdmOption.NPV();
365 if (std::fabs(x: calculated - expected) > tol) {
366 BOOST_FAIL("Failed to reproduce FDM vanilla option price "
367 << "\n time: " << maturity
368 << "\n strike: " << strike
369 << "\n calculated: " << calculated
370 << "\n expected: " << expected);
371 }
372}
373
374void NormalCLVModelTest::testMoustacheGraph() {
375 BOOST_TEST_MESSAGE(
376 "Testing double no-touch pricing with normal CLV model...");
377
378 /*
379 The comparison of Black-Scholes and normal CLV prices is derived
380 from figure 8.8 in Iain J. Clark's book,
381 Foreign Exchange Option Pricing: A Practitioner’s Guide
382 */
383
384 const DayCounter dc = ActualActual(ActualActual::ISDA);
385 const Date todaysDate(5, Aug, 2016);
386 const Date maturityDate = todaysDate + Period(1, Years);
387 const Time maturityTime = dc.yearFraction(d1: todaysDate, d2: maturityDate);
388
389 Settings::instance().evaluationDate() = todaysDate;
390
391 const Real s0 = 100;
392 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
393 const Rate r = 0.02;
394 const Rate q = 0.01;
395
396 // parameter of the "calibrated" Heston model
397 const Real kappa = 1.0;
398 const Real theta = 0.06;
399 const Real rho = -0.8;
400 const Real sigma = 0.8;
401 const Real v0 = 0.09;
402
403 const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
404 const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
405
406 const ext::shared_ptr<HestonModel> hestonModel(
407 ext::make_shared<HestonModel>(
408 args: ext::make_shared<HestonProcess>(
409 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
410
411 const Handle<BlackVolTermStructure> vTS(
412 ext::make_shared<HestonBlackVolSurface>(
413 args: Handle<HestonModel>(hestonModel)));
414
415 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess =
416 ext::make_shared<GeneralizedBlackScholesProcess>(
417 args: spot, args: qTS, args: rTS, args: vTS);
418
419 // Ornstein-Uhlenbeck
420 const Real speed = -0.80;
421 const Real level = 100;
422 const Real sigmaOU = 0.15;
423 const Real x0 = 100;
424
425 const ext::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess(
426 ext::make_shared<OrnsteinUhlenbeckProcess>(
427 args: speed, args: sigmaOU, args: x0, args: level));
428
429 const ext::shared_ptr<Exercise> europeanExercise(
430 ext::make_shared<EuropeanExercise>(args: maturityDate));
431
432 VanillaOption vanillaOption(
433 ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: s0),
434 europeanExercise);
435
436 vanillaOption.setPricingEngine(
437 ext::make_shared<AnalyticHestonEngine>(args: hestonModel));
438
439 const Volatility atmVol = vanillaOption.impliedVolatility(
440 price: vanillaOption.NPV(),
441 process: ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS,
442 args: Handle<BlackVolTermStructure>(flatVol(volatility: std::sqrt(x: theta), dc))));
443
444 const ext::shared_ptr<PricingEngine> analyticEngine(
445 ext::make_shared<AnalyticDoubleBarrierBinaryEngine>(
446 args: ext::make_shared<GeneralizedBlackScholesProcess>(
447 args: spot, args: qTS, args: rTS,
448 args: Handle<BlackVolTermStructure>(flatVol(volatility: atmVol, dc)))));
449
450
451 std::vector<Date> maturities(1, todaysDate + Period(2, Weeks));
452 while (maturities.back() < maturityDate)
453 maturities.push_back(x: maturities.back() + Period(2, Weeks));
454
455 const NormalCLVModel m(bsProcess, ouProcess, maturities, 8);
456 const ext::function<Real(Real, Real)> g = m.g();
457
458 const Size n = 18;
459 Array barrier_lo(n), barrier_hi(n), bsNPV(n);
460
461 const ext::shared_ptr<CashOrNothingPayoff> payoff =
462 ext::make_shared<CashOrNothingPayoff>(args: Option::Call, args: 0.0, args: 1.0);
463
464 for (Size i=0; i < n; ++i) {
465 const Real dist = 10.0+5.0*i;
466
467 barrier_lo[i] = std::max(a: s0 - dist, b: 1e-2);
468 barrier_hi[i] = s0 + dist;
469 DoubleBarrierOption doubleBarrier(
470 DoubleBarrier::KnockOut, barrier_lo[i], barrier_hi[i], 0.0,
471 payoff,
472 europeanExercise);
473
474 doubleBarrier.setPricingEngine(analyticEngine);
475 bsNPV[i] = doubleBarrier.NPV();
476 }
477
478 typedef SobolBrownianBridgeRsg rsg_type;
479 typedef PathGenerator<rsg_type>::sample_type sample_type;
480
481 const Size factors = 1;
482 const Size tSteps = 200;
483 const TimeGrid grid(maturityTime, tSteps);
484
485 const ext::shared_ptr<PathGenerator<rsg_type> > pathGenerator =
486 ext::make_shared<PathGenerator<rsg_type> >(
487 args: ouProcess, args: grid, args: rsg_type(factors, tSteps), args: false);
488
489 const Size nSims = 100000;
490 std::vector<GeneralStatistics> stats(n);
491 const DiscountFactor df = rTS->discount(d: maturityDate);
492
493 for (Size i=0; i < nSims; ++i) {
494 std::vector<bool> touch(n, false);
495
496 const sample_type& path = pathGenerator->next();
497
498 Real s;
499 for (Size j=1; j <= tSteps; ++j) {
500 const Time t = grid.at(i: j);
501 s = g(t, path.value.at(i: j));
502
503 for (Size u=0; u < n; ++u) {
504 if (s <= barrier_lo[u] || s >= barrier_hi[u]) {
505 touch[u] = true;
506 }
507 }
508 }
509 for (Size u=0; u < n; ++u) {
510 if (touch[u]) {
511 stats[u].add(value: 0.0);
512 }
513 else {
514 stats[u].add(value: df*(*payoff)(s));
515 }
516 }
517 }
518
519 const Real expected[] = {
520 0.00931214, 0.0901481, 0.138982, 0.112059, 0.0595901,
521 0.0167549, -0.00906787, -0.0206768, -0.0225628, -0.0203593,
522 -0.016036, -0.0116629, -0.00728792, -0.00328821,
523 -0.000158562, 0.00502041, 0.00347706, 0.00238216, };
524
525 const Real tol = 1e-5;
526 for (Size u=0; u < n; ++u) {
527 const Real calculated = stats[u].mean() - bsNPV[u];
528
529 if (std::fabs(x: calculated - expected[u]) > tol) {
530 BOOST_FAIL("Failed to reproduce Double no Touch prices"
531 << "\n time: " << maturityDate
532 << "\n barrier lower: " << barrier_lo[u]
533 << "\n barrier high: " << barrier_hi[u]
534 << "\n calculated: " << calculated
535 << "\n expected: " << expected[u]);
536 }
537 }
538}
539
540test_suite* NormalCLVModelTest::experimental(SpeedLevel speed) {
541 auto* suite = BOOST_TEST_SUITE("NormalCLVModel tests");
542
543 suite->add(QUANTLIB_TEST_CASE(
544 &NormalCLVModelTest::testBSCumlativeDistributionFunction));
545 suite->add(QUANTLIB_TEST_CASE(
546 &NormalCLVModelTest::testHestonCumlativeDistributionFunction));
547 suite->add(QUANTLIB_TEST_CASE(
548 &NormalCLVModelTest::testIllustrative1DExample));
549 suite->add(QUANTLIB_TEST_CASE(
550 &NormalCLVModelTest::testMonteCarloBSOptionPricing));
551
552 if (speed == Slow) {
553 suite->add(QUANTLIB_TEST_CASE(
554 &NormalCLVModelTest::testMoustacheGraph));
555 }
556 return suite;
557}
558

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