1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2015 Johannes Göttker-Schnetmann
5 Copyright (C) 2015, 2016 Klaus Spanderen
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 "riskneutraldensitycalculator.hpp"
22#include "utilities.hpp"
23#include <ql/instruments/vanillaoption.hpp>
24#include <ql/math/distributions/normaldistribution.hpp>
25#include <ql/math/integrals/gausslobattointegral.hpp>
26#include <ql/methods/finitedifferences/utilities/bsmrndcalculator.hpp>
27#include <ql/methods/finitedifferences/utilities/cevrndcalculator.hpp>
28#include <ql/methods/finitedifferences/utilities/gbsmrndcalculator.hpp>
29#include <ql/methods/finitedifferences/utilities/hestonrndcalculator.hpp>
30#include <ql/methods/finitedifferences/utilities/localvolrndcalculator.hpp>
31#include <ql/methods/finitedifferences/utilities/squarerootprocessrndcalculator.hpp>
32#include <ql/models/equity/hestonmodel.hpp>
33#include <ql/pricingengines/blackcalculator.hpp>
34#include <ql/pricingengines/vanilla/fdblackscholesvanillaengine.hpp>
35#include <ql/processes/blackscholesprocess.hpp>
36#include <ql/processes/hestonprocess.hpp>
37#include <ql/quotes/simplequote.hpp>
38#include <ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp>
39#include <ql/termstructures/volatility/equityfx/localconstantvol.hpp>
40#include <ql/termstructures/volatility/equityfx/noexceptlocalvolsurface.hpp>
41#include <ql/time/calendars/nullcalendar.hpp>
42#include <ql/timegrid.hpp>
43#include <ql/types.hpp>
44#include <utility>
45
46using namespace QuantLib;
47using namespace boost::unit_test_framework;
48
49void RiskNeutralDensityCalculatorTest::testDensityAgainstOptionPrices() {
50 BOOST_TEST_MESSAGE("Testing density against option prices...");
51
52 const DayCounter dayCounter = Actual365Fixed();
53 const Date todaysDate = Settings::instance().evaluationDate();
54
55 const Real s0 = 100;
56 const Handle<Quote> spot(
57 ext::make_shared<SimpleQuote>(args: s0));
58
59 const Rate r = 0.075;
60 const Rate q = 0.04;
61 const Volatility v = 0.27;
62
63 const Handle<YieldTermStructure> rTS(flatRate(today: todaysDate, forward: r, dc: dayCounter));
64
65 const Handle<YieldTermStructure> qTS(flatRate(today: todaysDate, forward: q, dc: dayCounter));
66
67 const ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
68 new BlackScholesMertonProcess(
69 spot, qTS, rTS,
70 Handle<BlackVolTermStructure>(flatVol(volatility: v, dc: dayCounter))));
71
72 const BSMRNDCalculator bsm(bsmProcess);
73 const Time times[] = { 0.5, 1.0, 2.0 };
74 const Real strikes[] = { 75.0, 100.0, 150.0 };
75
76 for (Real t : times) {
77 const Volatility stdDev = v * std::sqrt(x: t);
78 const DiscountFactor df = rTS->discount(t);
79 const Real fwd = s0*qTS->discount(t)/df;
80
81 for (Real strike : strikes) {
82 const Real xs = std::log(x: strike);
83 const BlackCalculator blackCalc(
84 Option::Put, strike, fwd, stdDev, df);
85
86 const Real tol = 10*std::sqrt(QL_EPSILON);
87 const Real calculatedCDF = bsm.cdf(x: xs, t);
88 const Real expectedCDF
89 = blackCalc.strikeSensitivity()/df;
90
91 if (std::fabs(x: calculatedCDF - expectedCDF) > tol) {
92 BOOST_FAIL("failed to reproduce Black-Scholes-Merton cdf"
93 << "\n calculated: " << calculatedCDF
94 << "\n expected: " << expectedCDF
95 << "\n diff: " << calculatedCDF - expectedCDF
96 << "\n tol: " << tol);
97 }
98
99 const Real deltaStrike = strike*std::sqrt(QL_EPSILON);
100
101 const Real calculatedPDF = bsm.pdf(x: xs, t);
102 const Real expectedPDF = strike/df*
103 ( BlackCalculator(Option::Put, strike+deltaStrike,
104 fwd, stdDev, df).strikeSensitivity()
105 - BlackCalculator(Option::Put, strike - deltaStrike,
106 fwd, stdDev, df).strikeSensitivity())/(2*deltaStrike);
107
108 if (std::fabs(x: calculatedPDF - expectedPDF) > tol) {
109 BOOST_FAIL("failed to reproduce Black-Scholes-Merton pdf"
110 << "\n calculated: " << calculatedPDF
111 << "\n expected: " << expectedPDF
112 << "\n diff: " << calculatedPDF - expectedPDF
113 << "\n tol: " << tol);
114 }
115 }
116 }
117}
118
119void RiskNeutralDensityCalculatorTest::testBSMagainstHestonRND() {
120 BOOST_TEST_MESSAGE("Testing Black-Scholes-Merton and Heston densities...");
121
122 const DayCounter dayCounter = Actual365Fixed();
123 const Date todaysDate = Settings::instance().evaluationDate();
124
125 const Real s0 = 10;
126 const Handle<Quote> spot(
127 ext::make_shared<SimpleQuote>(args: s0));
128
129 const Rate r = 0.155;
130 const Rate q = 0.0721;
131 const Volatility v = 0.27;
132
133 const Real kappa = 1.0;
134 const Real theta = v*v;
135 const Real rho = -0.75;
136 const Real v0 = v*v;
137 const Real sigma = 0.0001;
138
139 const Handle<YieldTermStructure> rTS(flatRate(today: todaysDate, forward: r, dc: dayCounter));
140
141 const Handle<YieldTermStructure> qTS(flatRate(today: todaysDate, forward: q, dc: dayCounter));
142
143 const ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
144 new BlackScholesMertonProcess(
145 spot, qTS, rTS,
146 Handle<BlackVolTermStructure>(flatVol(volatility: v, dc: dayCounter))));
147
148 const BSMRNDCalculator bsm(bsmProcess);
149 const HestonRNDCalculator heston(
150 ext::make_shared<HestonProcess>(
151 args: rTS, args: qTS, args: spot,
152 args: v0, args: kappa, args: theta, args: sigma, args: rho), 1e-8);
153
154 const Time times[] = { 0.5, 1.0, 2.0 };
155 const Real strikes[] = { 7.5, 10, 15 };
156 const Real probs[] = { 1e-6, 0.01, 0.5, 0.99, 1.0-1e-6 };
157
158 for (Real t : times) {
159 for (Real strike : strikes) {
160 const Real xs = std::log(x: strike);
161
162 const Real expectedPDF = bsm.pdf(x: xs, t);
163 const Real calculatedPDF = heston.pdf(x: xs, t);
164
165 const Real tol = 1e-4;
166 if (std::fabs(x: expectedPDF - calculatedPDF) > tol) {
167 BOOST_FAIL("failed to reproduce Black-Scholes-Merton pdf "
168 "with the Heston model"
169 << "\n calculated: " << calculatedPDF
170 << "\n expected: " << expectedPDF
171 << "\n diff: " << calculatedPDF - expectedPDF
172 << "\n tol: " << tol);
173 }
174
175 const Real expectedCDF = bsm.cdf(x: xs, t);
176 const Real calculatedCDF = heston.cdf(x: xs, t);
177
178 if (std::fabs(x: expectedCDF - calculatedCDF) > tol) {
179 BOOST_FAIL("failed to reproduce Black-Scholes-Merton cdf "
180 "with the Heston model"
181 << "\n calculated: " << calculatedCDF
182 << "\n expected: " << expectedCDF
183 << "\n diff: " << calculatedCDF - expectedCDF
184 << "\n tol: " << tol);
185 }
186 }
187
188 for (Real prob : probs) {
189 const Real expectedInvCDF = bsm.invcdf(q: prob, t);
190 const Real calculatedInvCDF = heston.invcdf(q: prob, t);
191
192 const Real tol = 1e-3;
193 if (std::fabs(x: expectedInvCDF - calculatedInvCDF) > tol) {
194 BOOST_FAIL("failed to reproduce Black-Scholes-Merton "
195 "inverse cdf with the Heston model"
196 << "\n calculated: " << calculatedInvCDF
197 << "\n expected: " << expectedInvCDF
198 << "\n diff: " << calculatedInvCDF - expectedInvCDF
199 << "\n tol: " << tol);
200 }
201 }
202 }
203}
204
205namespace {
206 // see Svetlana Borovkova, Ferry J. Permana
207 // Implied volatility in oil markets
208 // http://www.researchgate.net/publication/46493859_Implied_volatility_in_oil_markets
209 class DumasParametricVolSurface : public BlackVolatilityTermStructure {
210 public:
211 DumasParametricVolSurface(Real b1,
212 Real b2,
213 Real b3,
214 Real b4,
215 Real b5,
216 ext::shared_ptr<Quote> spot,
217 const ext::shared_ptr<YieldTermStructure>& rTS,
218 ext::shared_ptr<YieldTermStructure> qTS)
219 : BlackVolatilityTermStructure(0, NullCalendar(), Following, rTS->dayCounter()), b1_(b1),
220 b2_(b2), b3_(b3), b4_(b4), b5_(b5), spot_(std::move(spot)), rTS_(rTS),
221 qTS_(std::move(qTS)) {}
222
223 Date maxDate() const override { return Date::maxDate(); }
224 Rate minStrike() const override { return 0.0; }
225 Rate maxStrike() const override { return QL_MAX_REAL; }
226
227 protected:
228 Volatility blackVolImpl(Time t, Real strike) const override {
229 QL_REQUIRE(t >= 0.0, "t must be >= 0");
230
231 if (t < QL_EPSILON)
232 return b1_;
233
234 const Real fwd = spot_->value()*qTS_->discount(t)/rTS_->discount(t);
235 const Real mn = std::log(x: fwd/strike)/std::sqrt(x: t);
236
237 return b1_ + b2_*mn + b3_*mn*mn + b4_*t + b5_*mn*t;
238 }
239
240 private:
241 const Real b1_, b2_, b3_, b4_, b5_;
242 const ext::shared_ptr<Quote> spot_;
243 const ext::shared_ptr<YieldTermStructure> rTS_;
244 const ext::shared_ptr<YieldTermStructure> qTS_;
245 };
246
247 class ProbWeightedPayoff {
248 public:
249 ProbWeightedPayoff(Time t,
250 ext::shared_ptr<Payoff> payoff,
251 ext::shared_ptr<RiskNeutralDensityCalculator> calc)
252 : t_(t), payoff_(std::move(payoff)), calc_(std::move(calc)) {}
253
254 Real operator()(Real x) const {
255 return calc_->pdf(x, t: t_) * (*payoff_)(std::exp(x: x));
256 }
257
258 private:
259 const Real t_;
260 const ext::shared_ptr<Payoff> payoff_;
261 const ext::shared_ptr<RiskNeutralDensityCalculator> calc_;
262 };
263
264 std::vector<Time> adaptiveTimeGrid(
265 Size maxStepsPerYear, Size minStepsPerYear, Real decay, Time endTime) {
266 const Time maxDt = 1.0/maxStepsPerYear;
267 const Time minDt = 1.0/minStepsPerYear;
268
269 Time t=0.0;
270 std::vector<Time> times(1, t);
271 while (t < endTime) {
272 const Time dt = maxDt*std::exp(x: -decay*t)
273 + minDt*(1.0-std::exp(x: -decay*t));
274 t+=dt;
275 times.push_back(x: std::min(a: endTime, b: t));
276 }
277
278 return times;
279 }
280}
281
282void RiskNeutralDensityCalculatorTest::testLocalVolatilityRND() {
283 BOOST_TEST_MESSAGE("Testing Fokker-Planck forward equation "
284 "for local volatility process to calculate "
285 "risk neutral densities...");
286
287 const DayCounter dayCounter = Actual365Fixed();
288 const Date todaysDate = Date(28, Dec, 2012);
289 Settings::instance().evaluationDate() = todaysDate;
290
291 const Rate r = 0.015;
292 const Rate q = 0.025;
293 const Real s0 = 100;
294 const Volatility v = 0.25;
295
296 const ext::shared_ptr<Quote> spot(
297 ext::make_shared<SimpleQuote>(args: s0));
298 const ext::shared_ptr<YieldTermStructure> rTS(
299 flatRate(today: todaysDate, forward: r, dc: dayCounter));
300 const ext::shared_ptr<YieldTermStructure> qTS(
301 flatRate(today: todaysDate, forward: q, dc: dayCounter));
302
303 const ext::shared_ptr<TimeGrid> timeGrid(new TimeGrid(1.0, 101));
304
305 const ext::shared_ptr<LocalVolRNDCalculator> constVolCalc(
306 new LocalVolRNDCalculator(
307 spot, rTS, qTS,
308 ext::make_shared<LocalConstantVol>(args: todaysDate, args: v, args: dayCounter),
309 timeGrid, 201));
310
311 const Real rTol = 0.01, atol = 0.005;
312 for (Time t=0.1; t < 0.99; t+=0.015) {
313 const Volatility stdDev = v * std::sqrt(x: t);
314 const Real xm = - 0.5 * stdDev * stdDev +
315 std::log(x: s0 * qTS->discount(t)/rTS->discount(t));
316
317 const GaussianDistribution gaussianPDF(xm, stdDev);
318 const CumulativeNormalDistribution gaussianCDF(xm, stdDev);
319 const InverseCumulativeNormal gaussianInvCDF(xm, stdDev);
320
321 for (Real x = xm - 3*stdDev; x < xm + 3*stdDev; x+=0.05) {
322 const Real expectedPDF = gaussianPDF(x);
323 const Real calculatedPDF = constVolCalc->pdf(x, t);
324 const Real absDiffPDF = std::fabs(x: expectedPDF - calculatedPDF);
325
326 if (absDiffPDF > atol || absDiffPDF/expectedPDF > rTol) {
327 BOOST_FAIL("failed to reproduce forward probability density"
328 << "\n time: " << t
329 << "\n spot " << std::exp(x)
330 << "\n calculated: " << calculatedPDF
331 << "\n expected: " << expectedPDF
332 << "\n abs diff: " << absDiffPDF
333 << "\n rel diff: " << absDiffPDF/expectedPDF
334 << "\n abs tol: " << atol
335 << "\n rel tol: " << rTol);
336 }
337
338 const Real expectedCDF = gaussianCDF(x);
339 const Real calculatedCDF = constVolCalc->cdf(x, t);
340 const Real absDiffCDF = std::fabs(x: expectedCDF - calculatedCDF);
341
342 if (absDiffCDF > atol) {
343 BOOST_FAIL("failed to reproduce forward "
344 "cumulative probability density"
345 << "\n time: " << t
346 << "\n spot " << std::exp(x)
347 << "\n calculated: " << calculatedCDF
348 << "\n expected: " << expectedCDF
349 << "\n abs diff: " << absDiffCDF
350 << "\n abs tol: " << atol);
351 }
352
353 const Real expectedX = x;
354 const Real calculatedX = constVolCalc->invcdf(p: expectedCDF, t);
355 const Real absDiffX = std::fabs(x: expectedX - calculatedX);
356
357 if (absDiffX > atol || absDiffX/expectedX > rTol) {
358 BOOST_FAIL("failed to reproduce "
359 "inverse cumulative probability density"
360 << "\n time: " << t
361 << "\n spot " << std::exp(x)
362 << "\n calculated: " << calculatedX
363 << "\n expected: " << expectedX
364 << "\n abs diff: " << absDiffX
365 << "\n abs tol: " << atol);
366 }
367 }
368 }
369
370 const Time tl = timeGrid->at(i: timeGrid->size()-5);
371 const Real xl = constVolCalc->mesher(t: tl)->locations().front();
372 if (!( constVolCalc->pdf(x: xl+0.0001, t: tl) > 0.0
373 && constVolCalc->pdf(x: xl-0.0001, t: tl) == 0.0)) {
374 BOOST_FAIL("probability outside interpolation range is not zero");
375 }
376
377 const Real b1 = 0.25;
378 const Real b2 = 0.03;
379 const Real b3 = 0.005;
380 const Real b4 = -0.02;
381 const Real b5 = -0.005;
382
383 const ext::shared_ptr<DumasParametricVolSurface> dumasVolSurface(
384 new DumasParametricVolSurface(b1, b2, b3, b4, b5, spot, rTS, qTS));
385
386 const ext::shared_ptr<BlackScholesMertonProcess> bsmProcess(
387 new BlackScholesMertonProcess(
388 Handle<Quote>(spot),
389 Handle<YieldTermStructure>(qTS),
390 Handle<YieldTermStructure>(rTS),
391 Handle<BlackVolTermStructure>(dumasVolSurface)));
392
393 const ext::shared_ptr<LocalVolTermStructure> localVolSurface
394 = ext::make_shared<NoExceptLocalVolSurface>(
395 args: Handle<BlackVolTermStructure>(dumasVolSurface),
396 args: Handle<YieldTermStructure>(rTS),
397 args: Handle<YieldTermStructure>(qTS),
398 args: Handle<Quote>(spot), args: b1);
399
400 const std::vector<Time> adaptiveGrid
401 = adaptiveTimeGrid(maxStepsPerYear: 400, minStepsPerYear: 50, decay: 5.0, endTime: 3.0);
402
403 const ext::shared_ptr<TimeGrid> dumasTimeGrid(
404 new TimeGrid(adaptiveGrid.begin(), adaptiveGrid.end()));
405
406 const ext::shared_ptr<LocalVolRNDCalculator> dumasVolCalc(
407 new LocalVolRNDCalculator(
408 spot, rTS, qTS, localVolSurface, dumasTimeGrid, 401, 0.1, 1e-8));
409
410 const Real strikes[] = { 25, 50, 95, 100, 105, 150, 200, 400 };
411 const std::vector<Date> maturities = {
412 todaysDate + Period(1, Weeks), todaysDate + Period(1, Months),
413 todaysDate + Period(3, Months), todaysDate + Period(6, Months),
414 todaysDate + Period(12, Months), todaysDate + Period(18, Months),
415 todaysDate + Period(2, Years), todaysDate + Period(3, Years) };
416
417
418 for (auto maturity : maturities) {
419 const Time expiry
420 = rTS->dayCounter().yearFraction(d1: todaysDate, d2: maturity);
421
422 const ext::shared_ptr<PricingEngine> engine(
423 new FdBlackScholesVanillaEngine(
424 bsmProcess, std::max(a: Size(51), b: Size(expiry*101)),
425 201, 0, FdmSchemeDesc::Douglas(), true, b1));
426
427 const ext::shared_ptr<Exercise> exercise(new EuropeanExercise(maturity));
428
429 for (Real strike : strikes) {
430 const ext::shared_ptr<StrikedTypePayoff> payoff(new PlainVanillaPayoff(
431 (strike > spot->value()) ? Option::Call : Option::Put, strike));
432
433 VanillaOption option(payoff, exercise);
434 option.setPricingEngine(engine);
435 const Real expected = option.NPV();
436
437 const Time tx = std::max(a: dumasTimeGrid->at(i: 1),
438 b: dumasTimeGrid->closestTime(t: expiry));
439 const std::vector<Real> x = dumasVolCalc->mesher(t: tx)->locations();
440
441 const ProbWeightedPayoff probWeightedPayoff(
442 expiry, payoff, dumasVolCalc);
443
444 const DiscountFactor df = rTS->discount(t: expiry);
445 const Real calculated = GaussLobattoIntegral(10000, 1e-10)(
446 probWeightedPayoff, x.front(), x.back()) * df;
447
448 const Real absDiff = std::fabs(x: expected - calculated);
449
450 if (absDiff > 0.5*atol) {
451 BOOST_ERROR("failed to reproduce option prices for"
452 << "\n expiry: " << expiry
453 << "\n strike: " << strike
454 << "\n expected: " << expected
455 << "\n calculated: " << calculated
456 << "\n diff: " << absDiff
457 << "\n abs tol: " << atol);
458 }
459 }
460 }
461}
462
463void RiskNeutralDensityCalculatorTest::testSquareRootProcessRND() {
464 BOOST_TEST_MESSAGE("Testing probability density for a square root process...");
465
466 struct SquareRootProcessParams {
467 const Real v0, kappa, theta, sigma;
468 };
469
470 const SquareRootProcessParams params[]
471 = { { .v0: 0.17, .kappa: 1.0, .theta: 0.09, .sigma: 0.5 },
472 { .v0: 1.0, .kappa: 0.6, .theta: 0.1, .sigma: 0.75 },
473 { .v0: 0.005, .kappa: 0.6, .theta: 0.1, .sigma: 0.05 } };
474
475 for (const auto& param : params) {
476 const SquareRootProcessRNDCalculator rndCalculator(param.v0, param.kappa, param.theta,
477 param.sigma);
478
479 const Time t = 0.75;
480 const Time tInfty = 60.0 / param.kappa;
481
482 const Real tol = 1e-10;
483 for (Real v = 1e-5; v < 1.0; v += (v < param.theta) ? 0.005 : 0.1) {
484
485 const Real cdfCalculated = rndCalculator.cdf(v, t);
486 const Real cdfExpected = GaussLobattoIntegral(10000, 0.01*tol)(
487 [&](Real _x) { return rndCalculator.pdf(v: _x, t); }, 0, v);
488
489 if (std::fabs(x: cdfCalculated - cdfExpected) > tol) {
490 BOOST_FAIL("failed to calculate cdf"
491 << "\n t: " << t
492 << "\n v: " << v
493 << "\n calculated: " << cdfCalculated
494 << "\n expected: " << cdfExpected
495 << "\n diff: " << cdfCalculated - cdfExpected
496 << "\n tolerance: " << tol);
497 }
498
499 if (cdfExpected < (1-1e-6) && cdfExpected > 1e-6) {
500 const Real vCalculated = rndCalculator.invcdf(q: cdfCalculated, t);
501
502 if (std::fabs(x: v - vCalculated) > tol) {
503 BOOST_FAIL("failed to calculate round trip cdf <-> invcdf"
504 << "\n t: " << t
505 << "\n v: " << v
506 << "\n cdf: " << cdfExpected
507 << "\n calculated: " << vCalculated
508 << "\n diff: " << v - vCalculated
509 << "\n tolerance: " << tol);
510 }
511 }
512
513 const Real statPdfCalculated = rndCalculator.pdf(v, t: tInfty);
514 const Real statPdfExpected = rndCalculator.stationary_pdf(v);
515
516 if (std::fabs(x: statPdfCalculated - statPdfExpected) > tol) {
517 BOOST_FAIL("failed to calculate stationary pdf"
518 << "\n v: " << v
519 << "\n calculated: " << statPdfCalculated
520 << "\n expected: " << statPdfExpected
521 << "\n diff: " << statPdfCalculated - statPdfExpected
522 << "\n tolerance: " << tol);
523 }
524
525 const Real statCdfCalculated = rndCalculator.cdf(v, t: tInfty);
526 const Real statCdfExpected = rndCalculator.stationary_cdf(v);
527
528 if (std::fabs(x: statCdfCalculated - statCdfExpected) > tol) {
529 BOOST_FAIL("failed to calculate stationary cdf"
530 << "\n v: " << v
531 << "\n calculated: " << statCdfCalculated
532 << "\n expected: " << statCdfExpected
533 << "\n diff: " << statCdfCalculated - statCdfExpected
534 << "\n tolerance: " << tol);
535 }
536 }
537
538 for (Real q = 1e-5; q < 1.0; q+=0.001) {
539 const Real statInvCdfCalculated = rndCalculator.invcdf(q, t: tInfty);
540 const Real statInvCdfExpected = rndCalculator.stationary_invcdf(q);
541
542 if (std::fabs(x: statInvCdfCalculated - statInvCdfExpected) > tol) {
543 BOOST_FAIL("failed to calculate stationary inverse of cdf"
544 << "\n q: " << q
545 << "\n calculated: " << statInvCdfCalculated
546 << "\n expected: " << statInvCdfExpected
547 << "\n diff: " << statInvCdfCalculated - statInvCdfExpected
548 << "\n tolerance: " << tol);
549 }
550 }
551 }
552}
553
554void RiskNeutralDensityCalculatorTest::testBlackScholesWithSkew() {
555 BOOST_TEST_MESSAGE(
556 "Testing probability density for a BSM process "
557 "with strike dependent volatility vs local volatility...");
558
559 const Date todaysDate = Date(3, Oct, 2016);
560 Settings::instance().evaluationDate() = todaysDate;
561
562 const DayCounter dc = Actual365Fixed();
563 const Date maturityDate = todaysDate + Period(3, Months);
564 const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
565
566 // use Heston model to create volatility surface with skew
567 const Real r = 0.08;
568 const Real q = 0.03;
569 const Real s0 = 100;
570 const Real v0 = 0.06;
571 const Real kappa = 1.0;
572 const Real theta = 0.06;
573 const Real sigma = 0.4;
574 const Real rho = -0.75;
575
576 const Handle<YieldTermStructure> rTS(flatRate(today: todaysDate, forward: r, dc));
577 const Handle<YieldTermStructure> qTS(flatRate(today: todaysDate, forward: q, dc));
578 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
579
580 const ext::shared_ptr<HestonProcess> hestonProcess(
581 ext::make_shared<HestonProcess>(
582 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho));
583
584 const Handle<BlackVolTermStructure> hestonSurface(
585 ext::make_shared<HestonBlackVolSurface>(
586 args: Handle<HestonModel>(ext::make_shared<HestonModel>(args: hestonProcess)),
587 args: AnalyticHestonEngine::AndersenPiterbarg,
588 args: AnalyticHestonEngine::Integration::discreteTrapezoid(evaluation: 64)));
589
590 const ext::shared_ptr<TimeGrid> timeGrid(new TimeGrid(maturity, 51));
591
592 const ext::shared_ptr<LocalVolTermStructure> localVol(
593 ext::make_shared<NoExceptLocalVolSurface>(
594 args: hestonSurface, args: rTS, args: qTS, args: spot, args: std::sqrt(x: theta)));
595
596 const LocalVolRNDCalculator localVolCalc(
597 spot.currentLink(), rTS.currentLink(), qTS.currentLink(), localVol,
598 timeGrid, 151, 0.25);
599
600 const HestonRNDCalculator hestonCalc(hestonProcess);
601
602 const GBSMRNDCalculator gbsmCalc(
603 ext::make_shared<BlackScholesMertonProcess>(
604 args: spot, args: qTS, args: rTS, args: hestonSurface));
605
606 const Real strikes[] = { 85, 75, 90, 110, 125, 150 };
607
608 for (Real strike : strikes) {
609 const Real logStrike = std::log(x: strike);
610
611 const Real expected = hestonCalc.cdf(x: logStrike, t: maturity);
612 const Real calculatedGBSM = gbsmCalc.cdf(s: strike, t: maturity);
613
614 const Real gbsmTol = 1e-5;
615 if (std::fabs(x: expected - calculatedGBSM) > gbsmTol) {
616 BOOST_FAIL("failed to match Heston and GBSM cdf"
617 << "\n t: " << maturity
618 << "\n k: " << strike
619 << "\n calculated: " << calculatedGBSM
620 << "\n expected: " << expected
621 << "\n diff: " <<
622 std::fabs(calculatedGBSM - expected)
623 << "\n tolerance: " << gbsmTol);
624 }
625
626 const Real calculatedLocalVol = localVolCalc.cdf(x: logStrike, t: maturity);
627 const Real localVolTol = 1e-3;
628 if (std::fabs(x: expected - calculatedLocalVol) > localVolTol) {
629 BOOST_FAIL("failed to match Heston and local Volatility cdf"
630 << "\n t: " << maturity
631 << "\n k: " << strike
632 << "\n calculated: " << calculatedLocalVol
633 << "\n expected: " << expected
634 << "\n diff: " <<
635 std::fabs(calculatedLocalVol - expected)
636 << "\n tolerance: " << localVolTol);
637 }
638 }
639
640 for (Real strike : strikes) {
641 const Real logStrike = std::log(x: strike);
642
643 const Real expected = hestonCalc.pdf(x: logStrike, t: maturity)/strike;
644 const Real calculatedGBSM = gbsmCalc.pdf(s: strike, t: maturity);
645
646 const Real gbsmTol = 1e-5;
647 if (std::fabs(x: expected - calculatedGBSM) > gbsmTol) {
648 BOOST_FAIL("failed to match Heston and GBSM pdf"
649 << "\n t: " << maturity
650 << "\n k: " << strike
651 << "\n calculated: " << calculatedGBSM
652 << "\n expected: " << expected
653 << "\n diff: " <<
654 std::fabs(calculatedGBSM - expected)
655 << "\n tolerance: " << gbsmTol);
656 }
657
658 const Real calculatedLocalVol
659 = localVolCalc.pdf(x: logStrike, t: maturity)/strike;
660 const Real localVolTol = 1e-4;
661 if (std::fabs(x: expected - calculatedLocalVol) > localVolTol) {
662 BOOST_FAIL("failed to match Heston and local Volatility pdf"
663 << "\n t: " << maturity
664 << "\n k: " << strike
665 << "\n calculated: " << calculatedLocalVol
666 << "\n expected: " << expected
667 << "\n diff: " <<
668 std::fabs(calculatedLocalVol - expected)
669 << "\n tolerance: " << localVolTol);
670 }
671 }
672
673 const Real quantiles[] = { 0.05, 0.25, 0.5, 0.75, 0.95 };
674 for (Real quantile : quantiles) {
675 const Real expected = std::exp(x: hestonCalc.invcdf(q: quantile, t: maturity));
676 const Real calculatedGBSM = gbsmCalc.invcdf(q: quantile, t: maturity);
677
678 const Real gbsmTol = 1e-3;
679 if (std::fabs(x: expected - calculatedGBSM) > gbsmTol) {
680 BOOST_FAIL("failed to match Heston and GBSM invcdf"
681 << "\n t: " << maturity
682 << "\n quantile: " << quantile
683 << "\n calculated: " << calculatedGBSM
684 << "\n expected: " << expected
685 << "\n diff: " <<
686 std::fabs(calculatedGBSM - expected)
687 << "\n tolerance: " << gbsmTol);
688 }
689
690 const Real calculatedLocalVol
691 = std::exp(x: localVolCalc.invcdf(p: quantile, t: maturity));
692 const Real localVolTol = 0.1;
693 if (std::fabs(x: expected - calculatedLocalVol) > localVolTol) {
694 BOOST_FAIL("failed to match Heston and local Volatility invcdf"
695 << "\n t: " << maturity
696 << "\n k: " << quantile
697 << "\n calculated: " << calculatedLocalVol
698 << "\n expected: " << expected
699 << "\n diff: " <<
700 std::fabs(calculatedLocalVol - expected)
701 << "\n tolerance: " << localVolTol);
702 }
703 }
704}
705
706void RiskNeutralDensityCalculatorTest::testMassAtZeroCEVProcessRND() {
707 BOOST_TEST_MESSAGE("Testing the mass at zero for a "
708 "constant elasticity of variance (CEV) process...");
709
710 const Real f0 = 100.0;
711 const Time t = 2.75;
712
713 const std::pair<Real, Real> params[] = {
714 {0.1, 1.6},
715 {0.01, 2.0},
716 {10.0, 0.35},
717 {50.0, 0.1}
718 };
719
720 const Real tol = 1e-4;
721
722 for (const auto& param : params) {
723 const Real alpha = param.first;
724 const Real beta = param.second;
725
726 const ext::shared_ptr<CEVRNDCalculator> calculator =
727 ext::make_shared<CEVRNDCalculator>(args: f0, args: alpha, args: beta);
728
729 const Real ax = 15.0*std::sqrt(x: t)*alpha*std::pow(x: f0, y: beta);
730
731 const Real calculated = GaussLobattoIntegral(1000, 1e-8)(
732 [&](Real _x) -> Real { return calculator->pdf(f: _x, t); }, std::max(QL_EPSILON, b: f0-ax), f0+ax) +
733 calculator->massAtZero(t);
734
735 if (std::fabs(x: calculated - 1.0) > tol) {
736 BOOST_FAIL("failed to reproduce the total probability mass"
737 << "\n alpha: " << alpha
738 << "\n beta: " << beta
739 << "\n prob mass: " << calculated
740 << "\n tolerance: " << tol);
741 }
742 }
743}
744
745void RiskNeutralDensityCalculatorTest::testCEVCDF() {
746 BOOST_TEST_MESSAGE("Testing CDF for a "
747 "constant elasticity of variance (CEV) process...");
748
749 const Real f0 = 2.1;
750 const Time t = 0.75;
751
752 const Real alpha = 0.1;
753 const Real betas[] = { 0.45, 1.25 };
754
755 const Real tol = 1e-6;
756 for (Size i = 1; i < LENGTH(betas); ++i) {
757 const Real beta = betas[i];
758 const ext::shared_ptr<CEVRNDCalculator> calculator =
759 ext::make_shared<CEVRNDCalculator>(args: f0, args: alpha, args: beta);
760
761 for (Real x = 1.3; x < 3.1; x+=0.1) {
762
763 const Real cdfValue = calculator->cdf(f: x, t);
764 const Real calculated = calculator->invcdf(q: cdfValue, t);
765
766 if (std::fabs(x: x - calculated) > tol) {
767 BOOST_FAIL(
768 "failed to reproduce the inverse cumulative probability"
769 << "\n alpha: " << alpha
770 << "\n beta: " << beta
771 << "\n x: " << x
772 << "\n calculated:" << calculated
773 << "\n difference:" << x - calculated
774 << "\n tolerance: " << tol);
775 }
776 }
777 }
778}
779
780test_suite* RiskNeutralDensityCalculatorTest::experimental(SpeedLevel speed) {
781 auto* suite = BOOST_TEST_SUITE("Risk neutral density calculator tests");
782
783 suite->add(QUANTLIB_TEST_CASE(
784 &RiskNeutralDensityCalculatorTest::testDensityAgainstOptionPrices));
785 suite->add(QUANTLIB_TEST_CASE(
786 &RiskNeutralDensityCalculatorTest::testBSMagainstHestonRND));
787 suite->add(QUANTLIB_TEST_CASE(
788 &RiskNeutralDensityCalculatorTest::testLocalVolatilityRND));
789 suite->add(QUANTLIB_TEST_CASE(
790 &RiskNeutralDensityCalculatorTest::testSquareRootProcessRND));
791 suite->add(QUANTLIB_TEST_CASE(
792 &RiskNeutralDensityCalculatorTest::testMassAtZeroCEVProcessRND));
793 suite->add(QUANTLIB_TEST_CASE(
794 &RiskNeutralDensityCalculatorTest::testCEVCDF));
795
796 if (speed <= Fast) {
797 suite->add(QUANTLIB_TEST_CASE(
798 &RiskNeutralDensityCalculatorTest::testBlackScholesWithSkew));
799 }
800
801 return suite;
802}
803

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