I'm following Peter Jackel's book "Monte Carlo Methods in Finance", where an algorithm for transforming malformed correlation matrices into acceptable correlation matrices (positive semi-definite) to perform a Cholesky decomposition later on.
The basic approach is transforming the negative eigenvalues into 0 and then transforming the matrix to adjust to these new values. In http://www.mathworks.com/matlabcentral/fileexchange/26475-salvaging-a-linear-correlation-matrix/content/SpectralDecomp.m there's an implementation for matlab of this method.
The problem in this case is in floating point arithmetics. Suppose I have this matrix:
1 1 1 1 1 1 1 1 1
The eigenvalues for this matrix are:
-9.11271925974066E-17 1.11022302462516E-16 3
This matrix is not positive semi-definite, because of the first eigenvalue. However, when I try to adjust it to a new matrix that is positive semi-definite, I end up with the very same matrix as above! And the reason is quite clear: the first 2 eigenvalues are too close to zero already - if they're substituted by zero, nothing really changes at all. In fact, by inputting this matrix in a system with less precision (ie http://www.mathstools.com/section/main/system_equations_solver#.UYrN5rV-jK0), they do appear to be zero. In the final matrix multiplication, I end up with a sum that's not mathematically correct, but a consequence of floating point arithmetic:
3.57e-33 + 1 = 1
I'm using VBA for these calculations, with the main operations made mainly by the functions in the excellent Alglib library, but Excel only stores 15 precision digits for floating point numbers. Maybe I could try the hypershpere decomposition approach instead of the spectral decomposition, but I think I'll probably bump into the same issue. Do R, Matlab or any other mathematical software deal with this issue?
To reproduce the problem, try to multiply the following matrix by its own transpose:
0 7.45058059692383E-09 1 0 -7.45058059692383E-09 1 0 0 1
The result shouldn't be a matrix of 1's, but that's what's returned by the calculations.
Is there any way to circumvent this limitation? Or any other method to obtain a positive semi-definite matrix out of an arbitrary given correlation matrix? I need a Cholesky decomposition of this matrix to generate random samples for a copula, perhaps there's another method for this that does not involve the Cholesky decomposition and the consequent constraint of the matrix being positive semi-definite?