Let me preface this post by saying I am getting frustrated with R.  The syntax is not intuitive and the performance for matrix operations is slow.  Using Octave, a free Matlab clone, I can get over 6 Gflops on things that R is doing at less than 2.  After this post, I will focus on the statistical functions of R and will do any matrix operations in Octave, or maybe NumPy.


First, let’s look at a situation where you find your self with a non-positive semi-definite (PSD) matrix.  Remember that the default chol() requires a positive definite (PD) matrix and we created a function that can handle a PSD matrix.  Imagine if you have 3 time series with different sampling intervals — in our example we will have one that samples every day, another that samples M-W-F-Sa-Su, and one that samples T-R (Thursday) – Sa.  Further we only have a limited set of data – 3 weeks in our sample.  To create a correlation matrix, we would be forced to either only consider Saturdays (where all 3 series intersect), or to compute the correlations pairwise.  With only 3 intersection points, we would probably choose the pairwise method to attempt to get a better handle on the true underlying process.


We will use SAS IML to construct the series and PROC CORR from SAS to calculate the pairwise correlations.

proc iml;
/*True covariance of the system*/
C1 = {1 .5 0,
     .5 1 .25,
     0 .25 1};

r = root(C1)`;n=21;
out = J(n,3,.);

/*Generate 21 days worth of data*/

do i=1 to n;
   hold = J(1,3);
   do j=1 to 3;
      hold[j] = rannor(2);
   end;

   hold = r * hold`;
   out[i,] = hold`;

end; 

create out from out;
append from out;
close out;
quit;

/*Filter out the unobserved values*/
data out;
format date date9.;
set out;
retain date “04DEC2011″d;
if _n_ ^=1 then   date = date + 1;

if ^(weekday(date) in (2, 4, 6, 7, 1)) then   col2 = .;

if ^(weekday(date) in (3, 5, 7)) then   col3 = .;
run;

title “Observed Data”;
proc print data=out noobs;
run;

title “Pairwise Correlation, observed”;
proc corr data=out(drop=date)
          out=corr_out(where=(_type_ = “CORR”))
          ;run;

data corr_out;
set corr_out(drop=_type_ _name_);
run;

title “Correlation Eigenvalues”;
proc iml;
use corr_out;
read all into corr;
close corr_out; e = eigval(corr); print e;
quit;

Here is the output:

Observed Data

Pairwise Correlation, observed

Correlation Eigenvalues

The Matrix package from R has a function called nearPD().  This function has some nice features such as requiring the diagonal to remain constant (useful for a covariance matrix) and a flag to make sure the output is still a correlation matrix.  Not sure what the difference between those two would be if I used a correlation matrix as the input…  The output also include the new eigenvalues, and the Frobenius Norm (a measure of closeness for matrices).


Rebonato and Jackel (1999) define 2 methods in their paper for finding a PSD matrix.  The first is rather involved and the second, the Spectral Decomposition, is very compact.  The Spectral Decomposition is what we will emulate.  Take a second to go read that section (pp 7 – 9).  

Here it is in R:

nearPSD = function(c){
   n = dim(c)[1]
   e = eigen(c,sym=TRUE)
   val = e$values * (e$values > 0)
   vec = e$vectors
   T = sqrt(diag( as.vector(1/(vec^2 %*% val)),n,n))
   B = T %*% vec %*% diag(as.vector(sqrt(val)),n,n)
   out = B %*% t(B)
   return(out)
}

 Using the example from the paper, we see this duplicates their results:

> c
     [,1] [,2] [,3]
[1,]  1.0  0.9  0.7
[2,]  0.9  1.0  0.3
[3,]  0.7  0.3  1.0
> nearPSD(c)
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.8940244 0.6963191
[2,] 0.8940244 1.0000000 0.3009690
[3,] 0.6963191 0.3009690 1.0000000

n = 1000
c = matrix(.9,n,n)
for( i in 1:n){
  c[i,i] = 1
}
c[1,2] = 0
c[2,1] = 0
c[1:5,1:5]
     [,1] [,2] [,3] [,4] [,5]
[1,]  1.0  0.0  0.9  0.9  0.9
[2,]  0.0  1.0  0.9  0.9  0.9
[3,]  0.9  0.9  1.0  0.9  0.9
[4,]  0.9  0.9  0.9  1.0  0.9
[5,]  0.9  0.9  0.9  0.9  1.0

e = eigen(c,sym=TRUE)
min(e$values)
[1] -0.7982018

s = Sys.time()
pd = nearPD(c,corr=TRUE)
pd = as.matrix(pd$mat)
rpd = chol(pd)
e = Sys.time()
e – s
Time difference of 45.659 secs

pd[1:5,1:5]
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.7934317 0.8984058 0.8984058 0.8984058
[2,] 0.7934317 1.0000000 0.8984058 0.8984058 0.8984058
[3,] 0.8984058 0.8984058 1.0000000 0.9000032 0.9000032
[4,] 0.8984058 0.8984058 0.9000032 1.0000000 0.9000032
[5,] 0.8984058 0.8984058 0.9000032 0.9000032 1.0000000

s = Sys.time()
psd = nearPSD(c)
rpsd = my_chol(psd,50)
e = Sys.time()
e – s
Time difference of 7.031 secs

psd[1:5,1:5]
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.2848481 0.7604250 0.7604250 0.7604250
[2,] 0.2848481 1.0000000 0.7604250 0.7604250 0.7604250
[3,] 0.7604250 0.7604250 1.0000000 0.9000002 0.9000002
[4,] 0.7604250 0.7604250 0.9000002 1.0000000 0.9000002
[5,] 0.7604250 0.7604250 0.9000002 0.9000002 1.0000000

norm(c – pd ,type=”F”)
[1] 1.126598

norm(c – psd,type=”F”)
[1] 8.827865

So the nearPSD method is much much faster. This is due to the fact that the nearPD method is iterative.  However the “difference” between the output matrix is larger. Qualitatively, the nearPD drastically increases the correlation between 1 and 2 — from 0 to .79. The nearPSD keeps that correlation lower, .285, at the expense of the correlation between 1 and 2 and the rest of the values.