[go: up one dir, main page]
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2017 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 "utilities.hpp"
22#include "squarerootclvmodel.hpp"
23#include <ql/quotes/simplequote.hpp>
24#include <ql/time/daycounters/actualactual.hpp>
25#include <ql/time/daycounters/actual365fixed.hpp>
26#include <ql/time/calendars/nullcalendar.hpp>
27#include <ql/instruments/impliedvolatility.hpp>
28#include <ql/instruments/forwardvanillaoption.hpp>
29#include <ql/math/statistics/statistics.hpp>
30#include <ql/math/integrals/gausslobattointegral.hpp>
31#include <ql/math/randomnumbers/rngtraits.hpp>
32#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
33#include <ql/math/optimization/constraint.hpp>
34#include <ql/math/optimization/simplex.hpp>
35#include <ql/processes/squarerootprocess.hpp>
36#include <ql/methods/montecarlo/multipathgenerator.hpp>
37#include <ql/pricingengines/blackcalculator.hpp>
38#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
39#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
40#include <ql/pricingengines/forward/forwardengine.hpp>
41#include <ql/methods/montecarlo/pathgenerator.hpp>
42#include <ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp>
43#include <ql/termstructures/volatility/equityfx/noexceptlocalvolsurface.hpp>
44#include <ql/experimental/models/squarerootclvmodel.hpp>
45#include <ql/models/equity/hestonslvfdmmodel.hpp>
46#include <ql/processes/hestonslvprocess.hpp>
47#include <ql/pricingengines/barrier/fdhestondoublebarrierengine.hpp>
48#include <ql/pricingengines/barrier/analyticdoublebarrierbinaryengine.hpp>
49#include <ql/experimental/volatility/sabrvoltermstructure.hpp>
50
51#include <boost/math/distributions/non_central_chi_squared.hpp>
52
53#include <set>
54#include <utility>
55
56using namespace QuantLib;
57using namespace boost::unit_test_framework;
58
59
60namespace square_root_clv_model {
61 class CLVModelPayoff : public PlainVanillaPayoff {
62 public:
63 CLVModelPayoff(Option::Type type, Real strike, ext::function<Real(Real)> g)
64 : PlainVanillaPayoff(type, strike), g_(std::move(g)) {}
65
66 Real operator()(Real x) const override { return PlainVanillaPayoff::operator()(price: g_(x)); }
67
68 private:
69 const ext::function<Real(Real)> g_;
70 };
71
72 typedef boost::math::non_central_chi_squared_distribution<Real>
73 chi_squared_type;
74}
75
76
77void SquareRootCLVModelTest::testSquareRootCLVVanillaPricing() {
78 BOOST_TEST_MESSAGE(
79 "Testing vanilla option pricing with square-root kernel process...");
80
81 using namespace square_root_clv_model;
82
83 const Date todaysDate(5, Oct, 2016);
84 Settings::instance().evaluationDate() = todaysDate;
85
86 const DayCounter dc = ActualActual(ActualActual::ISDA);
87 const Date maturityDate = todaysDate + Period(3, Months);
88 const Time maturity = dc.yearFraction(d1: todaysDate, d2: maturityDate);
89
90 const Real s0 = 100;
91 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
92
93 const Rate r = 0.08;
94 const Rate q = 0.03;
95 const Volatility vol = 0.3;
96
97 const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
98 const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
99 const Handle<BlackVolTermStructure> volTS(flatVol(today: todaysDate, volatility: vol, dc));
100 const Real fwd = s0*qTS->discount(t: maturity)/rTS->discount(t: maturity);
101
102 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
103 ext::make_shared<GeneralizedBlackScholesProcess>(
104 args: spot, args: qTS, args: rTS, args: volTS));
105
106 const Real kappa = 1.0;
107 const Real theta = 0.06;
108 const Volatility sigma = 0.2;
109 const Real x0 = 0.09;
110
111 const ext::shared_ptr<SquareRootProcess> sqrtProcess(
112 ext::make_shared<SquareRootProcess>(args: theta, args: kappa, args: sigma, args: x0));
113
114 const std::vector<Date> maturityDates(1, maturityDate);
115
116 const SquareRootCLVModel model(
117 bsProcess, sqrtProcess, maturityDates, 14, 1-1e-14, 1e-14);
118
119 const Array x = model.collocationPointsX(d: maturityDate);
120 const Array y = model.collocationPointsY(d: maturityDate);
121
122 const LagrangeInterpolation g(x.begin(), x.end(), y.begin());
123
124 const Real df = 4*theta*kappa/(sigma*sigma);
125 const Real ncp = 4*kappa*std::exp(x: -kappa*maturity)
126 / (sigma*sigma*(1-std::exp(x: -kappa*maturity)))*sqrtProcess->x0();
127
128 const chi_squared_type dist(df, ncp);
129
130 const Real strikes[] = { 50, 75, 100, 125, 150, 200 };
131 for (Real strike : strikes) {
132 const Option::Type optionType = (strike > fwd) ? Option::Call : Option::Put;
133
134 const Real expected = BlackCalculator(
135 optionType, strike, fwd,
136 std::sqrt(x: volTS->blackVariance(t: maturity, strike)),
137 rTS->discount(t: maturity)).value();
138
139 const CLVModelPayoff clvModelPayoff(optionType, strike, g);
140
141 const ext::function<Real(Real)> f = [&](Real xi) -> Real {
142 return clvModelPayoff(xi) * boost::math::pdf(dist, x: xi);
143 };
144
145 const Real calculated = GaussLobattoIntegral(1000, 1e-6)(
146 f, x.front(), x.back()) * rTS->discount(t: maturity);
147
148 const Real tol = 5e-3;
149 if (std::fabs(x: expected - calculated) > tol) {
150 BOOST_FAIL("failed to reproduce option SquaredCLVMOdel prices"
151 << "\n time: " << maturityDate
152 << "\n strike: " << strike
153 << "\n expected: " << expected
154 << "\n calculated: " << calculated);
155 }
156 }
157}
158
159void SquareRootCLVModelTest::testSquareRootCLVMappingFunction() {
160 BOOST_TEST_MESSAGE(
161 "Testing mapping function of the square-root kernel process...");
162
163 using namespace square_root_clv_model;
164
165 const Date todaysDate(16, Oct, 2016);
166 Settings::instance().evaluationDate() = todaysDate;
167 const Date maturityDate = todaysDate + Period(1, Years);
168
169 const DayCounter dc = Actual365Fixed();
170
171 const Real s0 = 100;
172 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
173
174 const Rate r = 0.05;
175 const Rate q = 0.02;
176
177 const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
178 const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
179
180 //SABR
181 const Real beta = 0.95;
182 const Real alpha= 0.2;
183 const Real rho = -0.9;
184 const Real gamma= 0.8;
185
186 const Handle<BlackVolTermStructure> sabrVol(
187 ext::make_shared<SABRVolTermStructure>(
188 args: alpha, args: beta, args: gamma, args: rho, args: s0, args: r, args: todaysDate, args: dc));
189
190 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
191 ext::make_shared<GeneralizedBlackScholesProcess>(
192 args: spot, args: qTS, args: rTS, args: sabrVol));
193
194 std::vector<Date> calibrationDates(1, todaysDate + Period(3, Months));
195 calibrationDates.reserve(n: Size(daysBetween(d1: todaysDate, d2: maturityDate)/7 + 1));
196 while (calibrationDates.back() < maturityDate)
197 calibrationDates.push_back(x: calibrationDates.back() + Period(1, Weeks));
198
199 // sqrt process
200 const Real kappa = 1.0;
201 const Real theta = 0.09;
202 const Volatility sigma = 0.2;
203 const Real x0 = 0.09;
204
205 const ext::shared_ptr<SquareRootProcess> sqrtProcess(
206 ext::make_shared<SquareRootProcess>(args: theta, args: kappa, args: sigma, args: x0));
207
208 const SquareRootCLVModel model(
209 bsProcess, sqrtProcess, calibrationDates, 14, 1-1e-10, 1e-10);
210
211 const ext::function<Real(Time, Real)> g = model.g();
212
213 const Real strikes[] = { 80, 100, 120 };
214 const Size offsets[] = { 92, 182, 183, 184, 185, 186, 365 };
215 for (unsigned long offset : offsets) {
216 const Date m = todaysDate + Period(offset, Days);
217 const Time t = dc.yearFraction(d1: todaysDate, d2: m);
218
219 const Real df = 4*theta*kappa/(sigma*sigma);
220 const Real ncp = 4*kappa*std::exp(x: -kappa*t)
221 / (sigma*sigma*(1-std::exp(x: -kappa*t)))*sqrtProcess->x0();
222
223 const chi_squared_type dist(df, ncp);
224
225 const Real fwd = s0*qTS->discount(d: m)/rTS->discount(d: m);
226
227 for (Real strike : strikes) {
228 const Option::Type optionType = (strike > fwd) ? Option::Call : Option::Put;
229
230 const Real expected = BlackCalculator(
231 optionType, strike, fwd,
232 std::sqrt(x: sabrVol->blackVariance(d: m, strike)),
233 rTS->discount(d: m)).value();
234
235 const CLVModelPayoff clvModelPayoff(optionType, strike, [&](Real x) { return g(t, x); });
236
237 const ext::function<Real(Real)> f = [&](Real xi) -> Real {
238 return clvModelPayoff(xi) * boost::math::pdf(dist, x: xi);
239 };
240
241 const Array x = model.collocationPointsX(d: m);
242 const Real calculated = GaussLobattoIntegral(1000, 1e-3)(
243 f, x.front(), x.back()) * rTS->discount(d: m);
244
245 const Real tol = 0.075;
246
247 if (std::fabs(x: expected) > 0.01
248 && std::fabs(x: (calculated - expected)/calculated) > tol) {
249 BOOST_FAIL("failed to reproduce option SquaredCLVMOdel prices"
250 << "\n time: " << m
251 << "\n strike: " << strike
252 << "\n expected: " << expected
253 << "\n calculated: " << calculated);
254 }
255 }
256 }
257}
258
259namespace square_root_clv_model {
260 class SquareRootCLVCalibrationFunction : public CostFunction {
261 public:
262 SquareRootCLVCalibrationFunction(Array strikes,
263 const std::vector<Date>& resetDates,
264 const std::vector<Date>& maturityDates,
265 ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess,
266 Array refVols,
267 Size nScenarios = 10000)
268 : strikes_(std::move(strikes)), resetDates_(resetDates), maturityDates_(maturityDates),
269 bsProcess_(std::move(bsProcess)), refVols_(std::move(refVols)), nScenarios_(nScenarios) {
270 std::set<Date> c(resetDates.begin(), resetDates.end());
271 c.insert(first: maturityDates.begin(), last: maturityDates.end());
272 calibrationDates_.insert(
273 position: calibrationDates_.begin(), first: c.begin(), last: c.end());
274 }
275
276 Real value(const Array& params) const override {
277 const Array diff = values(params);
278
279 Real retVal = 0.0;
280 for (Real i : diff)
281 retVal += i * i;
282
283 return retVal;
284 }
285
286 Array values(const Array& params) const override {
287 const Real theta = params[0];
288 const Real kappa = params[1];
289 const Real sigma = params[2];
290 const Real x0 = params[3];
291
292 const ext::shared_ptr<SimpleQuote> vol(
293 ext::make_shared<SimpleQuote>(args: 0.1));
294
295 const Handle<YieldTermStructure> rTS(bsProcess_->riskFreeRate());
296 const Handle<YieldTermStructure> qTS(bsProcess_->dividendYield());
297 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(
298 args: bsProcess_->x0()));
299
300 const ext::shared_ptr<PricingEngine> fwdEngine(
301 ext::make_shared<ForwardVanillaEngine<AnalyticEuropeanEngine> >(
302 args: ext::make_shared<GeneralizedBlackScholesProcess>(
303 args: spot, args: qTS, args: rTS,
304 args: Handle<BlackVolTermStructure>(
305 flatVol(today: rTS->referenceDate(), volatility: vol,
306 dc: rTS->dayCounter())))));
307
308 const ext::shared_ptr<SquareRootProcess> sqrtProcess(
309 ext::make_shared<SquareRootProcess>(args: theta, args: kappa, args: sigma, args: x0));
310
311 const SquareRootCLVModel clvSqrtModel(
312 bsProcess_, sqrtProcess, calibrationDates_,
313 14, 1-1e-14, 1e-14);
314
315 const ext::function<Real(Time, Real)> gSqrt = clvSqrtModel.g();
316
317 Array retVal(resetDates_.size()*strikes_.size());
318
319 for (Size i=0, n=resetDates_.size(); i < n; ++i) {
320 const Date resetDate = resetDates_[i];
321 const Date maturityDate = maturityDates_[i];
322
323 const Time t0 = bsProcess_->time(resetDate);
324 const Time t1 = bsProcess_->time(maturityDate);
325
326 const Real df = 4*theta*kappa/(sigma*sigma);
327 const Real ncp = 4*kappa*std::exp(x: -kappa*t0)
328 / (sigma*sigma*(1-std::exp(x: -kappa*t0)))*x0;
329
330 typedef boost::math::non_central_chi_squared_distribution<Real>
331 chi_squared_type;
332
333 const chi_squared_type dist(df, ncp);
334
335 const Real ncp1 = 4*kappa*std::exp(x: -kappa*(t1-t0))
336 / (sigma*sigma*(1-std::exp(x: -kappa*(t1-t0))));
337
338 const LowDiscrepancy::ursg_type ursg = LowDiscrepancy::ursg_type(2, 1235UL);
339
340 std::vector<GeneralStatistics> stats(strikes_.size());
341
342 for (Size j=0; j < nScenarios_; ++j) {
343 const std::vector<Real>& path = ursg.nextSequence().value;
344
345 const Real x1 = boost::math::quantile(dist, p: path[0]);
346 const Real u1 =
347 sigma*sigma*(1-std::exp(x: -kappa*t0))/(4*kappa)*x1;
348
349 const Real x2 = boost::math::quantile(
350 dist: chi_squared_type(df, ncp1*u1), p: path[1]);
351 const Real u2 =
352 sigma*sigma*(1-std::exp(x: -kappa*(t1-t0)))/(4*kappa)*x2;
353 const Real X2 =
354 u2*4*kappa/(sigma*sigma*(1-std::exp(x: -kappa*t1)));
355
356 const Real s1 = gSqrt(t0, x1);
357 const Real s2 = gSqrt(t1, X2);
358
359 for (Size k=0; k < strikes_.size(); ++k) {
360 const Real strike = strikes_[k];
361
362 const Real payoff = (strike < 1.0)
363 ? Real(s1 * std::max(a: 0.0, b: strike - s2/s1))
364 : Real(s1 * std::max(a: 0.0, b: s2/s1 - strike));
365
366 stats[k].add(value: payoff);
367 }
368 }
369
370 const ext::shared_ptr<Exercise> exercise(
371 ext::make_shared<EuropeanExercise>(args: maturityDate));
372
373 const DiscountFactor dF(
374 bsProcess_->riskFreeRate()->discount(d: maturityDate));
375
376 for (Size k=0; k < strikes_.size(); ++k) {
377 const Real strike = strikes_[k];
378 const Real npv = stats[k].mean() * dF;
379
380 const ext::shared_ptr<StrikedTypePayoff> payoff(
381 ext::make_shared<PlainVanillaPayoff>(
382 args: (strike < 1.0) ? Option::Put : Option::Call, args: strike));
383
384 const ext::shared_ptr<ForwardVanillaOption> fwdOption(
385 ext::make_shared<ForwardVanillaOption>(
386 args: strike, args: resetDate, args: payoff, args: exercise));
387
388 const Volatility implVol =
389 QuantLib::detail::ImpliedVolatilityHelper::calculate(
390 instrument: *fwdOption, engine: *fwdEngine, volQuote&: *vol, targetValue: npv, accuracy: 1e-8, maxEvaluations: 200, minVol: 1e-4, maxVol: 2.0);
391
392 const Size idx = k + i*strikes_.size();
393 retVal[idx] = implVol - refVols_[idx];
394 }
395 }
396
397 return retVal;
398 }
399
400
401 private:
402 const Array strikes_;
403 const std::vector<Date> resetDates_, maturityDates_;
404 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess_;
405 const Array refVols_;
406 const Size nScenarios_;
407
408 std::vector<Date> calibrationDates_;
409 };
410
411 class NonZeroConstraint : public Constraint {
412 private:
413 class Impl : public Constraint::Impl {
414 public:
415 bool test(const Array& params) const override {
416 const Real theta = params[0];
417 const Real kappa = params[1];
418 const Real sigma = params[2];
419 const Real x0 = params[3];
420
421 return (sigma >= 0.001 && kappa > 1e-6 && theta > 0.001
422 && x0 > 1e-4);
423 }
424
425 Array upperBound(const Array& params) const override {
426 const Real upper[] = { 1.0, 1.0, 1.0, 2.0 };
427
428 return Array(upper, upper + 4);
429 }
430
431 Array lowerBound(const Array& params) const override {
432 const Real lower[] = { 0.001, 0.001, 0.001, 1e-4 };
433
434 return Array(lower, lower + 4);
435 }
436 };
437
438 public:
439 NonZeroConstraint()
440 : Constraint(ext::make_shared<NonZeroConstraint::Impl>()) {}
441 };
442}
443
444void SquareRootCLVModelTest::testForwardSkew() {
445 BOOST_TEST_MESSAGE(
446 "Testing forward skew dynamics with square-root kernel process...");
447
448 using namespace square_root_clv_model;
449
450 const Date todaysDate(16, Oct, 2016);
451 Settings::instance().evaluationDate() = todaysDate;
452 const Date endDate = todaysDate + Period(4, Years);
453
454 const DayCounter dc = Actual365Fixed();
455
456 // Heston model is used to generate an arbitrage free volatility surface
457 const Real s0 = 100;
458 const Real r = 0.1;
459 const Real q = 0.05;
460 const Real v0 = 0.09;
461 const Real kappa = 1.0;
462 const Real theta = 0.09;
463 const Real sigma = 0.3;
464 const Real rho = -0.75;
465
466 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: s0));
467 const Handle<YieldTermStructure> rTS(flatRate(forward: r, dc));
468 const Handle<YieldTermStructure> qTS(flatRate(forward: q, dc));
469
470 const ext::shared_ptr<HestonModel> hestonModel(
471 ext::make_shared<HestonModel>(
472 args: ext::make_shared<HestonProcess>(
473 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: sigma, args: rho)));
474
475 const Handle<BlackVolTermStructure> blackVol(
476 ext::make_shared<HestonBlackVolSurface>(
477 args: Handle<HestonModel>(hestonModel)));
478
479 const Handle<LocalVolTermStructure> localVol(
480 ext::make_shared<NoExceptLocalVolSurface>(
481 args: blackVol, args: rTS, args: qTS, args: spot, args: std::sqrt(x: theta)));
482
483 const Real sTheta = 0.389302;
484 const Real sKappa = 0.1101849;
485 const Real sSigma = 0.275368;
486 const Real sX0 = 0.466809;
487
488 const ext::shared_ptr<SquareRootProcess> sqrtProcess(
489 ext::make_shared<SquareRootProcess>(
490 args: sTheta, args: sKappa, args: sSigma, args: sX0));
491
492 const ext::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
493 ext::make_shared<GeneralizedBlackScholesProcess>(
494 args: spot, args: qTS, args: rTS, args: blackVol));
495
496 std::vector<Date> calibrationDates(1, todaysDate + Period(6, Months));
497 while (calibrationDates.back() < endDate)
498 calibrationDates.push_back(x: calibrationDates.back() + Period(3, Months));
499
500 std::set<Date> clvCalibrationDates(
501 calibrationDates.begin(), calibrationDates.end());
502
503 Date tmpDate = todaysDate + Period(1, Days);
504 while (tmpDate < todaysDate + Period(1, Years)) {
505 clvCalibrationDates.insert(x: tmpDate);
506 tmpDate += Period(1, Weeks);
507 }
508
509 const SquareRootCLVModel clvSqrtModel(
510 bsProcess,
511 sqrtProcess,
512 std::vector<Date>(
513 clvCalibrationDates.begin(), clvCalibrationDates.end()),
514 14, 1-1e-14, 1e-14);
515
516 const ext::function<Real(Time, Real)> gSqrt = clvSqrtModel.g();
517
518 const ext::shared_ptr<SimpleQuote> vol(
519 ext::make_shared<SimpleQuote>(args: 0.1));
520
521 const ext::shared_ptr<PricingEngine> fwdEngine(
522 ext::make_shared<ForwardVanillaEngine<AnalyticEuropeanEngine> >(
523 args: ext::make_shared<GeneralizedBlackScholesProcess>(
524 args: spot, args: qTS, args: rTS,
525 args: Handle<BlackVolTermStructure>(flatVol(today: todaysDate, volatility: vol, dc)))));
526
527
528 // forward skew of the Heston-SLV model
529 std::vector<Time> mandatoryTimes;
530 mandatoryTimes.reserve(n: calibrationDates.size());
531 for (auto& calibrationDate : calibrationDates)
532 mandatoryTimes.push_back(x: dc.yearFraction(d1: todaysDate, d2: calibrationDate));
533
534 const Size tSteps = 200;
535 const TimeGrid grid(mandatoryTimes.begin(), mandatoryTimes.end(), tSteps);
536
537 std::vector<Date> resetDates, maturityDates;
538 std::vector<Size> resetIndices, maturityIndices;
539 for (Size i=0, n = calibrationDates.size()-2; i < n; ++i) {
540 resetDates.push_back(x: calibrationDates[i]);
541 maturityDates.push_back(x: calibrationDates[i+2]);
542
543 const Time resetTime = mandatoryTimes[i];
544 const Time maturityTime = mandatoryTimes[i+2];
545
546 resetIndices.push_back(x: grid.closestIndex(t: resetTime)-1);
547 maturityIndices.push_back(x: grid.closestIndex(t: maturityTime)-1);
548 }
549
550 const Real strikes[] = {
551 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2,
552 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0
553 };
554
555 const Size nScenarios = 20000;
556 Array refVols(resetIndices.size()*LENGTH(strikes));
557
558 // finite difference calibration of Heston SLV model
559
560 // define Heston Stochastic Local Volatility model
561 const Real eta = 0.25;
562 const Real corr = -0.0;
563
564 const ext::shared_ptr<HestonProcess> hestonProcess4slv(
565 ext::make_shared<HestonProcess>(
566 args: rTS, args: qTS, args: spot, args: v0, args: kappa, args: theta, args: eta*sigma, args: corr));
567
568 const Handle<HestonModel> hestonModel4slv(
569 ext::make_shared<HestonModel>(args: hestonProcess4slv));
570
571 const HestonSLVFokkerPlanckFdmParams logParams = {
572 .xGrid: 301, .vGrid: 601, .tMaxStepsPerYear: 1000, .tMinStepsPerYear: 30, .tStepNumberDecay: 2.0, .nRannacherTimeSteps: 0, .predictionCorretionSteps: 2,
573 .x0Density: 0.1, .localVolEpsProb: 1e-4, .maxIntegrationIterations: 10000,
574 .vLowerEps: 1e-5, .vUpperEps: 1e-5, .vMin: 0.0000025, .v0Density: 1.0, .vLowerBoundDensity: 0.1, .vUpperBoundDensity: 0.9, .leverageFctPropEps: 1e-5,
575 .greensAlgorithm: FdmHestonGreensFct::Gaussian,
576 .trafoType: FdmSquareRootFwdOp::Log,
577 .schemeDesc: FdmSchemeDesc::ModifiedCraigSneyd()
578 };
579
580 const ext::shared_ptr<LocalVolTermStructure> leverageFctFDM =
581 HestonSLVFDMModel(localVol, hestonModel4slv, endDate, logParams).
582 leverageFunction();
583
584 // calibrating to forward volatility dynamics
585
586 const ext::shared_ptr<HestonSLVProcess> fdmSlvProcess(
587 ext::make_shared<HestonSLVProcess>(
588 args: hestonProcess4slv, args: leverageFctFDM));
589
590 std::vector<std::vector<GeneralStatistics> > slvStats(
591 calibrationDates.size()-2,
592 std::vector<GeneralStatistics>(LENGTH(strikes)));
593
594 typedef SobolBrownianBridgeRsg rsg_type;
595 typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
596
597 const Size factors = fdmSlvProcess->factors();
598
599 const ext::shared_ptr<MultiPathGenerator<rsg_type> > pathGen(
600 ext::make_shared<MultiPathGenerator<rsg_type> >(
601 args: fdmSlvProcess, args: grid, args: rsg_type(factors, grid.size()-1), args: false));
602
603 for (Size k=0; k < nScenarios; ++k) {
604 const sample_type& path = pathGen->next();
605
606 for (Size i=0, n=resetIndices.size(); i < n; ++i) {
607 const Real S_t1 = path.value[0][resetIndices[i]];
608 const Real S_T1 = path.value[0][maturityIndices[i]];
609
610 for (Size j=0; j < LENGTH(strikes); ++j) {
611 const Real strike = strikes[j];
612 slvStats[i][j].add(value: (strike < 1.0)
613 ? Real(S_t1 * std::max(a: 0.0, b: strike - S_T1/S_t1))
614 : Real(S_t1 * std::max(a: 0.0, b: S_T1/S_t1 - strike)));
615 }
616
617 }
618 }
619
620 for (Size i=0, n=resetIndices.size(); i < n; ++i) {
621 const Date resetDate = calibrationDates[i];
622 const Date maturityDate(calibrationDates[i+2]);
623 const DiscountFactor df = rTS->discount(d: maturityDate);
624
625 const ext::shared_ptr<Exercise> exercise(
626 ext::make_shared<EuropeanExercise>(args: maturityDate));
627
628 for (Size j=0; j < LENGTH(strikes); ++j) {
629 const Real strike = strikes[j];
630 const Real npv = slvStats[i][j].mean()*df;
631
632 const ext::shared_ptr<StrikedTypePayoff> payoff(
633 ext::make_shared<PlainVanillaPayoff>(
634 args: (strike < 1.0) ? Option::Put : Option::Call, args: strike));
635
636 const ext::shared_ptr<ForwardVanillaOption> fwdOption(
637 ext::make_shared<ForwardVanillaOption>(
638 args: strike, args: resetDate, args: payoff, args: exercise));
639
640 const Volatility implVol =
641 QuantLib::detail::ImpliedVolatilityHelper::calculate(
642 instrument: *fwdOption, engine: *fwdEngine, volQuote&: *vol, targetValue: npv, accuracy: 1e-8, maxEvaluations: 200, minVol: 1e-4, maxVol: 2.0);
643
644 const Size idx = j + i*LENGTH(strikes);
645 refVols[idx] = implVol;
646 }
647 }
648
649 SquareRootCLVCalibrationFunction costFunction(
650 Array(strikes, strikes+LENGTH(strikes)),
651 resetDates,
652 maturityDates,
653 bsProcess,
654 refVols,
655 nScenarios);
656
657 NonZeroConstraint nonZeroConstraint;
658
659 CompositeConstraint constraint(
660 nonZeroConstraint,
661 HestonModel::FellerConstraint());
662
663 Array params(4);
664 params[0] = sTheta; params[1] = sKappa;
665 params[2] = sSigma; params[3] = sX0;
666
667
668 // Optimization would take too long
669 //
670 // Problem prob(costFunction, nonZeroConstraint, params);
671 //
672 // Simplex simplex(0.05);
673 // simplex.minimize(prob, EndCriteria(400, 40, 1.0e-8, 1.0e-8, 1.0e-8));
674
675 const Real tol = 0.5;
676 const Real costValue = costFunction.value(params);
677
678 if (costValue > tol) {
679 BOOST_FAIL("failed to reproduce small cost function value"
680 << "\n value: " << costValue
681 << "\n tolerance: " << tol);
682 }
683
684 const Date maturityDate = todaysDate + Period(1, Years);
685 const Time maturityTime = bsProcess->time(maturityDate);
686
687 const ext::shared_ptr<Exercise> europeanExercise(
688 ext::make_shared<EuropeanExercise>(args: maturityDate));
689
690 VanillaOption vanillaATMOption(
691 ext::make_shared<PlainVanillaPayoff>(args: Option::Call,
692 args: s0*qTS->discount(d: maturityDate)/rTS->discount(d: maturityDate)),
693 europeanExercise);
694
695 vanillaATMOption.setPricingEngine(
696 ext::make_shared<AnalyticHestonEngine>(args: hestonModel));
697
698 const Volatility atmVol = vanillaATMOption.impliedVolatility(
699 price: vanillaATMOption.NPV(),
700 process: ext::make_shared<GeneralizedBlackScholesProcess>(args: spot, args: qTS, args: rTS,
701 args: Handle<BlackVolTermStructure>(flatVol(volatility: std::sqrt(x: theta), dc))));
702
703 const ext::shared_ptr<PricingEngine> analyticEngine(
704 ext::make_shared<AnalyticDoubleBarrierBinaryEngine>(
705 args: ext::make_shared<GeneralizedBlackScholesProcess>(
706 args: spot, args: qTS, args: rTS,
707 args: Handle<BlackVolTermStructure>(flatVol(volatility: atmVol, dc)))));
708
709 const ext::shared_ptr<PricingEngine> fdSLVEngine(
710 ext::make_shared<FdHestonDoubleBarrierEngine>(
711 args: hestonModel4slv.currentLink(),
712 args: 51, args: 201, args: 51, args: 1,
713 args: FdmSchemeDesc::Hundsdorfer(), args: leverageFctFDM));
714
715 const Size n = 16;
716 Array barrier_lo(n), barrier_hi(n), bsNPV(n), slvNPV(n);
717
718 const ext::shared_ptr<CashOrNothingPayoff> payoff =
719 ext::make_shared<CashOrNothingPayoff>(args: Option::Call, args: 0.0, args: 1.0);
720
721 for (Size i=0; i < n; ++i) {
722 const Real dist = 20.0+5.0*i;
723
724 barrier_lo[i] = std::max(a: s0 - dist, b: 1e-2);
725 barrier_hi[i] = s0 + dist;
726 DoubleBarrierOption doubleBarrier(
727 DoubleBarrier::KnockOut, barrier_lo[i], barrier_hi[i], 0.0,
728 payoff,
729 europeanExercise);
730
731 doubleBarrier.setPricingEngine(analyticEngine);
732 bsNPV[i] = doubleBarrier.NPV();
733
734 doubleBarrier.setPricingEngine(fdSLVEngine);
735 slvNPV[i] = doubleBarrier.NPV();
736 }
737
738
739 const TimeGrid bGrid(maturityTime, tSteps);
740
741 const PseudoRandom::ursg_type ursg = PseudoRandom::ursg_type(tSteps, 1235UL);
742
743 std::vector<GeneralStatistics> stats(n);
744
745 const Real df = 4*sTheta*sKappa/(sSigma*sSigma);
746
747 for (Size i=0; i < nScenarios; ++i) {
748 std::vector<bool> touch(n, false);
749
750 const std::vector<Real>& path = ursg.nextSequence().value;
751
752 Real x = sX0;
753
754 for (Size j=0; j < tSteps; ++j) {
755 const Time t0 = bGrid.at(i: j);
756 const Time t1 = bGrid.at(i: j+1);
757
758 const Real ncp = 4*sKappa*std::exp(x: -sKappa*(t1-t0))
759 / (sSigma*sSigma*(1-std::exp(x: -sKappa*(t1-t0))))*x;
760
761 const boost::math::non_central_chi_squared_distribution<Real>
762 dist(df, ncp);
763
764 const Real u = boost::math::quantile(dist, p: path[j]);
765
766 x = sSigma*sSigma*(1-std::exp(x: -sKappa*(t1-t0)))/(4*sKappa) * u;
767
768 const Real X = x*4*sKappa/(sSigma*sSigma*(1-std::exp(x: -sKappa*t1)));
769
770 const Real s = gSqrt(t1, X);
771
772 if (t1 > 0.05) {
773 for (Size u=0; u < n; ++u) {
774 if (s <= barrier_lo[u] || s >= barrier_hi[u]) {
775 touch[u] = true;
776 }
777 }
778 }
779 }
780 for (Size u=0; u < n; ++u) {
781 if (touch[u]) {
782 stats[u].add(value: 0.0);
783 }
784 else {
785 stats[u].add(value: rTS->discount(d: maturityDate));
786 }
787 }
788 }
789
790
791 for (Size u=0; u < n; ++u) {
792 const Real calculated = stats[u].mean();
793 const Real error = stats[u].errorEstimate();
794 const Real expected = slvNPV[u];
795
796 const Real tol = 2.35*error;
797
798 if (std::fabs(x: calculated-expected) > tol) {
799 BOOST_FAIL("failed to reproduce CLV double no touch barrier price"
800 << "\n CLV value: " << calculated
801 << "\n error : " << error
802 << "\n SLV value: " << expected);
803 }
804 }
805}
806
807
808test_suite* SquareRootCLVModelTest::experimental() {
809 auto* suite = BOOST_TEST_SUITE("SquareRootCLVModel tests");
810
811 suite->add(QUANTLIB_TEST_CASE(
812 &SquareRootCLVModelTest::testSquareRootCLVVanillaPricing));
813
814 suite->add(QUANTLIB_TEST_CASE(
815 &SquareRootCLVModelTest::testSquareRootCLVMappingFunction));
816
817// this test takes very long
818// suite->add(QUANTLIB_TEST_CASE(
819// &SquareRootCLVModelTest::testForwardSkew));
820
821 return suite;
822}
823

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