1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Andreas Gaida
5 Copyright (C) 2008 Ralph Schreyer
6 Copyright (C) 2008, 2009, 2010, 2015 Klaus Spanderen
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22#include "fdmlinearop.hpp"
23#include "utilities.hpp"
24
25#include <ql/quotes/simplequote.hpp>
26#include <ql/time/daycounters/actual360.hpp>
27#include <ql/time/daycounters/actual365fixed.hpp>
28#include <ql/processes/hestonprocess.hpp>
29#include <ql/processes/hullwhiteprocess.hpp>
30#include <ql/processes/blackscholesprocess.hpp>
31#include <ql/processes/hybridhestonhullwhiteprocess.hpp>
32#include <ql/math/interpolations/bilinearinterpolation.hpp>
33#include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
34#include <ql/math/interpolations/cubicinterpolation.hpp>
35#include <ql/math/integrals/discreteintegrals.hpp>
36#include <ql/math/randomnumbers/rngtraits.hpp>
37#include <ql/models/equity/hestonmodel.hpp>
38#include <ql/termstructures/yield/zerocurve.hpp>
39#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
40#include <ql/pricingengines/vanilla/mchestonhullwhiteengine.hpp>
41#include <ql/methods/finitedifferences/finitedifferencemodel.hpp>
42#include <ql/math/matrixutilities/gmres.hpp>
43#include <ql/math/matrixutilities/bicgstab.hpp>
44#include <ql/methods/finitedifferences/schemes/douglasscheme.hpp>
45#include <ql/methods/finitedifferences/schemes/hundsdorferscheme.hpp>
46#include <ql/methods/finitedifferences/schemes/craigsneydscheme.hpp>
47#include <ql/methods/finitedifferences/meshers/uniformgridmesher.hpp>
48#include <ql/methods/finitedifferences/meshers/uniform1dmesher.hpp>
49#include <ql/methods/finitedifferences/meshers/concentrating1dmesher.hpp>
50#include <ql/methods/finitedifferences/meshers/fdmblackscholesmesher.hpp>
51#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
52#include <ql/methods/finitedifferences/operators/fdmblackscholesop.hpp>
53#include <ql/methods/finitedifferences/utilities/fdmmesherintegral.hpp>
54#include <ql/methods/finitedifferences/utilities/fdminnervaluecalculator.hpp>
55#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
56#include <ql/methods/finitedifferences/operators/fdmlinearop.hpp>
57#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
58#include <ql/methods/finitedifferences/operators/fdmlinearopcomposite.hpp>
59#include <ql/methods/finitedifferences/operators/fdmhestonhullwhiteop.hpp>
60#include <ql/methods/finitedifferences/meshers/fdmhestonvariancemesher.hpp>
61#include <ql/methods/finitedifferences/operators/fdmhestonop.hpp>
62#include <ql/methods/finitedifferences/solvers/fdmhestonsolver.hpp>
63#include <ql/methods/finitedifferences/meshers/fdmmeshercomposite.hpp>
64#include <ql/methods/finitedifferences/solvers/fdmndimsolver.hpp>
65#include <ql/methods/finitedifferences/solvers/fdm3dimsolver.hpp>
66#include <ql/methods/finitedifferences/stepconditions/fdmamericanstepcondition.hpp>
67#include <ql/methods/finitedifferences/stepconditions/fdmstepconditioncomposite.hpp>
68#include <ql/methods/finitedifferences/utilities/fdmdividendhandler.hpp>
69#include <ql/methods/finitedifferences/operators/firstderivativeop.hpp>
70#include <ql/methods/finitedifferences/operators/secondderivativeop.hpp>
71#include <ql/methods/finitedifferences/operators/secondordermixedderivativeop.hpp>
72#include <ql/math/matrixutilities/sparseilupreconditioner.hpp>
73#include <ql/functional.hpp>
74
75#include <boost/numeric/ublas/vector.hpp>
76#include <boost/numeric/ublas/operation.hpp>
77
78#include <numeric>
79#include <utility>
80
81using namespace QuantLib;
82using namespace boost::unit_test_framework;
83
84namespace {
85
86 class FdmHestonExpressCondition : public StepCondition<Array> {
87 public:
88 FdmHestonExpressCondition(std::vector<Real> redemptions,
89 std::vector<Real> triggerLevels,
90 std::vector<Time> exerciseTimes,
91 ext::shared_ptr<FdmMesher> mesher)
92 : redemptions_(std::move(redemptions)), triggerLevels_(std::move(triggerLevels)),
93 exerciseTimes_(std::move(exerciseTimes)), mesher_(std::move(mesher)) {}
94
95 void applyTo(Array& a, Time t) const override {
96 auto iter = std::find(first: exerciseTimes_.begin(), last: exerciseTimes_.end(), val: t);
97
98 if (iter != exerciseTimes_.end()) {
99 Size index = std::distance(first: exerciseTimes_.begin(), last: iter);
100
101 for (const auto& iter : *mesher_->layout()) {
102 const Real s = std::exp(x: mesher_->location(iter, direction: 0));
103
104 if (s > triggerLevels_[index]) {
105 a[iter.index()] = redemptions_[index];
106 }
107 }
108 }
109 }
110
111 private:
112 const std::vector<Real> redemptions_;
113 const std::vector<Real> triggerLevels_;
114 const std::vector<Time> exerciseTimes_;
115 const ext::shared_ptr<FdmMesher> mesher_;
116 };
117
118 class ExpressPayoff : public Payoff {
119 public:
120 std::string name() const override { return "ExpressPayoff"; }
121 std::string description() const override { return "ExpressPayoff"; }
122
123 Real operator()(Real s) const override {
124 return ((s >= 100.0) ? 108.0 : 100.0)
125 - ((s <= 75.0) ? Real(100.0 - s) : 0.0);
126 }
127 };
128
129 template <class T, class U, class V>
130 struct multiplies {
131 V operator()(T t, U u) { return t*u;}
132 };
133
134}
135
136void FdmLinearOpTest::testFdmLinearOpLayout() {
137
138 BOOST_TEST_MESSAGE("Testing indexing of a linear operator...");
139
140 const std::vector<Size> dim = {5,7,8};
141
142 FdmLinearOpLayout layout = FdmLinearOpLayout(dim);
143
144 Size calculatedDim = layout.dim().size();
145 Size expectedDim = dim.size();
146 if (calculatedDim != expectedDim) {
147 BOOST_ERROR("index.dimensions() should be " << expectedDim
148 << ", but is " << calculatedDim);
149 }
150
151 Size calculatedSize = layout.size();
152 Size expectedSize = std::accumulate(first: dim.begin(), last: dim.end(), init: 1, binary_op: std::multiplies<>());
153
154 if (calculatedSize != expectedSize) {
155 BOOST_FAIL("index.size() should be "
156 << expectedSize << ", but is " << calculatedSize);
157 }
158
159 for (Size k=0; k < dim[0]; ++k) {
160 for (Size l=0; l < dim[1]; ++l) {
161 for (Size m=0; m < dim[2]; ++m) {
162 std::vector<Size> tmp(3);
163 tmp[0] = k; tmp[1] = l; tmp[2] = m;
164
165 Size calculatedIndex = layout.index(coordinates: tmp);
166 Size expectedIndex = k + l*dim[0] + m*dim[0]*dim[1];
167
168 if (expectedIndex != layout.index(coordinates: tmp)) {
169 BOOST_FAIL("index.size() should be " << expectedIndex
170 <<", but is " << calculatedIndex);
171 }
172 }
173 }
174 }
175
176 FdmLinearOpIterator iter = layout.begin();
177
178 for (Size m=0; m < dim[2]; ++m) {
179 for (Size l=0; l < dim[1]; ++l) {
180 for (Size k=0; k < dim[0]; ++k, ++iter) {
181 for (Size n=1; n < 4; ++n) {
182 Size nn = layout.neighbourhood(iterator: iter, i: 1, offset: n);
183 Size calculatedIndex = k + m*dim[0]*dim[1]
184 + ((l < dim[1]-n)? l+n
185 : dim[1]-1-(l+n-(dim[1]-1)))*dim[0];
186
187 if (nn != calculatedIndex) {
188 BOOST_FAIL("next neighbourhood index is " << nn
189 << " but should be " << calculatedIndex);
190 }
191 }
192
193 for (Size n=1; n < 7; ++n) {
194 Size nn = layout.neighbourhood(iterator: iter, i: 2, offset: -Integer(n));
195 Size calculatedIndex = k + l*dim[0]
196 + ((m < n) ? n-m : m-n)*dim[0]*dim[1];
197 if (nn != calculatedIndex) {
198 BOOST_FAIL("next neighbourhood index is " << nn
199 << " but should be " << calculatedIndex);
200 }
201 }
202 }
203 }
204 }
205}
206
207void FdmLinearOpTest::testUniformGridMesher() {
208
209 BOOST_TEST_MESSAGE("Testing uniform grid mesher...");
210
211 const std::vector<Size> dim = {5,7,8};
212
213 ext::shared_ptr<FdmLinearOpLayout> layout(new FdmLinearOpLayout(dim));
214 std::vector<std::pair<Real, Real> > boundaries = {{-5, 10}, {5, 100}, {10, 20}};
215
216 UniformGridMesher mesher(layout, boundaries);
217
218 const Real dx1 = 15.0/(dim[0]-1);
219 const Real dx2 = 95.0/(dim[1]-1);
220 const Real dx3 = 10.0/(dim[2]-1);
221
222 constexpr double tol = 100*QL_EPSILON;
223 if ( std::fabs(x: dx1-mesher.dminus(layout->begin(),direction: 0)) > tol
224 || std::fabs(x: dx1-mesher.dplus(layout->begin(),direction: 0)) > tol
225 || std::fabs(x: dx2-mesher.dminus(layout->begin(),direction: 1)) > tol
226 || std::fabs(x: dx2-mesher.dplus(layout->begin(),direction: 1)) > tol
227 || std::fabs(x: dx3-mesher.dminus(layout->begin(),direction: 2)) > tol
228 || std::fabs(x: dx3-mesher.dplus(layout->begin(),direction: 2)) > tol ) {
229 BOOST_FAIL("inconsistent uniform mesher object");
230 }
231}
232
233void FdmLinearOpTest::testFirstDerivativesMapApply() {
234
235 BOOST_TEST_MESSAGE("Testing application of first-derivatives map...");
236
237 const std::vector<Size> dim = {400, 100, 50};
238
239 ext::shared_ptr<FdmLinearOpLayout> index(new FdmLinearOpLayout(dim));
240
241 std::vector<std::pair<Real, Real> > boundaries = {{-5, 5}, {0, 10}, { 5, 15}};
242
243 ext::shared_ptr<FdmMesher> mesher(
244 new UniformGridMesher(index, boundaries));
245
246 FirstDerivativeOp map(2, mesher);
247
248 Array r(mesher->layout()->size());
249 for (const auto& iter : *index) {
250 r[iter.index()] = std::sin(x: mesher->location(iter, direction: 0))
251 + std::cos(x: mesher->location(iter, direction: 2));
252 }
253
254 Array t = map.apply(r);
255 const Real dz = (boundaries[2].second-boundaries[2].first)/(dim[2]-1);
256 for (const auto& iter : *index) {
257 const Size z = iter.coordinates()[2];
258
259 const Size z0 = (z > 0) ? z-1 : 1;
260 const Size z2 = (z < dim[2]-1) ? z+1 : dim[2]-2;
261 const Real lz0 = boundaries[2].first + z0*dz;
262 const Real lz2 = boundaries[2].first + z2*dz;
263
264 Real expected;
265 if (z == 0) {
266 expected = (std::cos(x: boundaries[2].first+dz)
267 - std::cos(x: boundaries[2].first))/dz;
268 }
269 else if (z == dim[2]-1) {
270 expected = (std::cos(x: boundaries[2].second)
271 - std::cos(x: boundaries[2].second-dz))/dz;
272 }
273 else {
274 expected = (std::cos(x: lz2)-std::cos(x: lz0))/(2*dz);
275 }
276
277 const Real calculated = t[iter.index()];
278 if (std::fabs(x: calculated - expected) > 1e-10) {
279 BOOST_FAIL("first derivative calculation failed."
280 << "\n calculated: " << calculated
281 << "\n expected: " << expected);
282 }
283 }
284
285
286}
287
288void FdmLinearOpTest::testSecondDerivativesMapApply() {
289
290 BOOST_TEST_MESSAGE("Testing application of second-derivatives map...");
291
292 const std::vector<Size> dim = {50, 50, 50};
293
294 ext::shared_ptr<FdmLinearOpLayout> index(new FdmLinearOpLayout(dim));
295
296 std::vector<std::pair<Real, Real> > boundaries = {{0, 0.5}, {0, 0.5}, {0, 0.5}};
297
298 ext::shared_ptr<FdmMesher> mesher(
299 new UniformGridMesher(index, boundaries));
300 Array r(mesher->layout()->size());
301 for (const auto& iter : *index) {
302 const Real x = mesher->location(iter, direction: 0);
303 const Real y = mesher->location(iter, direction: 1);
304 const Real z = mesher->location(iter, direction: 2);
305
306 r[iter.index()] = std::sin(x: x)*std::cos(x: y)*std::exp(x: z);
307 }
308
309 Array t = SecondDerivativeOp(0, mesher).apply(r);
310
311 const Real tol = 5e-2;
312 for (const auto& iter : *index) {
313 const Size i = iter.index();
314 const Real x = mesher->location(iter, direction: 0);
315 const Real y = mesher->location(iter, direction: 1);
316 const Real z = mesher->location(iter, direction: 2);
317
318 Real d = -std::sin(x: x)*std::cos(x: y)*std::exp(x: z);
319 if (iter.coordinates()[0] == 0 || iter.coordinates()[0] == dim[0]-1) {
320 d = 0;
321 }
322
323 if (std::fabs(x: d - t[i]) > tol) {
324 BOOST_FAIL("numerical derivative in dx^2 deviation is too big"
325 << "\n found at " << x << " " << y << " " << z);
326 }
327 }
328
329 t = SecondDerivativeOp(1, mesher).apply(r);
330 for (const auto& iter : *index) {
331 const Size i = iter.index();
332 const Real x = mesher->location(iter, direction: 0);
333 const Real y = mesher->location(iter, direction: 1);
334 const Real z = mesher->location(iter, direction: 2);
335
336 Real d = -std::sin(x: x)*std::cos(x: y)*std::exp(x: z);
337 if (iter.coordinates()[1] == 0 || iter.coordinates()[1] == dim[1]-1) {
338 d = 0;
339 }
340
341 if (std::fabs(x: d - t[i]) > tol) {
342 BOOST_FAIL("numerical derivative in dy^2 deviation is too big"
343 << "\n found at " << x << " " << y << " " << z);
344 }
345 }
346
347 t = SecondDerivativeOp(2, mesher).apply(r);
348 for (const auto& iter : *index) {
349 const Size i = iter.index();
350 const Real x = mesher->location(iter, direction: 0);
351 const Real y = mesher->location(iter, direction: 1);
352 const Real z = mesher->location(iter, direction: 2);
353
354 Real d = std::sin(x: x)*std::cos(x: y)*std::exp(x: z);
355 if (iter.coordinates()[2] == 0 || iter.coordinates()[2] == dim[2]-1) {
356 d = 0;
357 }
358
359 if (std::fabs(x: d - t[i]) > tol) {
360 BOOST_FAIL("numerical derivative in dz^2 deviation is too big"
361 << "\n found at " << x << " " << y << " " << z);
362 }
363 }
364
365
366}
367
368void FdmLinearOpTest::testDerivativeWeightsOnNonUniformGrids() {
369 BOOST_TEST_MESSAGE("Testing finite differences coefficients...");
370
371 const ext::shared_ptr<Fdm1dMesher> mesherX(
372 new Concentrating1dMesher(-2.0, 3.0, 50, std::make_pair(x: 0.5, y: 0.01)));
373 const ext::shared_ptr<Fdm1dMesher> mesherY(
374 new Concentrating1dMesher(0.5, 5.0, 25, std::make_pair(x: 0.5, y: 0.1)));
375 const ext::shared_ptr<Fdm1dMesher> mesherZ(
376 new Concentrating1dMesher(-1.0, 2.0, 31, std::make_pair(x: 1.5, y: 0.01)));
377
378 const ext::shared_ptr<FdmMesher> meshers(
379 new FdmMesherComposite(mesherX, mesherY, mesherZ));
380
381 const Real tol = 1e-13;
382 for (Size direction=0; direction < 3; ++direction) {
383
384 const SparseMatrix dfdx
385 = FirstDerivativeOp(direction, meshers).toMatrix();
386 const SparseMatrix d2fdx2
387 = SecondDerivativeOp(direction, meshers).toMatrix();
388
389 const Array gridPoints = meshers->locations(direction);
390
391 for (const auto& iter : *meshers->layout()) {
392
393 const Size c = iter.coordinates()[direction];
394 const Size index = iter.index();
395 const Size indexM1 = meshers->layout()->neighbourhood(iterator: iter,i: direction,offset: -1);
396 const Size indexP1 = meshers->layout()->neighbourhood(iterator: iter,i: direction,offset: +1);
397
398 // test only if not on the boundary
399 if (c == 0) {
400 Array twoPoints(2);
401 twoPoints[0] = 0.0;
402 twoPoints[1] = gridPoints.at(i: indexP1)-gridPoints.at(i: index);
403
404 const Array ndWeights1st =
405 NumericalDifferentiation({}, 1 , twoPoints).weights();
406
407 const Real beta1 = dfdx(index, index);
408 const Real gamma1 = dfdx(index, indexP1);
409 if ( std::fabs(x: (beta1 - ndWeights1st.at(i: 0))/beta1) > tol
410 || std::fabs(x: (gamma1 - ndWeights1st.at(i: 1))/gamma1) > tol) {
411 BOOST_FAIL("can not reproduce the weights of the "
412 "first order derivative operator "
413 "on the lower boundary"
414 << "\n expected beta: " << ndWeights1st.at(0)
415 << "\n calculated beta: " << beta1
416 << "\n difference beta: "
417 << beta1 - ndWeights1st.at(0)
418 << "\n expected gamma: " << ndWeights1st.at(1)
419 << "\n calculated gamma: " << gamma1
420 << "\n difference gamma: "
421 << gamma1 - ndWeights1st.at(1));
422 }
423
424 // free boundary condition by default
425 const Real beta2 = d2fdx2(index, index);
426 const Real gamma2 = d2fdx2(index, indexP1);
427
428 if ( std::fabs(x: beta2) > QL_EPSILON
429 || std::fabs(x: gamma2) > QL_EPSILON) {
430 BOOST_FAIL("can not reproduce the weights of the "
431 "second order derivative operator "
432 "on the lower boundary"
433 << "\n expected beta: " << 0.0
434 << "\n calculated beta: " << beta2
435 << "\n expected gamma: " << 0.0
436 << "\n calculated gamma: " << gamma2);
437 }
438 }
439 else if (c == meshers->layout()->dim()[direction]-1) {
440 Array twoPoints(2);
441 twoPoints[0] = gridPoints.at(i: indexM1)-gridPoints.at(i: index);
442 twoPoints[1] = 0.0;
443
444 const Array ndWeights1st =
445 NumericalDifferentiation({}, 1 , twoPoints).weights();
446
447 const Real alpha1 = dfdx(index, indexM1);
448 const Real beta1 = dfdx(index, index);
449 if ( std::fabs(x: (alpha1 - ndWeights1st.at(i: 0))/alpha1) > tol
450 || std::fabs(x: (beta1 - ndWeights1st.at(i: 1))/beta1) > tol) {
451 BOOST_FAIL("can not reproduce the weights of the "
452 "first order derivative operator "
453 "on the upper boundary"
454 << "\n expected alpha: " << ndWeights1st.at(0)
455 << "\n calculated alpha: " << alpha1
456 << "\n difference alpha: "
457 << alpha1 - ndWeights1st.at(0)
458 << "\n expected beta: " << ndWeights1st.at(1)
459 << "\n calculated beta: " << beta1
460 << "\n difference beta: "
461 << beta1 - ndWeights1st.at(1));
462 }
463
464 // free boundary condition by default
465 const Real alpha2 = d2fdx2(index, indexM1);
466 const Real beta2 = d2fdx2(index, index);
467
468 if ( std::fabs(x: alpha2) > QL_EPSILON
469 || std::fabs(x: beta2) > QL_EPSILON) {
470 BOOST_FAIL("can not reproduce the weights of the "
471 "second order derivative operator "
472 "on the upper boundary"
473 << "\n expected alpha: " << 0.0
474 << "\n calculated alpha: " << alpha2
475 << "\n expected beta: " << 0.0
476 << "\n calculated beta: " << beta2);
477 }
478 }
479 else {
480 Array threePoints(3);
481 threePoints[0] = gridPoints.at(i: indexM1)-gridPoints.at(i: index);
482 threePoints[1] = 0.0;
483 threePoints[2] = gridPoints.at(i: indexP1)-gridPoints.at(i: index);
484
485 const Array ndWeights1st =
486 NumericalDifferentiation({}, 1 , threePoints).weights();
487
488 const Real alpha1 = dfdx(index, indexM1);
489 const Real beta1 = dfdx(index, index);
490 const Real gamma1 = dfdx(index, indexP1);
491
492 if ( std::fabs(x: (alpha1 - ndWeights1st.at(i: 0))/alpha1) > tol
493 || std::fabs(x: (beta1 - ndWeights1st.at(i: 1))/beta1) > tol
494 || std::fabs(x: (gamma1 - ndWeights1st.at(i: 2))/gamma1) > tol) {
495 BOOST_FAIL("can not reproduce the weights of the "
496 "first order derivative operator"
497 << "\n expected alpha: " << ndWeights1st.at(0)
498 << "\n calculated alpha: " << alpha1
499 << "\n difference alpha: "
500 << alpha1 - ndWeights1st.at(0)
501 << "\n expected beta: " << ndWeights1st.at(1)
502 << "\n calculated beta: " << beta1
503 << "\n difference beta: "
504 << beta1 - ndWeights1st.at(1)
505 << "\n expected gamma: " << ndWeights1st.at(2)
506 << "\n calculated gamma: " << gamma1
507 << "\n difference gamma: "
508 << gamma1 - ndWeights1st.at(2));
509 }
510
511 const Array ndWeights2nd =
512 NumericalDifferentiation({}, 2 , threePoints).weights();
513
514 const Real alpha2 = d2fdx2(index, indexM1);
515 const Real beta2 = d2fdx2(index, index);
516 const Real gamma2 = d2fdx2(index, indexP1);
517 if ( std::fabs(x: (alpha2 - ndWeights2nd.at(i: 0))/alpha2) > tol
518 || std::fabs(x: (beta2 - ndWeights2nd.at(i: 1))/beta2) > tol
519 || std::fabs(x: (gamma2 - ndWeights2nd.at(i: 2))/gamma2) > tol) {
520 BOOST_FAIL("can not reproduce the weights of the "
521 "second order derivative operator"
522 << "\n expected alpha: " << ndWeights2nd.at(0)
523 << "\n calculated alpha: " << alpha2
524 << "\n difference alpha: "
525 << alpha2 - ndWeights2nd.at(0)
526 << "\n expected beta: " << ndWeights2nd.at(1)
527 << "\n calculated beta: " << beta2
528 << "\n difference beta: "
529 << beta2 - ndWeights2nd.at(1)
530 << "\n expected gamma: " << ndWeights2nd.at(2)
531 << "\n calculated gamma: " << gamma2
532 << "\n difference gamma: "
533 << gamma2 - ndWeights2nd.at(2));
534 }
535 }
536 }
537 }
538}
539
540void FdmLinearOpTest::testSecondOrderMixedDerivativesMapApply() {
541
542 BOOST_TEST_MESSAGE(
543 "Testing application of second-order mixed-derivatives map...");
544
545 const std::vector<Size> dim = {50, 50, 50};
546
547 ext::shared_ptr<FdmLinearOpLayout> index(new FdmLinearOpLayout(dim));
548
549 std::vector<std::pair<Real, Real> > boundaries = {{0, 0.5}, {0, 0.5}, {0, 0.5}};
550
551 ext::shared_ptr<FdmMesher> mesher(
552 new UniformGridMesher(index, boundaries));
553
554 Array r(mesher->layout()->size());
555
556 for (const auto& iter : *index) {
557 const Real x = mesher->location(iter, direction: 0);
558 const Real y = mesher->location(iter, direction: 1);
559 const Real z = mesher->location(iter, direction: 2);
560
561 r[iter.index()] = std::sin(x: x)*std::cos(x: y)*std::exp(x: z);
562 }
563
564 Array t = SecondOrderMixedDerivativeOp(0, 1, mesher).apply(r);
565 Array u = SecondOrderMixedDerivativeOp(1, 0, mesher).apply(r);
566
567 const Real tol = 5e-2;
568 for (const auto& iter : *index) {
569 const Size i = iter.index();
570 const Real x = mesher->location(iter, direction: 0);
571 const Real y = mesher->location(iter, direction: 1);
572 const Real z = mesher->location(iter, direction: 2);
573
574 const Real d = -std::cos(x: x)*std::sin(x: y)*std::exp(x: z);
575
576 if (std::fabs(x: d - t[i]) > tol) {
577 BOOST_FAIL("numerical derivative in dxdy deviation is too big"
578 << "\n found at " << x << " " << y << " " << z);
579 }
580
581 if (std::fabs(x: t[i]-u[i]) > 1e5*QL_EPSILON) {
582 BOOST_FAIL("numerical derivative in dxdy not equal to dydx"
583 << "\n found at " << x << " " << y << " " << z
584 << "\n value " << std::fabs(t[i]-u[i]));
585 }
586 }
587
588 t = SecondOrderMixedDerivativeOp(0, 2, mesher).apply(r);
589 u = SecondOrderMixedDerivativeOp(2, 0, mesher).apply(r);
590 for (const auto& iter : *index) {
591 const Size i = iter.index();
592 const Real x = mesher->location(iter, direction: 0);
593 const Real y = mesher->location(iter, direction: 1);
594 const Real z = mesher->location(iter, direction: 2);
595
596 const Real d = std::cos(x: x)*std::cos(x: y)*std::exp(x: z);
597
598 if (std::fabs(x: d - t[i]) > tol) {
599 BOOST_FAIL("numerical derivative in dxdy deviation is too big"
600 << "\n found at " << x << " " << y << " " << z);
601 }
602
603 if (std::fabs(x: t[i]-u[i]) > 1e5*QL_EPSILON) {
604 BOOST_FAIL("numerical derivative in dxdz not equal to dzdx"
605 << "\n found at " << x << " " << y << " " << z
606 << "\n value " << std::fabs(t[i]-u[i]));
607 }
608 }
609
610 t = SecondOrderMixedDerivativeOp(1, 2, mesher).apply(r);
611 u = SecondOrderMixedDerivativeOp(2, 1, mesher).apply(r);
612 for (const auto& iter : *index) {
613 const Size i = iter.index();
614 const Real x = mesher->location(iter, direction: 0);
615 const Real y = mesher->location(iter, direction: 1);
616 const Real z = mesher->location(iter, direction: 2);
617
618 const Real d = -std::sin(x: x)*std::sin(x: y)*std::exp(x: z);
619
620 if (std::fabs(x: d - t[i]) > tol) {
621 BOOST_FAIL("numerical derivative in dydz deviation is too big"
622 << "\n found at " << x << " " << y << " " << z);
623 }
624
625 if (std::fabs(x: t[i]-u[i]) > 1e5*QL_EPSILON) {
626 BOOST_FAIL("numerical derivative in dydz not equal to dzdy"
627 << "\n found at " << x << " " << y << " " << z
628 << "\n value " << std::fabs(t[i]-u[i]));
629 }
630 }
631
632
633}
634
635void FdmLinearOpTest::testTripleBandMapSolve() {
636
637 BOOST_TEST_MESSAGE("Testing triple-band map solution...");
638
639 const std::vector<Size> dim = {100, 400};
640
641 ext::shared_ptr<FdmLinearOpLayout> layout(new FdmLinearOpLayout(dim));
642
643 std::vector<std::pair<Real, Real> > boundaries = {{0, 1.0}, {0, 1.0}};
644
645 ext::shared_ptr<FdmMesher> mesher(
646 new UniformGridMesher(layout, boundaries));
647
648 FirstDerivativeOp dy(1, mesher);
649 dy.axpyb(a: Array(1, 2.0), x: dy, y: dy, b: Array(1, 1.0));
650
651 // check copy constructor
652 FirstDerivativeOp copyOfDy(dy);
653
654 Array u(layout->size());
655 for (Size i=0; i < layout->size(); ++i)
656 u[i] = std::sin(x: 0.1*i)+std::cos(x: 0.35*i);
657
658 Array t(dy.solve_splitting(r: copyOfDy.apply(r: u), a: 1.0, b: 0.0));
659 for (Size i=0; i < u.size(); ++i) {
660 if (std::fabs(x: u[i] - t[i]) > 1e-6) {
661 BOOST_FAIL("solve and apply are not consistent "
662 << "\n expected : " << u[i]
663 << "\n calculated : " << t[i]);
664 }
665 }
666
667 FirstDerivativeOp dx(0, mesher);
668 dx.axpyb(a: Array(), x: dx, y: dx, b: Array(1, 1.0));
669
670 FirstDerivativeOp copyOfDx(0, mesher);
671 // check assignment
672 copyOfDx = dx;
673
674 t = dx.solve_splitting(r: copyOfDx.apply(r: u), a: 1.0, b: 0.0);
675 for (Size i=0; i < u.size(); ++i) {
676 if (std::fabs(x: u[i] - t[i]) > 1e-6) {
677 BOOST_FAIL("solve and apply are not consistent "
678 << "\n expected : " << u[i]
679 << "\n calculated : " << t[i]);
680 }
681 }
682
683 SecondDerivativeOp dxx(0, mesher);
684 dxx.axpyb(a: Array(1, 0.5), x: dxx, y: dx, b: Array(1, 1.0));
685
686 // check of copy constructor
687 SecondDerivativeOp copyOfDxx(dxx);
688
689 t = dxx.solve_splitting(r: copyOfDxx.apply(r: u), a: 1.0, b: 0.0);
690
691 for (Size i=0; i < u.size(); ++i) {
692 if (std::fabs(x: u[i] - t[i]) > 1e-6) {
693 BOOST_FAIL("solve and apply are not consistent "
694 << "\n expected : " << u[i]
695 << "\n calculated : " << t[i]);
696 }
697 }
698
699 //check assignment operator
700 copyOfDxx.add(m: SecondDerivativeOp(1, mesher));
701 copyOfDxx = dxx;
702
703 t = dxx.solve_splitting(r: copyOfDxx.apply(r: u), a: 1.0, b: 0.0);
704
705 for (Size i=0; i < u.size(); ++i) {
706 if (std::fabs(x: u[i] - t[i]) > 1e-6) {
707 BOOST_FAIL("solve and apply are not consistent "
708 << "\n expected : " << u[i]
709 << "\n calculated : " << t[i]);
710 }
711 }
712}
713
714
715void FdmLinearOpTest::testFdmHestonBarrier() {
716
717 BOOST_TEST_MESSAGE("Testing FDM with barrier option in Heston model...");
718
719 const std::vector<Size> dim = {200, 100};
720
721 ext::shared_ptr<FdmLinearOpLayout> index(new FdmLinearOpLayout(dim));
722
723 std::vector<std::pair<Real, Real> > boundaries = {{3.8, 4.905274778}, {0.0, 1.0}};
724
725 ext::shared_ptr<FdmMesher> mesher(
726 new UniformGridMesher(index, boundaries));
727
728 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
729
730 Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc: Actual365Fixed()));
731 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual365Fixed()));
732
733 ext::shared_ptr<HestonProcess> hestonProcess(
734 new HestonProcess(rTS, qTS, s0, 0.04, 2.5, 0.04, 0.66, -0.8));
735
736 Settings::instance().evaluationDate() = Date(28, March, 2004);
737 Date exerciseDate(28, March, 2005);
738
739 ext::shared_ptr<FdmLinearOpComposite> hestonOp(
740 new FdmHestonOp(mesher, hestonProcess));
741
742 Array rhs(mesher->layout()->size());
743 for (const auto& iter : *mesher->layout()) {
744 rhs[iter.index()]=std::max(a: std::exp(x: mesher->location(iter,direction: 0))-100, b: 0.0);
745 }
746
747 FdmBoundaryConditionSet bcSet = {
748 ext::make_shared<FdmDirichletBoundary>(args&: mesher, args: 0.0, args: 0,
749 args: FdmDirichletBoundary::Upper)
750 };
751
752 const Real theta=0.5+std::sqrt(x: 3.0)/6.;
753 HundsdorferScheme hsEvolver(theta, 0.5, hestonOp, bcSet);
754 FiniteDifferenceModel<HundsdorferScheme> hsModel(hsEvolver);
755 hsModel.rollback(a&: rhs, from: 1.0, to: 0.0, steps: 50);
756
757 Matrix ret(dim[0], dim[1]);
758 for (Size i=0; i < dim[0]; ++i)
759 for (Size j=0; j < dim[1]; ++j)
760 ret[i][j] = rhs[i+j*dim[0]];
761
762 std::vector<Real> tx, ty;
763 for (const auto& iter : *mesher->layout()) {
764 if (iter.coordinates()[1] == 0) {
765 tx.push_back(x: mesher->location(iter, direction: 0));
766 }
767 if (iter.coordinates()[0] == 0) {
768 ty.push_back(x: mesher->location(iter, direction: 1));
769 }
770 }
771
772 BilinearInterpolation interpolate(ty.begin(), ty.end(),
773 tx.begin(), tx.end(), ret);
774
775 const Real x = 100;
776 const Real v0 = 0.04;
777
778 const Real npv = interpolate(v0, std::log(x: x));
779 const Real delta = 0.5*(
780 interpolate(v0, std::log(x: x+1))-interpolate(v0, std::log(x: x-1)));
781 const Real gamma = interpolate(v0, std::log(x: x+1))
782 + interpolate(v0, std::log(x: x-1)) - 2*npv;
783
784 const Real npvExpected = 9.049016;
785 const Real deltaExpected = 0.511285;
786 const Real gammaExpected = -0.034296;
787
788 if (std::fabs(x: npv - npvExpected) > 0.000001) {
789 BOOST_FAIL("Error in calculating PV for Heston barrier option");
790 }
791
792 if (std::fabs(x: delta - deltaExpected) > 0.000001) {
793 BOOST_FAIL("Error in calculating Delta for Heston barrier option");
794 }
795
796 if (std::fabs(x: gamma - gammaExpected) > 0.000001) {
797 BOOST_FAIL("Error in calculating Gamma for Heston barrier option");
798 }
799}
800
801void FdmLinearOpTest::testFdmHestonAmerican() {
802
803 BOOST_TEST_MESSAGE("Testing FDM with American option in Heston model...");
804
805 const std::vector<Size> dim = {200, 100};
806
807 ext::shared_ptr<FdmLinearOpLayout> index(new FdmLinearOpLayout(dim));
808
809 std::vector<std::pair<Real, Real> > boundaries = {{3.8, std::log(x: 220.0)}, {0.0, 1.0}};
810
811 ext::shared_ptr<FdmMesher> mesher(
812 new UniformGridMesher(index, boundaries));
813
814 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
815
816 Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc: Actual365Fixed()));
817 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual365Fixed()));
818
819 ext::shared_ptr<HestonProcess> hestonProcess(
820 new HestonProcess(rTS, qTS, s0, 0.04, 2.5, 0.04, 0.66, -0.8));
821
822 Settings::instance().evaluationDate() = Date(28, March, 2004);
823 Date exerciseDate(28, March, 2005);
824
825 ext::shared_ptr<FdmLinearOpComposite> LinearOp(
826 new FdmHestonOp(mesher, hestonProcess));
827
828 ext::shared_ptr<Payoff> payoff(new PlainVanillaPayoff(Option::Put, 100.0));
829 Array rhs(mesher->layout()->size());
830 for (const auto& iter : *mesher->layout()) {
831 rhs[iter.index()]
832 = payoff->operator ()(price: std::exp(x: mesher->location(iter, direction: 0)));
833 }
834
835 FdmAmericanStepCondition condition(mesher,
836 ext::shared_ptr<FdmInnerValueCalculator>(
837 new FdmLogInnerValue(payoff, mesher, 0)));
838 const Real theta=0.5+std::sqrt(x: 3.0)/6.;
839 HundsdorferScheme hsEvolver(theta, 0.5, LinearOp);
840 FiniteDifferenceModel<HundsdorferScheme> hsModel(hsEvolver);
841 hsModel.rollback(a&: rhs, from: 1.0, to: 0.0, steps: 50, condition);
842
843 Matrix ret(dim[0], dim[1]);
844 for (Size i=0; i < dim[0]; ++i)
845 for (Size j=0; j < dim[1]; ++j)
846 ret[i][j] = rhs[i+j*dim[0]];
847
848 std::vector<Real> tx, ty;
849 for (const auto& iter : *mesher->layout()) {
850 if (iter.coordinates()[1] == 0) {
851 tx.push_back(x: mesher->location(iter, direction: 0));
852 }
853 if (iter.coordinates()[0] == 0) {
854 ty.push_back(x: mesher->location(iter, direction: 1));
855 }
856 }
857
858 BilinearInterpolation interpolate(ty.begin(), ty.end(),
859 tx.begin(), tx.end(), ret);
860
861 const Real x = 100;
862 const Real v0 = 0.04;
863
864 const Real npv = interpolate(v0, std::log(x: x));
865 const Real npvExpected = 5.641648;
866
867 if (std::fabs(x: npv - npvExpected) > 0.000001) {
868 BOOST_FAIL("Error in calculating PV for Heston American Option");
869 }
870}
871
872void FdmLinearOpTest::testFdmHestonExpress() {
873
874 BOOST_TEST_MESSAGE("Testing FDM with express certificate in Heston model...");
875
876 const std::vector<Size> dim = {200, 100};
877
878 ext::shared_ptr<FdmLinearOpLayout> index(new FdmLinearOpLayout(dim));
879
880 std::vector<std::pair<Real, Real> > boundaries = {{3.8, std::log(x: 220.0)}, {0.0, 1.0}};
881
882 ext::shared_ptr<FdmMesher> mesher(
883 new UniformGridMesher(index, boundaries));
884
885 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
886
887 Handle<YieldTermStructure> rTS(flatRate(forward: 0.05, dc: Actual365Fixed()));
888 Handle<YieldTermStructure> qTS(flatRate(forward: 0.0 , dc: Actual365Fixed()));
889
890 Handle<HestonProcess> hestonProcess(ext::make_shared<HestonProcess> (
891 args&: rTS, args&: qTS, args&: s0, args: 0.04, args: 2.5, args: 0.04, args: 0.66, args: -0.8));
892
893 const Date exerciseDate(28, March, 2005);
894 const Date evaluationDate(28, March, 2004);
895 Settings::instance().evaluationDate() = evaluationDate;
896
897 std::vector<Real> triggerLevels(2);
898 triggerLevels[0] = triggerLevels[1] = 100.0;
899 std::vector<Real> redemptions(2);
900 redemptions[0] = redemptions[1] = 108.0;
901 std::vector<Time> exerciseTimes(2);
902 exerciseTimes[0] = 0.333; exerciseTimes[1] = 0.666;
903
904 DividendSchedule dividendSchedule(1, ext::shared_ptr<Dividend>(
905 new FixedDividend(2.5, evaluationDate + Period(6, Months))));
906 ext::shared_ptr<FdmDividendHandler> dividendCondition(
907 new FdmDividendHandler(dividendSchedule, mesher,
908 rTS->referenceDate(),
909 rTS->dayCounter(), 0));
910
911 ext::shared_ptr<StepCondition<Array> > expressCondition(
912 new FdmHestonExpressCondition(redemptions, triggerLevels,
913 exerciseTimes, mesher));
914
915 std::list<std::vector<Time>> stoppingTimes = {
916 exerciseTimes, dividendCondition->dividendTimes()
917 };
918
919 std::list<ext::shared_ptr<StepCondition<Array>>> conditions = {
920 expressCondition, dividendCondition
921 };
922
923 ext::shared_ptr<FdmStepConditionComposite> condition(
924 new FdmStepConditionComposite(stoppingTimes, conditions));
925
926 ext::shared_ptr<Payoff> payoff(new ExpressPayoff());
927
928 ext::shared_ptr<FdmInnerValueCalculator> calculator(
929 new FdmLogInnerValue(payoff, mesher, 0));
930
931 const FdmBoundaryConditionSet bcSet;
932 const FdmSolverDesc solverDesc = { .mesher: mesher, .bcSet: bcSet,
933 .condition: condition, .calculator: calculator, .maturity: 1.0, .timeSteps: 50, .dampingSteps: 0 };
934 FdmHestonSolver solver(hestonProcess, solverDesc);
935
936 const Real s = s0->value();
937 const Real v0 = 0.04;
938
939 if (std::fabs(x: solver.valueAt(s, v: v0) - 101.027) > 0.01) {
940 BOOST_FAIL("Error in calculating PV for Heston Express Certificate");
941 }
942
943 if (std::fabs(x: solver.deltaAt(s, v: v0) - 0.4181) > 0.001) {
944 BOOST_FAIL("Error in calculating Delta for Heston Express Certificate");
945 }
946
947 if (std::fabs(x: solver.gammaAt(s, v: v0) + 0.0400) > 0.001) {
948 BOOST_FAIL("Error in calculating Gamma for Heston Express Certificate");
949 }
950
951 if (std::fabs(x: solver.meanVarianceDeltaAt(s, v: v0) - 0.6602) > 0.001) {
952 BOOST_FAIL("Error in calculating mean variance Delta for "
953 "Heston Express Certificate");
954 }
955
956 if (std::fabs(x: solver.meanVarianceGammaAt(s, v: v0) + 0.0316) > 0.001) {
957 BOOST_FAIL("Error in calculating mean variance Delta for "
958 "Heston Express Certificate");
959 }
960}
961
962
963namespace {
964
965 ext::shared_ptr<HybridHestonHullWhiteProcess> createHestonHullWhite(
966 Time maturity) {
967
968 DayCounter dc = Actual365Fixed();
969 const Date today = Settings::instance().evaluationDate();
970 Handle<Quote> s0(ext::shared_ptr<Quote>(new SimpleQuote(100.0)));
971
972 std::vector<Date> dates;
973 std::vector<Rate> rates, divRates;
974
975 for (Size i=0; i <= 25; ++i) {
976 dates.push_back(x: today+Period(i, Years));
977 rates.push_back(x: 0.05);
978 divRates.push_back(x: 0.02);
979 }
980
981 const Handle<YieldTermStructure> rTS(
982 ext::shared_ptr<YieldTermStructure>(new ZeroCurve(dates, rates, dc)));
983 const Handle<YieldTermStructure> qTS(
984 ext::shared_ptr<YieldTermStructure>(
985 new ZeroCurve(dates, divRates, dc)));
986
987 const Real v0 = 0.04;
988 ext::shared_ptr<HestonProcess> hestonProcess(
989 new HestonProcess(rTS, qTS, s0, v0, 1.0, v0*0.75, 0.4, -0.7));
990
991 ext::shared_ptr<HullWhiteForwardProcess> hwFwdProcess(
992 new HullWhiteForwardProcess(rTS, 0.00883, 0.01));
993 hwFwdProcess->setForwardMeasureTime(maturity);
994
995 const Real equityShortRateCorr = -0.7;
996
997 return ext::make_shared<HybridHestonHullWhiteProcess>(
998 args&: hestonProcess, args&: hwFwdProcess,
999 args: equityShortRateCorr);
1000 }
1001
1002 FdmSolverDesc createSolverDesc(
1003 const std::vector<Size>& dim,
1004 const ext::shared_ptr<HybridHestonHullWhiteProcess>& process) {
1005
1006 const Time maturity
1007 = process->hullWhiteProcess()->getForwardMeasureTime();
1008
1009 ext::shared_ptr<FdmLinearOpLayout> layout(new FdmLinearOpLayout(dim));
1010
1011 std::vector<ext::shared_ptr<Fdm1dMesher> > mesher1d = {
1012 ext::shared_ptr<Fdm1dMesher>(
1013 new Uniform1dMesher(std::log(x: 22.0), std::log(x: 440.0), dim[0])),
1014 ext::shared_ptr<Fdm1dMesher>(
1015 new FdmHestonVarianceMesher(dim[1], process->hestonProcess(),
1016 maturity)),
1017 ext::shared_ptr<Fdm1dMesher>(
1018 new Uniform1dMesher(-0.15, 0.15, dim[2]))
1019 };
1020
1021 const ext::shared_ptr<FdmMesher> mesher(
1022 new FdmMesherComposite(mesher1d));
1023
1024 const FdmBoundaryConditionSet boundaries;
1025
1026 std::list<std::vector<Time> > stoppingTimes;
1027 std::list<ext::shared_ptr<StepCondition<Array> > > stepConditions;
1028
1029 ext::shared_ptr<FdmStepConditionComposite> conditions(
1030 new FdmStepConditionComposite(
1031 std::list<std::vector<Time> >(),
1032 FdmStepConditionComposite::Conditions()));
1033
1034 ext::shared_ptr<StrikedTypePayoff> payoff(
1035 new PlainVanillaPayoff(Option::Call, 160.0));
1036
1037 ext::shared_ptr<FdmInnerValueCalculator> calculator(
1038 new FdmLogInnerValue(payoff, mesher, 0));
1039
1040 const Size tGrid = 100;
1041 const Size dampingSteps = 0;
1042
1043 FdmSolverDesc desc = { .mesher: mesher, .bcSet: boundaries,
1044 .condition: conditions, .calculator: calculator,
1045 .maturity: maturity, .timeSteps: tGrid, .dampingSteps: dampingSteps };
1046
1047 return desc;
1048 }
1049}
1050
1051void FdmLinearOpTest::testFdmHestonHullWhiteOp() {
1052 BOOST_TEST_MESSAGE("Testing FDM with Heston Hull-White model...");
1053
1054 const Date today = Date(28, March, 2004);
1055 Settings::instance().evaluationDate() = today;
1056
1057 Date exerciseDate(28, March, 2012);
1058 const Time maturity = Actual365Fixed().yearFraction(d1: today, d2: exerciseDate);
1059
1060 const std::vector<Size> dim = {51, 31, 31};
1061
1062 ext::shared_ptr<HybridHestonHullWhiteProcess> jointProcess
1063 = createHestonHullWhite(maturity);
1064 FdmSolverDesc desc = createSolverDesc(dim, process: jointProcess);
1065 ext::shared_ptr<FdmMesher> mesher = desc.mesher;
1066
1067 ext::shared_ptr<HullWhiteForwardProcess> hwFwdProcess
1068 = jointProcess->hullWhiteProcess();
1069
1070 ext::shared_ptr<HullWhiteProcess> hwProcess(
1071 new HullWhiteProcess(jointProcess->hestonProcess()->riskFreeRate(),
1072 hwFwdProcess->a(), hwFwdProcess->sigma()));
1073
1074 ext::shared_ptr<FdmLinearOpComposite> linearOp(
1075 new FdmHestonHullWhiteOp(mesher,
1076 jointProcess->hestonProcess(),
1077 hwProcess,
1078 jointProcess->eta()));
1079
1080 Array rhs(mesher->layout()->size());
1081 for (const auto& iter : *mesher->layout()) {
1082 rhs[iter.index()] = desc.calculator->avgInnerValue(iter, t: maturity);
1083 }
1084
1085 const Real theta = 0.5+std::sqrt(x: 3.0)/6.;
1086 HundsdorferScheme hsEvolver(theta, 0.5, linearOp);
1087 FiniteDifferenceModel<HundsdorferScheme> hsModel(hsEvolver);
1088
1089 hsModel.rollback(a&: rhs, from: maturity, to: 0.0, steps: desc.timeSteps);
1090
1091 std::vector<Real> tx, ty, tr, y;
1092 for (const auto& iter : *mesher->layout()) {
1093 if (iter.coordinates()[1] == 0 && iter.coordinates()[2] == 0) {
1094 tx.push_back(x: mesher->location(iter, direction: 0));
1095 }
1096 if (iter.coordinates()[0] == 0 && iter.coordinates()[2] == 0) {
1097 ty.push_back(x: mesher->location(iter, direction: 1));
1098 }
1099 if (iter.coordinates()[0] == 0 && iter.coordinates()[1] == 0) {
1100 tr.push_back(x: mesher->location(iter, direction: 2));
1101 }
1102 }
1103
1104 const Real x0 = 100;
1105 const Real v0 = jointProcess->hestonProcess()->v0();
1106 const Real r0 = 0;
1107 for (Size k=0; k < dim[2]; ++k) {
1108 Matrix ret(dim[0], dim[1]);
1109 for (Size i=0; i < dim[0]; ++i)
1110 for (Size j=0; j < dim[1]; ++j)
1111 ret[i][j] = rhs[ i+j*dim[0]+k*dim[0]*dim[1] ];
1112
1113 y.push_back(x: BicubicSpline(ty.begin(), ty.end(),
1114 tx.begin(), tx.end(), ret)(v0, std::log(x: x0)));
1115 }
1116
1117 const Real directCalc
1118 = MonotonicCubicNaturalSpline(tr.begin(), tr.end(), y.begin())(r0);
1119
1120 std::vector<Real> x(3);
1121 x[0] = std::log(x: x0); x[1] = v0; x[2] = r0;
1122
1123 Fdm3DimSolver solver3d(desc, FdmSchemeDesc::Hundsdorfer(), linearOp);
1124 const Real solverCalc = solver3d.interpolateAt(x: x[0], y: x[1], z: x[2]);
1125 const Real solverTheta = solver3d.thetaAt(x: x[0], y: x[1], z: x[2]);
1126
1127 if (std::fabs(x: directCalc - solverCalc) > 1e-4) {
1128 BOOST_FAIL("Error in calculating PV for Heston Hull White Option");
1129 }
1130
1131
1132 FdmNdimSolver<3> solverNd(desc, FdmSchemeDesc::Hundsdorfer(), linearOp);
1133 const Real solverNdCalc = solverNd.interpolateAt(x);
1134 const Real solverNdTheta = solverNd.thetaAt(x);
1135
1136 if (std::fabs(x: solverNdCalc - solverCalc) > 1e-4) {
1137 BOOST_FAIL("Error in calculating PV for Heston Hull White Option");
1138 }
1139 if (std::fabs(x: solverNdTheta - solverTheta) > 1e-4) {
1140 BOOST_FAIL("Error in calculating PV for Heston Hull White Option");
1141 }
1142
1143 VanillaOption option(
1144 ext::shared_ptr<StrikedTypePayoff>(
1145 new PlainVanillaPayoff(Option::Call, 160.0)),
1146 ext::shared_ptr<Exercise>(new EuropeanExercise(exerciseDate)));
1147
1148 const Real tol = 0.025;
1149 option.setPricingEngine(
1150 MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
1151 .withSteps(steps: 200)
1152 .withAntitheticVariate()
1153 .withControlVariate()
1154 .withAbsoluteTolerance(tolerance: tol)
1155 .withSeed(seed: 42));
1156
1157 // the following takes far too long
1158 // const Real expected = option.NPV();
1159 // use precalculated value instead
1160 const Real expected = 4.73;
1161
1162 if (std::fabs(x: directCalc - expected) > 3*tol) {
1163 BOOST_FAIL("Error in calculating MC PV for Heston Hull White Option");
1164 }
1165}
1166
1167namespace {
1168 Array axpy(const boost::numeric::ublas::compressed_matrix<Real>& A,
1169 const Array& x) {
1170
1171 boost::numeric::ublas::vector<Real> tmpX(x.size()), y(x.size());
1172 std::copy(first: x.begin(), last: x.end(), result: tmpX.begin());
1173 boost::numeric::ublas::axpy_prod(e1: A, e2: tmpX, v&: y);
1174
1175 return Array(y.begin(), y.end());
1176 }
1177
1178 boost::numeric::ublas::compressed_matrix<Real> createTestMatrix(
1179 Size n, Size m, Real theta) {
1180
1181 boost::numeric::ublas::compressed_matrix<Real> a(n*m, n*m);
1182
1183 for (Size i=0; i < n; ++i) {
1184 for (Size j=0; j < m; ++j) {
1185 const Size k = i*m+j;
1186 a(k,k)=1.0;
1187
1188 if (i > 0 && j > 0 && i <n-1 && j < m-1) {
1189 const Size im1 = i-1;
1190 const Size ip1 = i+1;
1191 const Size jm1 = j-1;
1192 const Size jp1 = j+1;
1193 const Real delta = theta/((ip1-im1)*(jp1-jm1));
1194
1195 a(k,im1*m+jm1) = delta;
1196 a(k,im1*m+jp1) = -delta;
1197 a(k,ip1*m+jm1) = -delta;
1198 a(k,ip1*m+jp1) = delta;
1199 }
1200 }
1201 }
1202
1203 return a;
1204 }
1205}
1206
1207void FdmLinearOpTest::testBiCGstab() {
1208 BOOST_TEST_MESSAGE(
1209 "Testing bi-conjugated gradient stabilized algorithm...");
1210
1211 const Size n=41, m=21;
1212 const Real theta = 1.0;
1213 const boost::numeric::ublas::compressed_matrix<Real> a
1214 = createTestMatrix(n, m, theta);
1215
1216 const ext::function<Array(const Array&)> matmult
1217 = [&](const Array& _x) { return axpy(A: a, x: _x); };
1218
1219 SparseILUPreconditioner ilu(a, 4);
1220 ext::function<Array(const Array&)> precond
1221 = [&](const Array& _x) { return ilu.apply(b: _x); };
1222
1223 Array b(n*m);
1224 MersenneTwisterUniformRng rng(1234);
1225 for (Real& i : b) {
1226 i = rng.next().value;
1227 }
1228
1229 const Real tol = 1e-10;
1230
1231 const BiCGstab biCGstab(matmult, n*m, tol, precond);
1232 const Array x = biCGstab.solve(b).x;
1233
1234 const Real error = std::sqrt(x: DotProduct(v1: b-axpy(A: a, x),
1235 v2: b-axpy(A: a, x))/DotProduct(v1: b,v2: b));
1236
1237 if (error > tol) {
1238 BOOST_FAIL("Error calculating the inverse using BiCGstab" <<
1239 "\n tolerance: " << tol <<
1240 "\n error: " << error);
1241 }
1242}
1243
1244void FdmLinearOpTest::testGMRES() {
1245 BOOST_TEST_MESSAGE("Testing GMRES algorithm...");
1246
1247 const Size n=41, m=21;
1248 const Real theta = 1.0;
1249 const boost::numeric::ublas::compressed_matrix<Real> a
1250 = createTestMatrix(n, m, theta);
1251
1252 const ext::function<Array(const Array&)> matmult
1253 = [&](const Array& _x) { return axpy(A: a, x: _x); };
1254
1255 SparseILUPreconditioner ilu(a, 4);
1256 ext::function<Array(const Array&)> precond
1257 = [&](const Array& _x) { return ilu.apply(b: _x); };
1258
1259 Array b(n*m);
1260 MersenneTwisterUniformRng rng(1234);
1261 for (Real& i : b) {
1262 i = rng.next().value;
1263 }
1264
1265 const Real tol = 1e-10;
1266
1267 const GMRES gmres(matmult, n*m, tol, precond);
1268 const GMRESResult result = gmres.solve(b, x0: b);
1269 const Array x = result.x;
1270 const Real errorCalculated = result.errors.back();
1271
1272 const Real error = std::sqrt(x: DotProduct(v1: b-axpy(A: a, x),
1273 v2: b-axpy(A: a, x))/DotProduct(v1: b,v2: b));
1274
1275 if (error > tol) {
1276 BOOST_FAIL("Error calculating the inverse using GMRES" <<
1277 "\n tolerance: " << tol <<
1278 "\n error: " << error);
1279 }
1280
1281 if (std::fabs(x: error - errorCalculated) > 10*QL_EPSILON) {
1282 BOOST_FAIL("Calculation if the error in GMRES went wrong" <<
1283 "\n calculated: " << errorCalculated <<
1284 "\n error: " << error);
1285
1286 }
1287
1288 const GMRES gmresRestart(matmult, 5, tol, precond);
1289 const GMRESResult resultRestart = gmresRestart.solveWithRestart(restart: 5, b, x0: b);
1290 const Real errorWithRestart = resultRestart.errors.back();
1291
1292 if (errorWithRestart > tol) {
1293 BOOST_FAIL("Error calculating the inverse using "
1294 "GMRES with restarts" <<
1295 "\n tolerance: " << tol <<
1296 "\n error: " << errorWithRestart);
1297 }
1298}
1299
1300void FdmLinearOpTest::testCrankNicolsonWithDamping() {
1301
1302 BOOST_TEST_MESSAGE("Testing Crank-Nicolson with initial implicit damping steps "
1303 "for a digital option...");
1304
1305 DayCounter dc = Actual360();
1306 Date today = Date::todaysDate();
1307
1308 ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(100.0));
1309 ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, forward: 0.06, dc);
1310 ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: 0.06, dc);
1311 ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, volatility: 0.35, dc);
1312
1313 ext::shared_ptr<StrikedTypePayoff> payoff(
1314 new CashOrNothingPayoff(Option::Put, 100, 10.0));
1315
1316 Time maturity = 0.75;
1317 Date exDate = today + timeToDays(t: maturity);
1318 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(exDate));
1319
1320 ext::shared_ptr<BlackScholesMertonProcess> process(new
1321 BlackScholesMertonProcess(Handle<Quote>(spot),
1322 Handle<YieldTermStructure>(qTS),
1323 Handle<YieldTermStructure>(rTS),
1324 Handle<BlackVolTermStructure>(volTS)));
1325 ext::shared_ptr<PricingEngine> engine(
1326 new AnalyticEuropeanEngine(process));
1327
1328 VanillaOption opt(payoff, exercise);
1329 opt.setPricingEngine(engine);
1330 Real expectedPV = opt.NPV();
1331 Real expectedGamma = opt.gamma();
1332
1333 // fd pricing using implicit damping steps and Crank Nicolson
1334 const Size csSteps = 25, dampingSteps = 3, xGrid = 400;
1335 const std::vector<Size> dim(1, xGrid);
1336
1337 ext::shared_ptr<FdmLinearOpLayout> layout(new FdmLinearOpLayout(dim));
1338 const ext::shared_ptr<Fdm1dMesher> equityMesher(
1339 new FdmBlackScholesMesher(
1340 dim[0], process, maturity, payoff->strike(),
1341 Null<Real>(), Null<Real>(), 0.0001, 1.5,
1342 std::pair<Real, Real>(payoff->strike(), 0.01)));
1343
1344 const ext::shared_ptr<FdmMesher> mesher (
1345 new FdmMesherComposite(equityMesher));
1346
1347 ext::shared_ptr<FdmBlackScholesOp> map(
1348 new FdmBlackScholesOp(mesher, process, payoff->strike()));
1349
1350 ext::shared_ptr<FdmInnerValueCalculator> calculator(
1351 new FdmLogInnerValue(payoff, mesher, 0));
1352
1353 Array rhs(layout->size()), x(layout->size());
1354
1355 for (const auto& iter : *layout) {
1356 rhs[iter.index()] = calculator->avgInnerValue(iter, t: maturity);
1357 x[iter.index()] = mesher->location(iter, direction: 0);
1358 }
1359
1360 FdmBackwardSolver solver(map, FdmBoundaryConditionSet(),
1361 ext::shared_ptr<FdmStepConditionComposite>(),
1362 FdmSchemeDesc::Douglas());
1363 solver.rollback(a&: rhs, from: maturity, to: 0.0, steps: csSteps, dampingSteps);
1364
1365 MonotonicCubicNaturalSpline spline(x.begin(), x.end(), rhs.begin());
1366
1367 Real s = spot->value();
1368 Real calculatedPV = spline(std::log(x: s));
1369 Real calculatedGamma = (spline.secondDerivative(x: std::log(x: s))
1370 - spline.derivative(x: std::log(x: s)) )/(s*s);
1371
1372 Real relTol = 2e-3;
1373
1374 if (std::fabs(x: calculatedPV - expectedPV) > relTol*expectedPV) {
1375 BOOST_FAIL("Error calculating the PV of the digital option" <<
1376 "\n rel. tolerance: " << relTol <<
1377 "\n expected: " << expectedPV <<
1378 "\n calculated: " << calculatedPV);
1379 }
1380 if (std::fabs(x: calculatedGamma - expectedGamma) > relTol*expectedGamma) {
1381 BOOST_FAIL("Error calculating the Gamma of the digital option" <<
1382 "\n rel. tolerance: " << relTol <<
1383 "\n expected: " << expectedGamma <<
1384 "\n calculated: " << calculatedGamma);
1385 }
1386}
1387
1388void FdmLinearOpTest::testSpareMatrixReference() {
1389 BOOST_TEST_MESSAGE("Testing SparseMatrixReference type...");
1390
1391 const Size rows = 10;
1392 const Size columns = 10;
1393 const Size nMatrices = 5;
1394 const Size nElements = 50;
1395
1396 PseudoRandom::urng_type rng(1234UL);
1397
1398 SparseMatrix expected(rows, columns);
1399 std::vector<SparseMatrix> v(nMatrices, SparseMatrix(rows, columns));
1400 std::vector<SparseMatrixReference> refs;
1401
1402 for (auto& i : v) {
1403 SparseMatrixReference m(i);
1404 for (Size j=0; j < nElements; ++j) {
1405 const Size row = Size(rng.next().value*rows);
1406 const Size column = Size(rng.next().value*columns);
1407
1408 const Real value = rng.next().value;
1409 m(row, column) += value;
1410 expected(row, column) += value;
1411 }
1412
1413 refs.push_back(x: m);
1414 }
1415
1416 SparseMatrix calculated = std::accumulate(first: refs.begin()+1, last: refs.end(),
1417 init: SparseMatrix(refs.front()));
1418
1419 for (Size i=0; i < rows; ++i) {
1420 for (Size j=0; j < columns; ++j) {
1421 if (std::fabs(x: Real(calculated(i,j)) - Real(expected(i,j))) > 100*QL_EPSILON) {
1422 BOOST_FAIL("Error using sparse matrix references in " <<
1423 "Element (" << i << ", " << j << ")" <<
1424 "\n expected : " << Real(expected(i, j)) <<
1425 "\n calculated: " << Real(calculated(i, j)));
1426 }
1427 }
1428 }
1429}
1430
1431namespace {
1432
1433 Size nrElementsOfSparseMatrix(const SparseMatrix& m) {
1434 Size retVal = 0;
1435 for (SparseMatrix::const_iterator1 i1 = m.begin1();
1436 i1 != m.end1(); ++i1) {
1437 retVal+=std::distance(first: i1.begin(), last: i1.end());
1438 }
1439 return retVal;
1440 }
1441
1442}
1443
1444void FdmLinearOpTest::testSparseMatrixZeroAssignment() {
1445 BOOST_TEST_MESSAGE("Testing assignment to zero in sparse matrix...");
1446
1447 SparseMatrix m(5,5);
1448 if (nrElementsOfSparseMatrix(m) != 0U) {
1449 BOOST_FAIL("non zero return for an emtpy matrix");
1450 }
1451 m(0, 0) = 0.0; m(1, 2) = 0.0;
1452 if (nrElementsOfSparseMatrix(m) != 2) {
1453 BOOST_FAIL("two elements expected");
1454 }
1455 m(1, 3) = 1.0;
1456 if (nrElementsOfSparseMatrix(m) != 3) {
1457 BOOST_FAIL("three elements expected");
1458 }
1459 m(1, 3) = 0.0;
1460 if (nrElementsOfSparseMatrix(m) != 3) {
1461 BOOST_FAIL("three elements expected");
1462 }
1463}
1464
1465void FdmLinearOpTest::testFdmMesherIntegral() {
1466 BOOST_TEST_MESSAGE("Testing integrals over meshers functions...");
1467
1468 const ext::shared_ptr<FdmMesherComposite> mesher(
1469 new FdmMesherComposite(
1470 ext::shared_ptr<Fdm1dMesher>(new Concentrating1dMesher(
1471 -1, 1.6, 21, std::pair<Real, Real>(0, 0.1))),
1472 ext::shared_ptr<Fdm1dMesher>(new Concentrating1dMesher(
1473 -3, 4, 11, std::pair<Real, Real>(1, 0.01))),
1474 ext::shared_ptr<Fdm1dMesher>(new Concentrating1dMesher(
1475 -2, 1, 5, std::pair<Real, Real>(0.5, 0.1)))));
1476
1477 Array f(mesher->layout()->size());
1478 for (const auto& iter : *mesher->layout()) {
1479 const Real x = mesher->location(iter, direction: 0);
1480 const Real y = mesher->location(iter, direction: 1);
1481 const Real z = mesher->location(iter, direction: 2);
1482
1483 f[iter.index()] = x*x + 3*y*y - 3*z*z
1484 + 2*x*y - x*z - 3*y*z
1485 + 4*x - y - 3*z + 2 ;
1486 }
1487
1488 const Real tol = 1e-12;
1489
1490 // Simpson's rule has to be exact here, Mathematica code gives
1491 // Integrate[x*x+3*y*y-3*z*z+2*x*y-x*z-3*y*z+4*x-y-3*z+2,
1492 // {x, -1, 16/10}, {y, -3, 4}, {z, -2, 1}]
1493 const Real expectedSimpson = 876.512;
1494 const Real calculatedSimpson
1495 = FdmMesherIntegral(mesher, DiscreteSimpsonIntegral()).integrate(f);
1496
1497 if (std::fabs(x: calculatedSimpson - expectedSimpson) > tol*expectedSimpson) {
1498 BOOST_FAIL(std::setprecision(16)
1499 << "discrete mesher integration using Simpson's rule failed: "
1500 << "\n calculated: " << calculatedSimpson
1501 << "\n expected: " << expectedSimpson);
1502 }
1503
1504 const Real expectedTrapezoid = 917.0148209153263;
1505 const Real calculatedTrapezoid
1506 = FdmMesherIntegral(mesher, DiscreteTrapezoidIntegral()).integrate(f);
1507
1508 if (std::fabs(x: calculatedTrapezoid - expectedTrapezoid)
1509 > tol*expectedTrapezoid) {
1510 BOOST_FAIL(std::setprecision(16)
1511 << "discrete mesher integration using Trapezoid rule failed: "
1512 << "\n calculated: " << calculatedTrapezoid
1513 << "\n expected: " << expectedTrapezoid);
1514 }
1515}
1516
1517void FdmLinearOpTest::testHighInterestRateBlackScholesMesher() {
1518 BOOST_TEST_MESSAGE("Testing Black-Scholes mesher in a "
1519 "high interest rate scenario...");
1520
1521 const DayCounter dc = Actual365Fixed();
1522 const Date today = Date(11, February, 2018);
1523
1524 const Real spot = 100;
1525 const Rate r = 0.21;
1526 const Rate q = 0.02;
1527 const Volatility v = 0.25;
1528
1529 const ext::shared_ptr<GeneralizedBlackScholesProcess> process =
1530 ext::make_shared<GeneralizedBlackScholesProcess>(
1531 args: Handle<Quote>(ext::make_shared<SimpleQuote>(args: spot)),
1532 args: Handle<YieldTermStructure>(flatRate(today, forward: q, dc)),
1533 args: Handle<YieldTermStructure>(flatRate(today, forward: r, dc)),
1534 args: Handle<BlackVolTermStructure>(flatVol(today, volatility: v, dc)));
1535
1536 const Size size = 10;
1537 const Time maturity = 2.0;
1538 const Real strike = 100;
1539 const Real eps = 0.05;
1540 const Real normInvEps = 1.64485363;
1541 const Real scaleFactor = 2.5;
1542
1543 const std::vector<Real> loc = FdmBlackScholesMesher(
1544 size, process, maturity, strike,
1545 Null<Real>(), Null<Real>(), eps, scaleFactor).locations();
1546
1547 const Real calculatedMin = std::exp(x: loc.front());
1548 const Real calculatedMax = std::exp(x: loc.back());
1549
1550 const Real minimum = spot
1551 * std::exp(x: -normInvEps*scaleFactor*v*std::sqrt(x: maturity));
1552 const Real maximum = spot
1553 / process->riskFreeRate()->discount(t: maturity)
1554 * process->dividendYield()->discount(t: maturity)
1555 * std::exp( x: normInvEps*scaleFactor*v*std::sqrt(x: maturity));
1556
1557 const Real relTol = 1e-7;
1558
1559 const Real maxDiff = std::fabs(x: calculatedMax - maximum);
1560 if (maxDiff > relTol*maximum) {
1561 BOOST_FAIL("Upper bound for Black-Scholes mesher failed: "
1562 << "\n calculated: " << calculatedMax
1563 << "\n expected: " << maximum
1564 << std::scientific
1565 << "\n difference: " << maxDiff
1566 << "\n tolerance: " << relTol*maximum);
1567 }
1568
1569 const Real minDiff = std::fabs(x: calculatedMin - minimum);
1570 if (minDiff > relTol*minimum) {
1571 BOOST_FAIL("Lower bound for Black-Scholes mesher failed: "
1572 << "\n calculated: " << calculatedMin
1573 << "\n expected: " << minimum
1574 << std::scientific
1575 << "\n difference: " << minDiff
1576 << "\n tolerance: " << relTol*minimum);
1577 }
1578}
1579
1580void FdmLinearOpTest::testLowVolatilityHighDiscreteDividendBlackScholesMesher() {
1581 BOOST_TEST_MESSAGE("Testing Black-Scholes mesher in a low volatility and "
1582 "high discrete dividend scenario...");
1583
1584 const DayCounter dc = Actual365Fixed();
1585 const Date today = Date(28, January, 2018);
1586
1587 const Handle<Quote> spot(ext::make_shared<SimpleQuote>(args: 100.0));
1588 const Handle<YieldTermStructure> qTS(flatRate(today, forward: 0.07, dc));
1589 const Handle<YieldTermStructure> rTS(flatRate(today, forward: 0.16, dc));
1590 const Handle<BlackVolTermStructure> volTS(flatVol(today, volatility: 0.0, dc));
1591
1592 const ext::shared_ptr<GeneralizedBlackScholesProcess> process =
1593 ext::make_shared<GeneralizedBlackScholesProcess>(
1594 args: spot, args: qTS, args: rTS, args: volTS);
1595
1596 const Date firstDivDate = today + Period(7, Months);
1597 const Real firstDivAmount = 10.0;
1598 const Date secondDivDate = today + Period(11, Months);
1599 const Real secondDivAmount = 5.0;
1600
1601 DividendSchedule divSchedule = {
1602 ext::make_shared<FixedDividend>(args: firstDivAmount, args: firstDivDate),
1603 ext::make_shared<FixedDividend>(args: secondDivAmount, args: secondDivDate)
1604 };
1605
1606 const Size size = 5;
1607 const Time maturity = 1.0;
1608 const Real strike = 100;
1609 const Real eps = 0.0001;
1610 const Real scaleFactor = 1.5;
1611
1612 const std::vector<Real> loc = FdmBlackScholesMesher(
1613 size,
1614 process,
1615 maturity, strike,
1616 Null<Real>(), Null<Real>(),
1617 eps, scaleFactor,
1618 std::make_pair(x: Null<Real>(), y: Null<Real>()),
1619 divSchedule).locations();
1620
1621 const Real maximum = spot->value() *
1622 qTS->discount(d: firstDivDate)/rTS->discount(d: firstDivDate);
1623
1624 const Real minimum = (1 - firstDivAmount
1625 /(spot->value()*qTS->discount(d: firstDivDate)/rTS->discount(d: firstDivDate)))
1626 * spot->value()*qTS->discount(d: secondDivDate)/rTS->discount(d: secondDivDate)
1627 - secondDivAmount;
1628
1629 const Real calculatedMax = std::exp(x: loc.back());
1630 const Real calculatedMin = std::exp(x: loc.front());
1631
1632
1633 constexpr double relTol = 1e5*QL_EPSILON;
1634
1635 const Real maxDiff = std::fabs(x: calculatedMax - maximum);
1636 if (maxDiff > relTol*maximum) {
1637 BOOST_FAIL("Upper bound for Black-Scholes mesher failed: "
1638 << "\n calculated: " << calculatedMax
1639 << "\n expected: " << maximum
1640 << "\n difference: " << maxDiff
1641 << "\n tolerance: " << relTol*maximum);
1642 }
1643
1644 const Real minDiff = std::fabs(x: calculatedMin - minimum);
1645 if (minDiff > relTol*minimum) {
1646 BOOST_FAIL("Lower bound for Black-Scholes mesher failed: "
1647 << "\n calculated: " << calculatedMin
1648 << "\n expected: " << minimum
1649 << "\n difference: " << minDiff
1650 << "\n tolerance: " << relTol*minimum);
1651 }
1652}
1653
1654test_suite* FdmLinearOpTest::suite(SpeedLevel speed) {
1655 auto* suite = BOOST_TEST_SUITE("linear operator tests");
1656
1657 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFdmLinearOpLayout));
1658 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testUniformGridMesher));
1659 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFirstDerivativesMapApply));
1660 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testSecondDerivativesMapApply));
1661 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testDerivativeWeightsOnNonUniformGrids));
1662 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testSecondOrderMixedDerivativesMapApply));
1663 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testTripleBandMapSolve));
1664 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFdmHestonBarrier));
1665 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFdmHestonAmerican));
1666 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFdmHestonExpress));
1667 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testBiCGstab));
1668 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testGMRES));
1669 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testCrankNicolsonWithDamping));
1670 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testSpareMatrixReference));
1671 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testSparseMatrixZeroAssignment));
1672 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFdmMesherIntegral));
1673 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testHighInterestRateBlackScholesMesher));
1674 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testLowVolatilityHighDiscreteDividendBlackScholesMesher));
1675
1676 if (speed <= Fast) {
1677 suite->add(QUANTLIB_TEST_CASE(&FdmLinearOpTest::testFdmHestonHullWhiteOp));
1678 }
1679
1680 return suite;
1681}
1682

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