QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
factorreduction.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2009 Jose Aparicio
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
22#include <vector>
23
24namespace QuantLib {
25
26 std::vector<Real> factorReduction(Matrix mtrx,
27 Size maxIters) {
28 static Real tolerance = 1.e-6;
29
30 QL_REQUIRE(mtrx.rows() == mtrx.columns(),
31 "Input matrix is not square");
32
33 const Size n = mtrx.columns();
34
35 #if defined(QL_EXTRA_SAFETY_CHECKS)
36 // check symmetry
37 for (Size iRow=0; iRow<mtrx.rows(); iRow++)
38 for (Size iCol=0; iCol<iRow; iCol++)
39 QL_REQUIRE(mtrx[iRow][iCol] == mtrx[iCol][iRow],
40 "input matrix is not symmetric");
41 QL_REQUIRE(*std::max_element(mtrx.begin(), mtrx.end()) <=1 &&
42 *std::min_element(mtrx.begin(), mtrx.end()) >=-1,
43 "input matrix data is not correlation data");
44 #endif
45
46 // Initial guess value
47 std::vector<Real> previousCorrels(n, 0.);
48 for(Size iCol=0; iCol<n; iCol++) {
49 for(Size iRow=0; iRow<n; iRow++)
50 previousCorrels[iCol] +=
51 mtrx[iRow][iCol] * mtrx[iRow][iCol];
52 previousCorrels[iCol] =
53 std::sqrt((previousCorrels[iCol]-1.)/(n-1.));
54 }
55
56 // iterative solution
57 Size iteration = 0;
58 Real distance;
59 do {
60 // patch Matrix diagonal
61 for(Size iCol=0; iCol<n; iCol++)
62 mtrx[iCol][iCol] =
63 previousCorrels[iCol];
64 // compute eigenvector decomposition
66 //const Matrix& eigenVect = ssDec.eigenvectors();
67 const Array& eigenVals = ssDec.eigenvalues();
68 Array::const_iterator itMaxEval =
69 std::max_element(eigenVals.begin(), eigenVals.end());
70 // We do not need the max value, only the position of the
71 // corresponding eigenvector
72 Size iMax = std::distance(eigenVals.begin(), itMaxEval);
73 std::vector<Real> newCorrels, distances;
74 for(Size iCol=0; iCol<n; iCol++) {
75 Real thisCorrel = mtrx[iMax][iCol];
76 newCorrels.push_back(thisCorrel);
77 // strictly is:
78 // abs(\sqrt{\rho}- \sqrt{\rho_{old}})/\sqrt{\rho_{old}}
79 distances.push_back(
80 std::abs(thisCorrel - previousCorrels[iCol])/
81 previousCorrels[iCol]);
82 }
83 previousCorrels = newCorrels;
84 distance = *std::max_element(distances.begin(), distances.end());
85 }while(distance > tolerance && ++iteration <= maxIters );
86
87 // test it did not go up to the max iteration and the matrix can
88 // be reduced to one factor.
89 QL_REQUIRE(iteration < maxIters,
90 "convergence not reached after " <<
91 iteration << " iterations");
92
93 #if defined(QL_EXTRA_SAFETY_CHECKS)
94 QL_REQUIRE(*std::max_element(previousCorrels.begin(),
95 previousCorrels.end()) <=1 &&
96 *std::min_element(previousCorrels.begin(),
97 previousCorrels.end()) >=-1,
98 "matrix can not be decomposed to a single factor dependence");
99 #endif
100
101 return previousCorrels;
102 }
103
104}
1-D array used in linear algebra.
Definition: array.hpp:52
const Real * const_iterator
Definition: array.hpp:124
const_iterator end() const
Definition: array.hpp:511
const_iterator begin() const
Definition: array.hpp:503
Matrix used in linear algebra.
Definition: matrix.hpp:41
const_iterator begin() const
Definition: matrix.hpp:327
Size rows() const
Definition: matrix.hpp:504
const_iterator end() const
Definition: matrix.hpp:335
Size columns() const
Definition: matrix.hpp:508
symmetric threshold Jacobi algorithm.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Single factor correlation reduction.
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
std::vector< Real > factorReduction(Matrix mtrx, Size maxIters)
Eigenvalues/eigenvectors of a real symmetric matrix.