[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) 2003, 2004 Ferdinando Ametrano
5 Copyright (C) 2003 RiskMap srl
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 "covariance.hpp"
22#include "utilities.hpp"
23#include <ql/math/matrixutilities/getcovariance.hpp>
24#include <ql/math/matrixutilities/pseudosqrt.hpp>
25#include <ql/math/statistics/sequencestatistics.hpp>
26
27using namespace QuantLib;
28using namespace boost::unit_test_framework;
29
30namespace covariance_test {
31
32 Real norm(const Matrix& m) {
33 Real sum = 0.0;
34 for (Size i=0; i<m.rows(); i++)
35 for (Size j=0; j<m.columns(); j++)
36 sum += m[i][j]*m[i][j];
37 return std::sqrt(x: sum);
38 }
39
40}
41
42
43void CovarianceTest::testRankReduction() {
44
45 BOOST_TEST_MESSAGE("Testing matrix rank reduction salvaging algorithms...");
46
47 using namespace covariance_test;
48
49 Real expected, calculated;
50
51 Size n = 3;
52
53 Matrix badCorr(n, n);
54 badCorr[0][0] = 1.0; badCorr[0][1] = 0.9; badCorr[0][2] = 0.7;
55 badCorr[1][0] = 0.9; badCorr[1][1] = 1.0; badCorr[1][2] = 0.3;
56 badCorr[2][0] = 0.7; badCorr[2][1] = 0.3; badCorr[2][2] = 1.0;
57
58 Matrix goodCorr(n, n);
59 goodCorr[0][0] = goodCorr[1][1] = goodCorr[2][2] = 1.00000000000;
60 goodCorr[0][1] = goodCorr[1][0] = 0.894024408508599;
61 goodCorr[0][2] = goodCorr[2][0] = 0.696319066114392;
62 goodCorr[1][2] = goodCorr[2][1] = 0.300969036104592;
63
64 Matrix b = rankReducedSqrt(badCorr, maxRank: 3, componentRetainedPercentage: 1.0, SalvagingAlgorithm::Spectral);
65 Matrix calcCorr = b * transpose(m: b);
66
67 for (Size i=0; i<n; i++) {
68 for (Size j=0; j<n; j++) {
69 expected = goodCorr[i][j];
70 calculated = calcCorr[i][j];
71 if (std::fabs(x: calculated-expected) > 1.0e-10)
72 BOOST_ERROR("Salvaging correlation with spectral alg "
73 "through rankReducedSqrt "
74 << "cor[" << i << "][" << j << "]:\n"
75 << std::setprecision(10)
76 << " calculated: " << calculated << "\n"
77 << " expected: " << expected);
78 }
79 }
80
81 Matrix badCov(n, n);
82 badCov[0][0] = 0.04000; badCov[0][1] = 0.03240; badCov[0][2] = 0.02240;
83 badCov[1][0] = 0.03240; badCov[1][1] = 0.03240; badCov[1][2] = 0.00864;
84 badCov[2][0] = 0.02240; badCov[2][1] = 0.00864; badCov[2][2] = 0.02560;
85
86 b = pseudoSqrt(badCov, SalvagingAlgorithm::Spectral);
87 b = rankReducedSqrt(badCov, maxRank: 3, componentRetainedPercentage: 1.0, SalvagingAlgorithm::Spectral);
88 Matrix goodCov = b * transpose(m: b);
89
90 Real error = norm(m: goodCov-badCov);
91 if (error > 4.0e-4)
92 BOOST_ERROR(
93 std::scientific << error
94 << " error while salvaging covariance matrix with spectral alg "
95 "through rankReducedSqrt\n"
96 << std::fixed
97 << "input matrix:\n" << badCov
98 << "salvaged matrix:\n" << goodCov);
99}
100
101void CovarianceTest::testSalvagingMatrix() {
102
103 BOOST_TEST_MESSAGE("Testing positive semi-definiteness salvaging "
104 "algorithms...");
105
106 using namespace covariance_test;
107
108 Real expected, calculated;
109
110 Size n = 3;
111
112 Matrix badCorr(n, n);
113 badCorr[0][0] = 1.0; badCorr[0][1] = 0.9; badCorr[0][2] = 0.7;
114 badCorr[1][0] = 0.9; badCorr[1][1] = 1.0; badCorr[1][2] = 0.3;
115 badCorr[2][0] = 0.7; badCorr[2][1] = 0.3; badCorr[2][2] = 1.0;
116
117 Matrix goodCorr(n, n);
118 goodCorr[0][0] = goodCorr[1][1] = goodCorr[2][2] = 1.00000000000;
119 goodCorr[0][1] = goodCorr[1][0] = 0.894024408508599;
120 goodCorr[0][2] = goodCorr[2][0] = 0.696319066114392;
121 goodCorr[1][2] = goodCorr[2][1] = 0.300969036104592;
122
123 Matrix b = pseudoSqrt(badCorr, SalvagingAlgorithm::Spectral);
124// Matrix b = pseudoSqrt(badCorr, Hypersphere);
125 Matrix calcCorr = b * transpose(m: b);
126
127 for (Size i=0; i<n; i++) {
128 for (Size j=0; j<n; j++) {
129 expected = goodCorr[i][j];
130 calculated = calcCorr[i][j];
131 if (std::fabs(x: calculated-expected) > 1.0e-10)
132 BOOST_ERROR("SalvagingCorrelation with spectral alg "
133 << "cor[" << i << "][" << j << "]:\n"
134 << std::setprecision(10)
135 << " calculated: " << calculated << "\n"
136 << " expected: " << expected);
137 }
138 }
139
140 Matrix badCov(n, n);
141 badCov[0][0] = 0.04000; badCov[0][1] = 0.03240; badCov[0][2] = 0.02240;
142 badCov[1][0] = 0.03240; badCov[1][1] = 0.03240; badCov[1][2] = 0.00864;
143 badCov[2][0] = 0.02240; badCov[2][1] = 0.00864; badCov[2][2] = 0.02560;
144
145 b = pseudoSqrt(badCov, SalvagingAlgorithm::Spectral);
146 Matrix goodCov = b * transpose(m: b);
147
148 Real error = norm(m: goodCov-badCov);
149 if (error > 4.0e-4)
150 BOOST_ERROR(
151 std::scientific << error
152 << " error while salvaging covariance matrix with spectral alg\n"
153 << std::fixed
154 << "input matrix:\n" << badCov
155 << "salvaged matrix:\n" << goodCov);
156}
157
158void CovarianceTest::testCovariance() {
159
160 BOOST_TEST_MESSAGE("Testing covariance and correlation calculations...");
161
162 std::vector<std::vector<Real>> data = {
163 { 3.0, 9.0 },
164 { 2.0, 7.0 },
165 { 4.0, 12.0 },
166 { 5.0, 15.0 },
167 { 6.0, 17.0 }
168 };
169 std::vector<Real> weights(data.size(), 1.0);
170
171 Size i, j, n = data[0].size();
172
173 Matrix expCor(n, n);
174 expCor[0][0] = 1.0000000000000000; expCor[0][1] = 0.9970544855015813;
175 expCor[1][0] = 0.9970544855015813; expCor[1][1] = 1.0000000000000000;
176
177 SequenceStatistics s(n);
178 std::vector<Real> temp(n);
179
180 for (i = 0; i<data.size(); i++) {
181 for (j=0; j<n; j++) {
182 temp[j]= data[i][j];
183 }
184 s.add(sample: temp, weight: weights[i]);
185 }
186
187 std::vector<Real> std = s.standardDeviation();
188 Matrix calcCov = s.covariance();
189 Matrix calcCor = s.correlation();
190
191 Matrix expCov(n, n);
192 for (i=0; i<n; i++) {
193 expCov[i][i] = std[i]*std[i];
194 for (j=0; j<i; j++) {
195 expCov[i][j] = expCov[j][i] = expCor[i][j]*std[i]*std[j];
196 }
197 }
198
199 Real expected, calculated;
200 for (i=0; i<n; i++) {
201 for (j=0; j<n; j++) {
202 expected = expCor[i][j];
203 calculated = calcCor[i][j];
204 if (std::fabs(x: calculated-expected) > 1.0e-10)
205 BOOST_ERROR("SequenceStatistics "
206 << "cor[" << i << "][" << j << "]:\n"
207 << std::setprecision(10)
208 << " calculated: " << calculated << "\n"
209 << " expected: " << expected);
210
211 expected = expCov[i][j];
212 calculated = calcCov[i][j];
213 if (std::fabs(x: calculated-expected) > 1.0e-10)
214 BOOST_ERROR("SequenceStatistics "
215 << "cov[" << i << "][" << j << "]:\n"
216 << std::setprecision(10)
217 << " calculated: " << calculated << "\n"
218 << " expected: " << expected);
219 }
220 }
221
222 calcCov = getCovariance(stdDevBegin: std.begin(), stdDevEnd: std.end(), corr: expCor);
223
224 for (i=0; i<n; i++) {
225 for (j=0; j<n; j++) {
226 Real calculated = calcCov[i][j],
227 expected = expCov[i][j];
228 if (std::fabs(x: calculated-expected) > 1.0e-10) {
229 BOOST_ERROR("getCovariance "
230 << "cov[" << i << "][" << j << "]:\n"
231 << std::setprecision(10)
232 << " calculated: " << calculated << "\n"
233 << " expected: " << expected);
234 }
235 }
236 }
237
238
239
240
241 CovarianceDecomposition covDecomposition(expCov);
242 calcCor = covDecomposition.correlationMatrix();
243 Array calcStd = covDecomposition.standardDeviations();
244
245 for (i=0; i<n; i++) {
246 calculated = calcStd[i];
247 expected = std[i];
248 if (std::fabs(x: calculated-expected) > 1.0e-16) {
249 BOOST_ERROR("CovarianceDecomposition "
250 << "standardDev[" << i << "]:\n"
251 << std::setprecision(16) << std::scientific
252 << " calculated: " << calculated << "\n"
253 << " expected: " << expected);
254 }
255 for (j=0; j<n; j++) {
256 calculated = calcCor[i][j];
257 expected = expCor[i][j];
258 if (std::fabs(x: calculated-expected) > 1.0e-14) {
259 BOOST_ERROR("\nCovarianceDecomposition "
260 << "corr[" << i << "][" << j << "]:\n"
261 << std::setprecision(14) << std::scientific
262 << " calculated: " << calculated << "\n"
263 << " expected: " << expected);
264 }
265 }
266 }
267
268
269
270}
271
272
273test_suite* CovarianceTest::suite() {
274 auto* suite = BOOST_TEST_SUITE("Covariance and correlation tests");
275 suite->add(QUANTLIB_TEST_CASE(&CovarianceTest::testCovariance));
276 suite->add(QUANTLIB_TEST_CASE(&CovarianceTest::testSalvagingMatrix));
277 suite->add(QUANTLIB_TEST_CASE(&CovarianceTest::testRankReduction));
278 return suite;
279}
280
281

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