14
$\begingroup$

Does Mathematica have a way to "fix" a correlation matrix that is not positive semi-definite?

I looked through the documentation and search the internet but could not find anything.

$\endgroup$
14
$\begingroup$

Update: the previous method continues the iteration until all eigenvalues are at or above the threshold, that is, it does not check and stop if the current matrix is PD. The following is a version that tests for positive definiteness using PositiveDefiniteMatrixQ at every iteration:

ClearAll[covBending2];
covBending2[mat_, tol_: 1/10000] := If[PositiveDefiniteMatrixQ[mat], mat,
   NestWhile[ (Eigensystem[#][[2]].DiagonalMatrix[
     Max[#, tol] & /@ (Eigensystem[#][[1]])].Transpose[
     Eigensystem[#][[2]]]) &,  N@mat, !PositiveDefiniteMatrixQ[#]&]]

Original answer:

Here is simple (unweighted) Mma version of the Matlab implementation of Covariance Bending.

ClearAll[covBending];
covBending[mat_, tol_: 1/10000]:=If[PositiveDefiniteMatrixQ[mat], mat, 
NestWhile[(Eigensystem[#][[2]].DiagonalMatrix[
  Max[#, tol] & /@(Eigensystem[#][[1]])].Transpose[
  Eigensystem[#][[2]]]) &, N@mat,  Min[Eigensystem[#][[1]]] < tol &]]

examples:

 matrices = ((1/2) (# + Transpose[#]) & /@  RandomReal[{0, 1}, {10, 4, 4}]);
{MatrixForm[#], PositiveDefiniteMatrixQ[#], MatrixForm[covBending[#]],
 PositiveDefiniteMatrixQ[covBending[#]]} & /@  matrices // TableForm

enter image description here

$\endgroup$
  • $\begingroup$ Why not have the test be ! PositiveDefiniteMatrixQ[mat]? $\endgroup$ – J. M. will be back soon Aug 31 '12 at 7:59
  • $\begingroup$ @J.M. yes indeed, checking if mat is positive definite is actually much faster than checking if minimum eigenvalue of mat less than a threshold. $\endgroup$ – kglr Aug 31 '12 at 8:15
  • $\begingroup$ Does this ensure that the diagonal has only ones ? $\endgroup$ – faysou Aug 31 '12 at 12:20
  • $\begingroup$ @Faysal, no. This is for a generic square matrix. For a covariance matrix symmetry is preserved, I think, unless you use a an arbitrary weight vector/matrix. For correlation matrices you need to add renormalization. $\endgroup$ – kglr Aug 31 '12 at 12:29
  • $\begingroup$ @kguler Maybe I'm wrong, but isn't this a fundamentally different problem than "fixing" non psd correlation/covariance matrices? This looks like it belongs to the class of shrinkage estimators and is converting a general square symmetric matrix to a psd matrix. One could argue that if a matrix is demonstrably far from psd (in some sense) other than due to minor numerical errors, then it was never a covariance matrix to begin with. Of course, OP was very vague about what they meant by "fixing" and I do recognize that yours could've very well been their implied meaning... $\endgroup$ – rm -rf Aug 31 '12 at 21:04
11
$\begingroup$

Generally, the reason why matrices that were supposed to be positive semi-definite but are not, is because the constraint of working in a finite precision world often introduces a wee bit of perturbation in the lowest eigenvalues of the matrix, making it either negative or complex. These errors are generally of the order of machine precision, but is enough to return False if you use standard algorithms to check for positive semi-definiteness. This can be easily fixed with Chop. For example, something like the following:

psdMat = #1.Chop@#2.#3\[ConjugateTranspose] & @@ SingularValueDecomposition[mat];

Now if this isn't sufficient to make your matrix positive semi-definite, you should go back and take a closer look at your problem to see if there are other reasons to not expect it to be positive semi-definite.

$\endgroup$
  • 5
    $\begingroup$ The last sentence cannot be emphasized enough; in fact, before applying any "correction" method, you should try to figure out why things are not coming out to be positive (semi)definite in the first place. "Correction" is more of a Band-Aid than anything... $\endgroup$ – J. M. will be back soon Aug 31 '12 at 8:22
  • $\begingroup$ When calculating the correlation of multiple entities, there can be missing values between entity 1 and entity 3, but not between 2 and 3. If there are many of these, you can produce a correlation matrix that is not positive semi-definite. The easy answer is to just use the data where it exists for all 3 entities, but if you want to use all the available data this can be a problem. Another example of this, is that you may have some of the correlation matrix be known but other values are not and the user needs to enter what they think is the correct correlation amount. $\endgroup$ – Stephen Lien Aug 31 '12 at 13:53
8
$\begingroup$

I like the approach in this paper: The most general methodology to create a valid correlation matrix for risk management and option pricing purposes. Below is a starting point for an implementation.

MakeGoodCorrelationMatrix[badCorrMat_,flag_]:=
Which[
    flag==1,
        Module[{LocalN=Length[badCorrMat],LocalTheta,LocalB,LocalOutput,LocalSolution},
            LocalTheta=Outer[Unique[] &, Range[LocalN], Range[LocalN-1], 1, 1];
            LocalB = Outer[If[#2 < LocalN, Cos[LocalTheta[[#1, #2]]] (Times @@ Sin[LocalTheta[[#1, 1 ;; #2 - 1]]]),
                Times @@ Sin[LocalTheta[[#1, 1 ;; #2 - 1]]]] &, Range[LocalN], Range[LocalN], 1, 1] //Chop;
            LocalOutput = Dot[LocalB, Transpose[LocalB]] //Chop;
            LocalSolution = Minimize[Total[Flatten[(LocalOutput - badCorrMat)^2]], Flatten[LocalTheta]]//Chop;
            LocalOutput /. LocalSolution[[2]] 

            ],  
    True,
        Module[{LocalLambdap = (Max[0.000001, #] & /@ Eigenvalues[badCorrMat]) , LocalS = Eigenvectors[badCorrMat] // Transpose, LocalT, LocalB, LocalOutput}, 
            LocalT = DiagonalMatrix[Dot[#^2,LocalLambdap]^(-1) & /@ LocalS] ; 
            LocalB = Dot[Sqrt[LocalT],LocalS,Sqrt[DiagonalMatrix[LocalLambdap]]]; 

            LocalOutput = Dot[LocalB,Transpose[LocalB]] ;
            LocalOutput = 1/2 (LocalOutput + Transpose[LocalOutput]) (* just to make sure the output is symmetric *)

            ]
]

Example :

badMat = {{1, 0.05, 0.9}, {0.05, 1, 0.8}, {0.9, 0.8, 1}};
Eigenvalues[badMat]
(* {2.22926, 0.950346, -0.179603} *)

goodMat1 = MakeGoodCorrelationMatrix[badMat, 1];
Eigenvalues[goodMat1]
(* {2.11798, 0.882022, 8.32667*10^-17} *)

goodMat2 = MakeGoodCorrelationMatrix[badMat, 2];
Eigenvalues[goodMat2]
(* {2.08911, 0.910888, 9.35756*10^-7} *)
$\endgroup$
  • $\begingroup$ Your answer describes the code as "...a starting point for an implementation". Have you fully implemented the Rebonato Formula? Just starting to look at implementing it myself, translating from MATLAB code to Mathematica. Any thoughts or suggestions appreciated. $\endgroup$ – Jagra Mar 6 '18 at 1:36
  • $\begingroup$ @Jagra I think the above works but it's far from being "production grade". The first approach (flag==1) was quite slow for matrices above 5x5 (maybe use NMinimize ?), the second approach shouldn't have a hard-coded cutoff and I did not check it will give reasonable results for large matrices. Do post your own version when ready please. $\endgroup$ – b.gates.you.know.what Mar 6 '18 at 8:28
  • $\begingroup$ Will do. Very interesting area. Thx $\endgroup$ – Jagra Mar 6 '18 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.