/*
Copyright (C) 2009 Ueli Hofstetter
This source code is release under the BSD License.
This file is part of JQuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://jquantlib.org/
JQuantLib is free software: you can redistribute it and/or modify it
under the terms of the JQuantLib license. You should have received a
copy of the license along with this program; if not, please email
. The license is also available online at
.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
JQuantLib is based on QuantLib. http://quantlib.org/
When applicable, the original copyright notice follows this notice.
*/
/*
Copyright (C) 2003, 2004, 2007 Ferdinando Ametrano
Copyright (C) 2006 Yiping Chen
Copyright (C) 2007 Neil Firth
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license. You should have received a
copy of the license along with this program; if not, please email
. The license is also available online at
.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the license for more details.
*/
package org.jquantlib.math.matrixutilities;
import org.jquantlib.QL;
import org.jquantlib.lang.exceptions.LibraryException;
import org.jquantlib.math.Closeness;
// TODO: code review :: please verify against QL/C++ code
/**
* @author Ueli Hofstetter
*/
public class PseudoSqrt {
private final static String unknown_salvaging_algorithm = "unknown salvaging algorithm";
public enum SalvagingAlgorithm {
None, Spectral, Hypersphere, LowerDiagonal, Higham;
}
//! Returns the pseudo square root of a real symmetric matrix
/*! Given a matrix \f$ M \f$, the result \f$ S \f$ is defined
as the matrix such that \f$ S S^T = M. \f$
If the matrix is not positive semi definite, it can
return an approximation of the pseudo square root
using a (user selected) salvaging algorithm.
For more information see: "The most general methodology to create
a valid correlation matrix for risk management and option pricing
purposes", by R. Rebonato and P. Jaeckel.
The Journal of Risk, 2(2), Winter 1999/2000
http://www.rebonato.com/correlationmatrix.pdf
Revised and extended in "Monte Carlo Methods in Finance",
by Peter Jaeckel, Chapter 6.
\pre the given matrix must be symmetric.
\relates Matrix
\warning Higham algorithm only works for correlation matrices.
\test
- the correctness of the results is tested by reproducing
known good data.
- the correctness of the results is tested by checking
returned values against numerical calculations.
*/
public Matrix pseudoSqrt( final Matrix matrix){
return pseudoSqrt(matrix, SalvagingAlgorithm.None);
}
//! Returns the rank-reduced pseudo square root of a real symmetric matrix
/*! The result matrix has rank<=maxRank. If maxRank>=size, then the
specified percentage of eigenvalues out of the eigenvalues' sum is
retained.
If the input matrix is not positive semi definite, it can return an
approximation of the pseudo square root using a (user selected)
salvaging algorithm.
\pre the given matrix must be symmetric.
\relates Matrix
*/
public static Matrix rankReducedSqrt(final Matrix matrix,
final int maxRank,
final int componentRetainedPercentage,
final SalvagingAlgorithm sa){
QL.validateExperimentalMode();
QL.require(matrix.rows == matrix.columns(), Cells.MATRIX_MUST_BE_SQUARE); // QA:[RG]::verified
QL.require(checkSymmetry(matrix), Cells.MATRIX_MUST_BE_SYMMETRIC); // QA:[RG]::verified
QL.require(componentRetainedPercentage>0.0, "no eigenvalues retained"); // TODO: message
QL.require(componentRetainedPercentage<=1.0, "percentage to be retained > 100%"); // TODO: message
QL.require(maxRank>=1, "max rank required < 1"); // TODO: message
final int size = matrix.rows;
// spectral (a.k.a Principal Component) analysis
SymmetricSchurDecomposition jd = new SymmetricSchurDecomposition(matrix);
Array eigenValues = jd.eigenvalues();
// salvaging algorithm
switch (sa) {
case None:
// eigenvalues are sorted in decreasing order
if(eigenValues.get(size-1)<-1e-16)
throw new IllegalArgumentException("negative eigenvalue(s) ("
+ eigenValues.get(size-1)+")");
break;
case Spectral:
// negative eigenvalues set to zero
for (int i=0; i0.0) {
final double normAdj = Math.sqrt(matrix.get(i,i)/norm);
for (int j=0; j1-eps)
theta[i*(i-1)/2+j]=1-eps;
if (theta[i*(i-1)/2+j]<-1+eps)
theta[i*(i-1)/2+j]=-1+eps;
for (k=0; k1-eps)
theta[i*(i-1)/2+j]=1-eps;
if (theta[i*(i-1)/2+j]<-1+eps)
theta[i*(i-1)/2+j]=-1+eps;
}
theta[i*(i-1)/2+j] = std::acos(theta[i*(i-1)/2+j]);
if (j==i-1) {
if (result[i][i]<0)
theta[i*(i-1)/2+j]=-theta[i*(i-1)/2+j];
}
}
}
Problem p(costFunction, constraint, theta);
optimize.minimize(p, endCriteria);
theta = p.currentValue();
std::fill(result.begin(),result.end(),1.0);
for (i=0; ii) {
result[i][k]=0;
} else {
for (j=0; j<=k; j++) {
if (j == k && k!=i)
result[i][k] *=
std::cos(theta[i*(i-1)/2+j]);
else if (j!=i)
result[i][k] *=
std::sin(theta[i*(i-1)/2+j]);
}
}
}
}
} else {
Array theta(size * (size-1));
const Real eps=1e-16;
for (i=0; i1-eps)
theta[j*size+i]=1-eps;
if (theta[j*size+i]<-1+eps)
theta[j*size+i]=-1+eps;
for (k=0;k1-eps)
theta[j*size+i]=1-eps;
if (theta[j*size+i]<-1+eps)
theta[j*size+i]=-1+eps;
}
theta[j*size+i] = std::acos(theta[j*size+i]);
if (j==size-2) {
if (result[i][j+1]<0)
theta[j*size+i]=-theta[j*size+i];
}
}
}
Problem p(costFunction, constraint, theta);
optimize.minimize(p, endCriteria);
theta=p.currentValue();
std::fill(result.begin(),result.end(),1.0);
for (i=0; i
Real normInf(const Matrix& M) {
Size rows = M.rows();
Size cols = M.columns();
Real norm = 0.0;
for (Size i=0; i
projectToUnitDiagonalMatrix(const Matrix& M) {
Size size = M.rows();
QL_REQUIRE(size == M.columns(),
"matrix not square");
Matrix result(M);
for (Size i=0; i
projectToPositiveSemidefiniteMatrix(Matrix& M) {
Size size = M.rows();
QL_REQUIRE(size == M.columns(),
"matrix not square");
Matrix diagonal(size, size, 0.0);
SymmetricSchurDecomposition jd(M);
for (Size i=0; i(jd.eigenvalues()[i], 0.0);
Matrix result =
jd.eigenvectors()*diagonal*transpose(jd.eigenvectors());
return result;
}
// implementation of the Higham algorithm to find the nearest
// correlation matrix.
const Disposable
highamImplementation(const Matrix& A,
const Size maxIterations,
const Real& tolerance) {
Size size = A.rows();
Matrix R, Y(A), X(A), deltaS(size, size, 0.0);
Matrix lastX(X);
Matrix lastY(Y);
for (Size i=0; i rankReducedSqrt(const Matrix& matrix,
Size maxRank,
Real componentRetainedPercentage,
SalvagingAlgorithm::Type sa) {
Size size = matrix.rows();
#if defined(QL_EXTRA_SAFETY_CHECKS)
checkSymmetry(matrix);
#else
QL_REQUIRE(size == matrix.columns(),
"non square matrix: " << size << " rows, " <<
matrix.columns() << " columns");
#endif
QL_REQUIRE(componentRetainedPercentage>0.0,
"no eigenvalues retained");
QL_REQUIRE(componentRetainedPercentage<=1.0,
"percentage to be retained > 100%");
QL_REQUIRE(maxRank>=1,
"max rank required < 1");
// spectral (a.k.a Principal Component) analysis
SymmetricSchurDecomposition jd(matrix);
Array eigenValues = jd.eigenvalues();
// salvaging algorithm
switch (sa) {
case SalvagingAlgorithm::None:
// eigenvalues are sorted in decreasing order
QL_REQUIRE(eigenValues[size-1]>=-1e-16,
"negative eigenvalue(s) ("
<< std::scientific << eigenValues[size-1]
<< ")");
break;
case SalvagingAlgorithm::Spectral:
// negative eigenvalues set to zero
for (Size i=0; i(eigenValues[i], 0.0);
break;
case SalvagingAlgorithm::Higham:
{
int maxIterations = 40;
Real tolerance = 1e-6;
Matrix adjustedMatrix = highamImplementation(matrix, maxIterations, tolerance);
jd = SymmetricSchurDecomposition(adjustedMatrix);
eigenValues = jd.eigenvalues();
}
break;
default:
QL_FAIL("unknown or invalid salvaging algorithm");
}
// factor reduction
Real enough = componentRetainedPercentage *
std::accumulate(eigenValues.begin(),
eigenValues.end(), 0.0);
if (componentRetainedPercentage == 1.0) {
// numerical glitches might cause some factors to be discarded
enough *= 1.1;
}
// retain at least one factor
Real components = eigenValues[0];
Size retainedFactors = 1;
for (Size i=1; components
Ads help maintain this website.