1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2019 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 "fdsabr.hpp"
21#include "utilities.hpp"
22#include <ql/functional.hpp>
23#include <ql/instruments/vanillaoption.hpp>
24#include <ql/math/comparison.hpp>
25#include <ql/math/randomnumbers/rngtraits.hpp>
26#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
27#include <ql/math/richardsonextrapolation.hpp>
28#include <ql/math/statistics/generalstatistics.hpp>
29#include <ql/methods/finitedifferences/utilities/cevrndcalculator.hpp>
30#include <ql/pricingengines/vanilla/analyticcevengine.hpp>
31#include <ql/pricingengines/vanilla/fdsabrvanillaengine.hpp>
32#include <ql/processes/blackscholesprocess.hpp>
33#include <ql/quotes/simplequote.hpp>
34#include <ql/termstructures/volatility/sabr.hpp>
35#include <ql/shared_ptr.hpp>
36#include <utility>
37
38using namespace QuantLib;
39using boost::unit_test_framework::test_suite;
40
41namespace {
42 class SabrMonteCarloPricer {
43 public:
44 SabrMonteCarloPricer(Real f0,
45 Time maturity,
46 ext::shared_ptr<Payoff> payoff,
47 Real alpha,
48 Real beta,
49 Real nu,
50 Real rho)
51 : f0_(f0), maturity_(maturity), payoff_(std::move(payoff)), alpha_(alpha), beta_(beta),
52 nu_(nu), rho_(rho) {}
53
54 Real operator()(Real dt) const {
55 const Size nSims = 64*1024;
56
57 const Real timeStepsPerYear = 1./dt;
58 const Size timeSteps = Size(maturity_*timeStepsPerYear+1e-8);
59
60 const Real sqrtDt = std::sqrt(x: dt);
61 const Real w = std::sqrt(x: 1.0-rho_*rho_);
62
63 const Real logAlpha = std::log(x: alpha_);
64
65 SobolBrownianBridgeRsg rsg(2, timeSteps, SobolBrownianGenerator::Diagonal, 12345U);
66
67 GeneralStatistics stats;
68
69 for (Size i=0; i < nSims; ++i) {
70 Real f = f0_;
71 Real a = logAlpha;
72
73 const std::vector<Real> n = rsg.nextSequence().value;
74
75 for (Size j=0; j < timeSteps && f > 0.0; ++j) {
76
77 const Real r1 = n[j];
78 const Real r2 = rho_*r1 + n[j+timeSteps]*w;
79
80 //Sample CEV distribution: accurate but slow
81 //
82 //const CEVRNDCalculator calc(f, std::exp(a), beta_);
83 //const Real u = CumulativeNormalDistribution()(r1);
84 //f = calc.invcdf(u, dt);
85
86 // simple Euler method
87 f += std::exp(x: a)*std::pow(x: f, y: beta_)*r1*sqrtDt;
88 a += - 0.5*nu_*nu_*dt + nu_*r2*sqrtDt;
89 }
90 f = std::max(a: 0.0, b: f);
91 stats.add(value: (*payoff_)(f));
92 }
93
94 return stats.mean();
95 }
96
97 private:
98 const Real f0_;
99 const Time maturity_;
100 const ext::shared_ptr<Payoff> payoff_;
101 const Real alpha_, beta_, nu_, rho_;
102 };
103
104}
105
106
107void FdSabrTest::testFdmSabrOp() {
108 BOOST_TEST_MESSAGE("Testing FDM SABR operator...");
109
110 const Date today = Date(22, February, 2018);
111 const DayCounter dc = Actual365Fixed();
112 Settings::instance().evaluationDate() = today;
113
114 const Date maturityDate = today + Period(2, Years);
115 const Time maturityTime = dc.yearFraction(d1: today, d2: maturityDate);
116
117 const Real strike = 1.5;
118
119 const ext::shared_ptr<Exercise> exercise =
120 ext::make_shared<EuropeanExercise>(args: maturityDate);
121
122 const ext::shared_ptr<PlainVanillaPayoff> putPayoff =
123 ext::make_shared<PlainVanillaPayoff>(args: Option::Put, args: strike);
124 const ext::shared_ptr<PlainVanillaPayoff> callPayoff =
125 ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike);
126
127 VanillaOption optionPut(putPayoff, exercise);
128 VanillaOption optionCall(callPayoff, exercise);
129
130 const Handle<YieldTermStructure> rTS =
131 Handle<YieldTermStructure>(flatRate(today, forward: 0.0, dc));
132
133 const Real f0 = 1.0;
134 const Real alpha = 0.35;
135 const Real nu = 1.0;
136 const Real rho = 0.25;
137
138 const Real betas[] = { 0.25, 0.6 };
139
140 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess =
141 ext::make_shared<GeneralizedBlackScholesProcess>(
142 args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: f0)),
143 args: rTS, args: rTS, args: Handle<BlackVolTermStructure>(flatVol(volatility: 0.2, dc)));
144
145 for (Real beta : betas) {
146
147 const ext::shared_ptr<PricingEngine> pdeEngine =
148 ext::make_shared<FdSabrVanillaEngine>(args: f0, args: alpha, args&: beta, args: nu, args: rho, args: rTS, args: 100, args: 400, args: 100);
149
150 optionPut.setPricingEngine(pdeEngine);
151 const Real pdePut = optionPut.NPV();
152
153 // check put/call parity
154 optionCall.setPricingEngine(pdeEngine);
155 const Real pdeCall = optionCall.NPV();
156
157 const Real pdeFwd = pdeCall - pdePut;
158
159 const Real parityDiff = std::fabs(x: pdeFwd - (f0 - strike));
160 const Real parityTol = 1e-4;
161 if (parityDiff > parityTol) {
162 BOOST_ERROR(
163 "failed to validate the call/put parity"
164 << "\n beta : " << beta
165 << "\n strike : " << strike
166 << "\n fwd (call/put) : " << pdeFwd
167 << "\n fwd (f0-strike): " << f0-strike
168 << "\n diff : " << parityDiff
169 << "\n tol : " << parityTol);
170 }
171
172 const Real putPdeImplVol =
173 optionPut.impliedVolatility(price: optionPut.NPV(), process: bsProcess, accuracy: 1e-6);
174
175 const ext::function<Real(Real)> mcSabr(
176 SabrMonteCarloPricer(f0, maturityTime, putPayoff,
177 alpha, beta, nu, rho));
178
179 const Real mcNPV = RichardsonExtrapolation(
180 mcSabr, 1/4.0)(4.0, 2.0);
181
182 const Real putMcImplVol =
183 optionPut.impliedVolatility(price: mcNPV, process: bsProcess, accuracy: 1e-6);
184
185 const Real volDiff = std::fabs(x: putPdeImplVol - putMcImplVol);
186
187 const Real volTol = 5e-3;
188 if (volDiff > volTol) {
189 BOOST_ERROR(
190 "failed to validate PDE against MC implied volatility"
191 << "\n beta : " << beta
192 << "\n strike : " << strike
193 << "\n PDE impl vol : " << putPdeImplVol
194 << "\n MC impl vol : " << putMcImplVol
195 << "\n diff : " << volDiff
196 << "\n tol : " << volTol);
197 }
198 }
199}
200
201void FdSabrTest::testFdmSabrCevPricing() {
202 BOOST_TEST_MESSAGE("Testing FDM CEV pricing with trivial SABR model...");
203
204 const Date today = Date(3, January, 2019);
205 const DayCounter dc = Actual365Fixed();
206 Settings::instance().evaluationDate() = today;
207
208 const Date maturityDate = today + Period(12, Months);
209
210 const Real betas[] = { 0.1, 0.9 };
211 const Real strikes[] = { 0.9, 1.5 };
212
213 const Real f0 = 1.2;
214 const Real alpha = 0.35;
215 const Real nu = 1e-3;
216 const Real rho = 0.25;
217
218 const Handle<YieldTermStructure> rTS = Handle<YieldTermStructure>(
219 flatRate(today, forward: 0.05, dc));
220
221 const ext::shared_ptr<Exercise> exercise =
222 ext::make_shared<EuropeanExercise>(args: maturityDate);
223
224 const Option::Type optionTypes[] = {Option::Put, Option::Call};
225
226 const Real tol = 5e-5;
227
228 for (auto optionType : optionTypes) {
229 for (Real strike : strikes) {
230 const ext::shared_ptr<PlainVanillaPayoff> payoff =
231 ext::make_shared<PlainVanillaPayoff>(args&: optionType, args&: strike);
232
233 VanillaOption option(payoff, exercise);
234
235 for (Real beta : betas) {
236 option.setPricingEngine(ext::make_shared<FdSabrVanillaEngine>(
237 args: f0, args: alpha, args&: beta, args: nu, args: rho, args: rTS, args: 100, args: 400, args: 3));
238
239 const Real calculated = option.NPV();
240
241 option.setPricingEngine(ext::make_shared<AnalyticCEVEngine>(
242 args: f0, args: alpha, args&: beta, args: rTS));
243
244 const Real expected = option.NPV();
245
246 if (std::fabs(x: expected-calculated) > tol) {
247 BOOST_ERROR(
248 "failed to calculate vanilla CEV option prices"
249 << "\n beta : " << beta
250 << "\n strike : " << strike
251 << "\n option type : "
252 << ((payoff->optionType() == Option::Call) ? "Call" : "Put")
253 << "\n analytic npv : " << expected
254 << "\n pde npv : " << calculated
255 << "\n npv difference : "
256 << std::fabs(expected - calculated)
257 << "\n tolerance : " << tol);
258 }
259 }
260 }
261 }
262}
263
264void FdSabrTest::testFdmSabrVsVolApproximation() {
265 BOOST_TEST_MESSAGE("Testing FDM SABR vs approximations...");
266
267 const Date today = Date(8, January, 2019);
268 const DayCounter dc = Actual365Fixed();
269 Settings::instance().evaluationDate() = today;
270
271 const Date maturityDate = today + Period(6, Months);
272 const Time maturityTime = dc.yearFraction(d1: today, d2: maturityDate);
273
274 const Handle<YieldTermStructure> rTS = Handle<YieldTermStructure>(
275 flatRate(today, forward: 0.05, dc));
276
277 const Real f0 = 100;
278
279 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess =
280 ext::make_shared<GeneralizedBlackScholesProcess>(
281 args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: f0)),
282 args: rTS, args: rTS, args: Handle<BlackVolTermStructure>(flatVol(volatility: 0.2, dc)));
283
284 const Real alpha = 0.35;
285 const Real beta = 0.85;
286 const Real nu = 0.75;
287 const Real rho = 0.85;
288
289 const Real strikes[] = { 90, 100, 110};
290 const Option::Type optionTypes[] = {Option::Put, Option::Call};
291
292 const Real tol = 2.5e-3;
293 for (auto optionType : optionTypes) {
294 for (Real strike : strikes) {
295 VanillaOption option(ext::make_shared<PlainVanillaPayoff>(args&: optionType, args&: strike),
296 ext::make_shared<EuropeanExercise>(args: maturityDate));
297
298 option.setPricingEngine(ext::make_shared<FdSabrVanillaEngine>(
299 args: f0, args: alpha, args: beta, args: nu, args: rho, args: rTS, args: 25, args: 100, args: 50));
300
301 const Volatility fdmVol =
302 option.impliedVolatility(price: option.NPV(), process: bsProcess);
303
304 const Real hagenVol = sabrVolatility(
305 strike, forward: f0, expiryTime: maturityTime, alpha, beta, nu, rho);
306
307 const Real diff = std::fabs(x: fdmVol - hagenVol);
308
309 if (std::fabs(x: fdmVol-hagenVol) > tol) {
310 BOOST_ERROR(
311 "large difference between Hagen formula and FDM"
312 << "\n strike : " << strike
313 << "\n option type : "
314 << ((optionType == Option::Call) ? "Call" : "Put")
315 << "\n Hagen vol : " << hagenVol
316 << "\n pde vol : " << fdmVol
317 << "\n vol difference : " << diff
318 << "\n tolerance : " << tol);
319 }
320 }
321 }
322}
323
324
325namespace {
326 /*
327 * Example and reference values are taken from
328 * B. Chen, C.W. Oosterlee, H. Weide,
329 * Efficient unbiased simulation scheme for the SABR stochastic volatility model.
330 * https://http://ta.twi.tudelft.nl/mf/users/oosterle/oosterlee/SABRMC.pdf
331 */
332
333 class OsterleeReferenceResults {
334 public:
335 explicit OsterleeReferenceResults(Size i) : i_(i) { }
336
337 Real operator()(Real t) const {
338 Size i;
339 if (close_enough(x: t, y: 1/16.))
340 i = 0;
341 else if (close_enough(x: t, y: 1/32.))
342 i = 1;
343 else
344 QL_FAIL("unmatched reference result lookup");
345
346 return data_[i_][i];
347 }
348
349 private:
350 const Size i_;
351 static Real data_[9][3];
352 };
353
354 Real OsterleeReferenceResults::data_[9][3] = {
355 { 0.0610, 0.0604 }, { 0.0468, 0.0463 }, { 0.0347, 0.0343 },
356 { 0.0632, 0.0625 }, { 0.0512, 0.0506 }, { 0.0406, 0.0400 },
357 { 0.0635, 0.0630 }, { 0.0523, 0.0520 }, { 0.0422, 0.0421 }
358 };
359}
360
361void FdSabrTest::testOosterleeTestCaseIV() {
362 BOOST_TEST_MESSAGE("Testing Chen, Oosterlee and Weide test case IV...");
363
364 const Date today = Date(8, January, 2019);
365 const DayCounter dc = Actual365Fixed();
366 Settings::instance().evaluationDate() = today;
367
368 const Handle<YieldTermStructure> rTS =
369 Handle<YieldTermStructure>(flatRate(today, forward: 0.0, dc));
370
371 const Real f0 = 0.07;
372 const Real alpha = 0.4;
373 const Real nu = 0.8;
374 const Real beta = 0.4;
375 const Real rho = -0.6;
376
377 const Period maturities[] = {
378 Period(2, Years), Period(5, Years), Period(10, Years)
379 };
380
381 const Real strikes[] = { 0.4*f0, f0, 1.6*f0 };
382
383 const Real tol = 0.00035;
384 for (Size i=0; i < LENGTH(maturities); ++i) {
385 const Date maturityDate = today + maturities[i];
386 const Time maturityTime = dc.yearFraction(d1: today, d2: maturityDate);
387
388 const Size timeSteps = Size(5*maturityTime);
389
390 const ext::shared_ptr<PricingEngine> engine =
391 ext::make_shared<FdSabrVanillaEngine>(
392 args: f0, args: alpha, args: beta, args: nu, args: rho, args: rTS, args: timeSteps, args: 200, args: 21);
393
394 const ext::shared_ptr<Exercise> exercise =
395 ext::make_shared<EuropeanExercise>(args: maturityDate);
396
397 for (Size j=0; j < LENGTH(strikes); ++j) {
398 const ext::shared_ptr<StrikedTypePayoff> payoff =
399 ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strikes[j]);
400
401 VanillaOption option(payoff, exercise);
402 option.setPricingEngine(engine);
403
404 const Real calculated = option.NPV();
405
406 const OsterleeReferenceResults referenceResuts(i*3+j);
407
408 const Real expected = RichardsonExtrapolation(
409 ext::function<Real(Real)>(referenceResuts), 1/16., 1)(2.);
410
411 const Real diff = std::fabs(x: calculated - expected);
412 if (diff > tol) {
413 BOOST_ERROR(
414 "can not reproduce reference values from Monte-Carlo"
415 << "\n strike : " << payoff->strike()
416 << "\n maturity : " << maturityDate
417 << "\n reference : " << expected
418 << "\n calculated : " << calculated
419 << "\n difference : " << diff
420 << "\n tolerance : " << tol);
421 }
422 }
423 }
424}
425
426void FdSabrTest::testBenchOpSabrCase() {
427 BOOST_TEST_MESSAGE("Testing SABR BenchOp problem...");
428
429 /*
430 * von Sydow, L, Milovanović, S, Larsson, E, In't Hout, K,
431 * Wiktorsson, M, Oosterlee, C.W, Shcherbakov, V, Wyns, M,
432 * Leitao Rodriguez, A, Jain, S, et al. (2018)
433 * BENCHOP–SLV: the BENCHmarking project in Option
434 * Pricing–Stochastic and Local Volatility problems
435 * https://ir.cwi.nl/pub/28249
436 */
437
438 const Date today = Date(8, January, 2019);
439 const DayCounter dc = Actual365Fixed();
440 Settings::instance().evaluationDate() = today;
441
442 const Handle<YieldTermStructure> rTS =
443 Handle<YieldTermStructure>(flatRate(today, forward: 0.0, dc));
444
445 const Size maturityInYears[] = { 2, 10 };
446
447 const Real f0s[] = { 0.5, 0.07 };
448 const Real alphas[] = { 0.5, 0.4 };
449 const Real nus[] = { 0.4, 0.8 };
450 const Real betas[] = { 0.5, 0.5 };
451 const Real rhos[] = { 0.0, -0.6 };
452
453 const Real expected[2][3] = {
454 { 0.221383196830866, 0.193836689413803, 0.166240814653231 },
455 { 0.052450313614407, 0.046585753491306, 0.039291470612989 }
456 };
457
458 const Size gridX = 400;
459 const Size gridY = 25;
460 const Size gridT = 10;
461
462 const Real factor = 2;
463
464 const Real tol = 2e-4;
465
466 for (Size i=0; i < LENGTH(f0s); ++i) {
467
468 const Date maturity = today + Period(maturityInYears[i]*365, Days);
469 const Time T = dc.yearFraction(d1: today, d2: maturity);
470
471 const Real f0 = f0s[i];
472 const Real alpha = alphas[i];
473 const Real nu = nus[i];
474 const Real beta = betas[i];
475 const Real rho = rhos[i];
476
477 const Real strikes[] = {
478 f0*std::exp(x: -0.1*std::sqrt(x: T)), f0, f0*std::exp(x: 0.1*std::sqrt(x: T))
479 };
480
481 for (Size j=0; j < LENGTH(strikes); ++j) {
482 const Real strike = strikes[j];
483
484 VanillaOption option(
485 ext::make_shared<PlainVanillaPayoff>(args: Option::Call, args: strike),
486 ext::make_shared<EuropeanExercise>(args: maturity));
487
488 option.setPricingEngine(ext::make_shared<FdSabrVanillaEngine>(
489 args: f0, args: alpha, args: beta, args: nu, args: rho, args: rTS,
490 args: Size(gridT*factor),
491 args: Size(gridX*factor),
492 args: Size(gridY*std::sqrt(x: factor))));
493
494 const Real calculated = option.NPV();
495 const Real diff = std::fabs(x: calculated - expected[i][j]);
496
497 if (diff > tol) {
498 BOOST_ERROR(
499 "failed to reproduce reference values"
500 << "\n strike : " << strike
501 << "\n maturity : " << maturity
502 << "\n reference : " << expected[i][j]
503 << "\n calculated : " << calculated
504 << "\n difference : " << diff
505 << "\n tolerance : " << tol);
506 }
507 }
508 }
509}
510
511test_suite* FdSabrTest::suite(SpeedLevel speed) {
512 auto* suite = BOOST_TEST_SUITE("Finite Difference SABR tests");
513
514 suite->add(QUANTLIB_TEST_CASE(&FdSabrTest::testFdmSabrCevPricing));
515 suite->add(QUANTLIB_TEST_CASE(&FdSabrTest::testFdmSabrVsVolApproximation));
516 suite->add(QUANTLIB_TEST_CASE(&FdSabrTest::testOosterleeTestCaseIV));
517 suite->add(QUANTLIB_TEST_CASE(&FdSabrTest::testBenchOpSabrCase));
518
519 if (speed <= Fast) {
520 suite->add(QUANTLIB_TEST_CASE(&FdSabrTest::testFdmSabrOp));
521 }
522
523 return suite;
524}
525

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