1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 RiskMap srl
5 Copyright (C) 2011 Ferdinando Ametrano
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include "operators.hpp"
22#include "utilities.hpp"
23#include <ql/math/distributions/normaldistribution.hpp>
24#include <ql/methods/finitedifferences/dzero.hpp>
25#include <ql/methods/finitedifferences/dplusdminus.hpp>
26#include <ql/methods/finitedifferences/bsmoperator.hpp>
27#include <ql/methods/finitedifferences/bsmtermoperator.hpp>
28#include <ql/time/daycounters/actual360.hpp>
29#include <ql/quotes/simplequote.hpp>
30#include <ql/utilities/dataformatters.hpp>
31
32using namespace QuantLib;
33using namespace boost::unit_test_framework;
34
35
36void OperatorTest::testTridiagonal() {
37
38 BOOST_TEST_MESSAGE("Testing tridiagonal operator...");
39
40 Size n = 8; // can use 3 for easier debugging
41
42 TridiagonalOperator T(n);
43 T.setFirstRow(valB: 1.0, valC: 2.0);
44 T.setMidRows( valA: 0.0, valB: 2.0, valC: 0.0);
45 T.setLastRow( valA: 2.0, valB: 1.0);
46
47 Array original(n, 1.0);
48
49 Array intermediate = T.applyTo(v: original);
50
51 Array final(intermediate);
52 T.solveFor(rhs: final, result&: final);
53 for (Size i=0; i<n; ++i) {
54 if (final[i]!=original[i])
55 BOOST_FAIL("\n applyTo + solveFor does not equal identity:"
56 "\n original vector: " << original <<
57 "\n transformed vector: " << intermediate <<
58 "\n inverse transformed vector: " << final);
59 }
60
61 final = Array(n, 0.0);
62 Array temp(intermediate);
63 T.solveFor(rhs: temp, result&: final);
64 for (Size i=0; i<n; ++i) {
65 if (temp[i]!=intermediate[i])
66 BOOST_FAIL("\n solveFor altered rhs:"
67 "\n original vector: " << original <<
68 "\n transformed vector: " << intermediate <<
69 "\n altered transformed vector: " << temp <<
70 "\n inverse transformed vector: " << final);
71 }
72 for (Size i=0; i<n; ++i) {
73 if (final[i]!=original[i])
74 BOOST_FAIL("\n applyTo + solveFor does not equal identity:"
75 "\n original vector: " << original <<
76 "\n transformed vector: " << intermediate <<
77 "\n inverse transformed vector: " << final);
78 }
79
80 final = T.solveFor(rhs: temp);
81 for (Size i=0; i<n; ++i) {
82 if (temp[i]!=intermediate[i])
83 BOOST_FAIL("\n solveFor altered rhs:"
84 "\n original vector: " << original <<
85 "\n transformed vector: " << intermediate <<
86 "\n altered transformed vector: " << temp <<
87 "\n inverse transformed vector: " << final);
88 }
89 for (Size i=0; i<n; ++i) {
90 if (final[i]!=original[i])
91 BOOST_FAIL("\n applyTo + solveFor does not equal identity:"
92 "\n original vector: " << original <<
93 "\n transformed vector: " << intermediate <<
94 "\n inverse transformed vector: " << final);
95 }
96
97 Real delta, error = 0.0, tolerance = 1e-9;
98 final = T.SOR(rhs: temp, tol: tolerance);
99 for (Size i=0; i<n; ++i) {
100 delta = final[i]-original[i];
101 error += delta * delta;
102 if (temp[i]!=intermediate[i])
103 BOOST_FAIL("\n SOR altered rhs:"
104 "\n original vector: " << original <<
105 "\n transformed vector: " << intermediate <<
106 "\n altered transformed vector: " << temp <<
107 "\n inverse transformed vector: " << final);
108 }
109 if (error>tolerance)
110 BOOST_FAIL("\n applyTo + SOR does not equal identity:"
111 "\n original vector: " << original <<
112 "\n transformed vector: " << intermediate <<
113 "\n inverse transformed vector: " << final <<
114 "\n error: " << error <<
115 "\n tolerance: " << tolerance);
116}
117
118void OperatorTest::testConsistency() {
119
120 BOOST_TEST_MESSAGE("Testing differential operators...");
121
122 Real average = 0.0, sigma = 1.0;
123
124 NormalDistribution normal(average,sigma);
125 CumulativeNormalDistribution cum(average,sigma);
126
127 Real xMin = average - 4*sigma,
128 xMax = average + 4*sigma;
129 Size N = 10001;
130 Real h = (xMax-xMin)/(N-1);
131
132 Array x(N), y(N), yi(N), yd(N), temp(N), diff(N);
133
134 Size i;
135 for (i=0; i<N; i++)
136 x[i] = xMin+h*i;
137 std::transform(first: x.begin(),last: x.end(),result: y.begin(),unary_op: normal);
138 std::transform(first: x.begin(),last: x.end(),result: yi.begin(),unary_op: cum);
139 for (i=0; i<x.size(); i++)
140 yd[i] = normal.derivative(x: x[i]);
141
142 // define the differential operators
143 DZero D(N,h);
144 DPlusDMinus D2(N,h);
145
146 // check that the derivative of cum is Gaussian
147 temp = D.applyTo(v: yi);
148 std::transform(first1: y.begin(), last1: y.end(), first2: temp.begin(), result: diff.begin(), binary_op: std::minus<>());
149 Real e = norm(begin: diff.begin(), end: diff.end(), h);
150 if (e > 1.0e-6) {
151 BOOST_FAIL("norm of 1st derivative of cum minus Gaussian: " << e
152 << "\ntolerance exceeded");
153 }
154
155 // check that the second derivative of cum is normal.derivative
156 temp = D2.applyTo(v: yi);
157 std::transform(first1: yd.begin(), last1: yd.end(), first2: temp.begin(), result: diff.begin(), binary_op: std::minus<>());
158 e = norm(begin: diff.begin(), end: diff.end(), h);
159 if (e > 1.0e-4) {
160 BOOST_FAIL("norm of 2nd derivative of cum minus Gaussian derivative: "
161 << e << "\ntolerance exceeded");
162 }
163}
164
165void OperatorTest::testBSMOperatorConsistency() {
166 BOOST_TEST_MESSAGE("Testing consistency of BSM operators...");
167
168 Array grid(10);
169 Real price = 20.0;
170 Real factor = 1.1;
171 Size i;
172 for (i = 0; i < grid.size(); i++) {
173 grid[i] = price;
174 price *= factor;
175 }
176 Real dx = std::log(x: factor);
177 Rate r = 0.05;
178 Rate q = 0.01;
179 Volatility sigma = 0.5;
180
181 BSMOperator ref(grid.size(), dx, r, q, sigma);
182
183 DayCounter dc = Actual360();
184 Date today = Date::todaysDate();
185 Date exercise = today + 2*Years;
186 Time residualTime = dc.yearFraction(d1: today,d2: exercise);
187
188 ext::shared_ptr<SimpleQuote> spot(new SimpleQuote(0.0));
189 ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, forward: q, dc);
190 ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: r, dc);
191 ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, volatility: sigma, dc);
192 ext::shared_ptr<GeneralizedBlackScholesProcess> stochProcess(
193 new GeneralizedBlackScholesProcess(
194 Handle<Quote>(spot),
195 Handle<YieldTermStructure>(qTS),
196 Handle<YieldTermStructure>(rTS),
197 Handle<BlackVolTermStructure>(volTS)));
198 PdeOperator<PdeBSM> op2(grid, stochProcess, residualTime);
199
200 Real tolerance = 1.0e-6;
201
202 Array lderror = ref.lowerDiagonal() - op2.lowerDiagonal();
203 Array derror = ref.diagonal() - op2.diagonal();
204 Array uderror = ref.upperDiagonal() - op2.upperDiagonal();
205
206 for (i=2; i<grid.size()-2; i++) {
207 if (std::fabs(x: lderror[i]) > tolerance ||
208 std::fabs(x: derror[i]) > tolerance ||
209 std::fabs(x: uderror[i]) > tolerance) {
210 BOOST_FAIL("inconsistency between BSM operators:\n"
211 << io::ordinal(i) << " row:\n"
212 << "expected: "
213 << ref.lowerDiagonal()[i] << ", "
214 << ref.diagonal()[i] << ", "
215 << ref.upperDiagonal()[i] << "\n"
216 << "calculated: "
217 << op2.lowerDiagonal()[i] << ", "
218 << op2.diagonal()[i] << ", "
219 << op2.upperDiagonal()[i]);
220 }
221 }
222}
223
224
225test_suite* OperatorTest::suite() {
226 auto* suite = BOOST_TEST_SUITE("Operator tests");
227 suite->add(QUANTLIB_TEST_CASE(&OperatorTest::testTridiagonal));
228 suite->add(QUANTLIB_TEST_CASE(&OperatorTest::testConsistency));
229 suite->add(QUANTLIB_TEST_CASE(&OperatorTest::testBSMOperatorConsistency));
230 return suite;
231}
232
233

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