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
27 using namespace QuantLib;
28 using namespace boost::unit_test_framework;
29
30 namespace covariance_test {
31
norm(const Matrix & m)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(sum);
38 }
39
40 }
41
42
testRankReduction()43 void 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, 3, 1.0, SalvagingAlgorithm::Spectral);
65 Matrix calcCorr = b * transpose(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(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, 3, 1.0, SalvagingAlgorithm::Spectral);
88 Matrix goodCov = b * transpose(b);
89
90 Real error = norm(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
testSalvagingMatrix()101 void 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(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(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(b);
147
148 Real error = norm(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
testCovariance()158 void CovarianceTest::testCovariance() {
159
160 BOOST_TEST_MESSAGE("Testing covariance and correlation calculations...");
161
162 Real data00[] = { 3.0, 9.0 };
163 Real data01[] = { 2.0, 7.0 };
164 Real data02[] = { 4.0, 12.0 };
165 Real data03[] = { 5.0, 15.0 };
166 Real data04[] = { 6.0, 17.0 };
167 Real* data[5] = { data00, data01, data02, data03, data04 };
168 std::vector<Real> weights(LENGTH(data), 1.0);
169
170 Size i, j, n = LENGTH(data00);
171
172 Matrix expCor(n, n);
173 expCor[0][0] = 1.0000000000000000; expCor[0][1] = 0.9970544855015813;
174 expCor[1][0] = 0.9970544855015813; expCor[1][1] = 1.0000000000000000;
175
176 SequenceStatistics s(n);
177 std::vector<Real> temp(n);
178
179 for (i = 0; i<LENGTH(data); i++) {
180 for (j=0; j<n; j++) {
181 temp[j]= data[i][j];
182 }
183 s.add(temp, weights[i]);
184 }
185
186 std::vector<Real> std = s.standardDeviation();
187 Matrix calcCov = s.covariance();
188 Matrix calcCor = s.correlation();
189
190 Matrix expCov(n, n);
191 for (i=0; i<n; i++) {
192 expCov[i][i] = std[i]*std[i];
193 for (j=0; j<i; j++) {
194 expCov[i][j] = expCov[j][i] = expCor[i][j]*std[i]*std[j];
195 }
196 }
197
198 Real expected, calculated;
199 for (i=0; i<n; i++) {
200 for (j=0; j<n; j++) {
201 expected = expCor[i][j];
202 calculated = calcCor[i][j];
203 if (std::fabs(calculated-expected) > 1.0e-10)
204 BOOST_ERROR("SequenceStatistics "
205 << "cor[" << i << "][" << j << "]:\n"
206 << std::setprecision(10)
207 << " calculated: " << calculated << "\n"
208 << " expected: " << expected);
209
210 expected = expCov[i][j];
211 calculated = calcCov[i][j];
212 if (std::fabs(calculated-expected) > 1.0e-10)
213 BOOST_ERROR("SequenceStatistics "
214 << "cov[" << i << "][" << j << "]:\n"
215 << std::setprecision(10)
216 << " calculated: " << calculated << "\n"
217 << " expected: " << expected);
218 }
219 }
220
221 calcCov = getCovariance(std.begin(), std.end(), expCor);
222
223 for (i=0; i<n; i++) {
224 for (j=0; j<n; j++) {
225 Real calculated = calcCov[i][j],
226 expected = expCov[i][j];
227 if (std::fabs(calculated-expected) > 1.0e-10) {
228 BOOST_ERROR("getCovariance "
229 << "cov[" << i << "][" << j << "]:\n"
230 << std::setprecision(10)
231 << " calculated: " << calculated << "\n"
232 << " expected: " << expected);
233 }
234 }
235 }
236
237
238
239
240 CovarianceDecomposition covDecomposition(expCov);
241 calcCor = covDecomposition.correlationMatrix();
242 Array calcStd = covDecomposition.standardDeviations();
243
244 for (i=0; i<n; i++) {
245 calculated = calcStd[i];
246 expected = std[i];
247 if (std::fabs(calculated-expected) > 1.0e-16) {
248 BOOST_ERROR("CovarianceDecomposition "
249 << "standardDev[" << i << "]:\n"
250 << std::setprecision(16) << std::scientific
251 << " calculated: " << calculated << "\n"
252 << " expected: " << expected);
253 }
254 for (j=0; j<n; j++) {
255 calculated = calcCor[i][j];
256 expected = expCor[i][j];
257 if (std::fabs(calculated-expected) > 1.0e-14) {
258 BOOST_ERROR("\nCovarianceDecomposition "
259 << "corr[" << i << "][" << j << "]:\n"
260 << std::setprecision(14) << std::scientific
261 << " calculated: " << calculated << "\n"
262 << " expected: " << expected);
263 }
264 }
265 }
266
267
268
269 }
270
271
suite()272 test_suite* CovarianceTest::suite() {
273 test_suite* suite = BOOST_TEST_SUITE("Covariance and correlation tests");
274 suite->add(QUANTLIB_TEST_CASE(&CovarianceTest::testCovariance));
275 suite->add(QUANTLIB_TEST_CASE(&CovarianceTest::testSalvagingMatrix));
276 suite->add(QUANTLIB_TEST_CASE(&CovarianceTest::testRankReduction));
277 return suite;
278 }
279
280