%AnalyzeCorrs 
% This program goes through the steps of checking whether a 
% correlation matrix has impossible combinations of values
% due to missing data or small numbers of observations, 
% and then correcting these problems. 

clear all;
%*****************Simulation Parameters***********************
%Define the correlation matrix: you can use a 1/2 matrix as below
% (with zeros in the upper triangle) and fill in the other 1/2
% (see below). This is the matrix for the desert tortoise 
corrmx = [...
1	0	0	0	0	0	0	0	0	0	0;
0.469	1	0	0	0	0	0	0	0	0	0;
0.382	0.995	1	0	0	0	0	0	0	0	0;
-0.597	0.429	0.514	1	0	0	0	0	0	0	0;
-0.014	0.877	0.919	0.811	1	0	0	0	0	0	0;
-1	-0.496	-0.41	0.571	-0.017	1	0	0	0	0	0;
-0.704	-0.957	-0.925	-0.149	-0.7	0.726	1	0	0	0	0;
0.898	0.81	0.75	-0.182	0.428	-0.911	-0.945	1	0	0	0;
0.956	0.189	0.094	-0.806	-0.307	-0.946	-0.465	0.729	1	0	0;
-0.514	0.516	0.597	0.995	0.865	0.487	-0.247	-0.083	-0.743	1	0;
-0.994	-0.563	-0.481	0.505	-0.096	0.997	0.778	-0.941	-0.918	0.417	1];
%*************************************************************
np = max(size(corrmx)); % how many vital rates?
corrmx = corrmx + corrmx' - eye(np); % make a full, symmetrical matrix
[W,D] = eig(corrmx);  % find the eigenvalues and eigenvectors

disp('this is the eigenvalue matrix: if any elements are')
disp('negative, then some adjustment is needed')
disp(D)
disp('the corresponding eigenvectors are:')
disp(W)
disp('  ');
disp('the C12 matrix is:')  
C12 = W*(sqrt(abs(D)))*W'

disp('for this correlation matrix, based on only three') 
disp('time periods, all but the last 2 eigenvalues are')
disp('small, so we should redo the calculation of the') 
disp('correlation matrix, eliminating all the others:')
newD = D;
newD(:,1:(np-2)) = zeros(np,(np-2));
newcorrmx = W*newD*W';

disp('but this newcorrmx must now be corrected to be')
disp('a proper correlation matrix, not a covariance matrix')
for ii=1:np   
 for jj=1:np
  newcorrmx(ii,jj) = ...    
     newcorrmx(ii,jj)/((newcorrmx(ii,ii)*newcorrmx(jj,jj))^0.5);
 end;
end;

[W,D] = eig(newcorrmx); % corrected eigenvectors and values
C12 = W*(sqrt(abs(D)))*W';
disp('the new, corrected, correlation matrix is: ')
disp(newcorrmx)
disp('and the new C12 is: ')
disp(C12)
disp('finally, the proportional differences between the old')
disp('correlations and the new ones are:')
disp([(corrmx - newcorrmx)./corrmx])