% This function performs Principal Components Analysis
% Input List: This function takes in 1 input: A, the matrix of observations
% Output List: This function has three outputs: the loadings matrix, the percentage of
% variance explained by each factor, and the cumulative variance explained
% by the number of factors until that point

function CumulativePercVarExplainedVarCov = PCAfil(A)

[nobs,n] = size(A);

% Compute Variance/Covariance and Correlation Matrices 
[VarCovMat,CorMat] = VarianceCovarianceCorrelationMatrix(A);

% Compute EigenValues and EigenVectors for both VarianceCovariance and Correlation Matrices
[EigenVectorsVarCovMat,EigenValuesVarCovMat] = eig(VarCovMat);
[EigenVectorsCorMat,EigenValuesCorMat] = eig(CorMat);

% Initialize Explained Variance Vectors
PercVarExplainedVarCov = zeros(n,1);
PercVarExplainedCor = zeros(n,1);
CumulativePercVarExplainedVarCov = zeros(n,1);
CumulativePercVarExplainedCor = zeros(n,1);


% Calculating Explained Variance
for i=n:-1:1
   PercVarExplainedVarCov(i) = EigenValuesVarCovMat(i,i)/trace(EigenValuesVarCovMat); 
   PercVarExplainedCor(i) = EigenValuesCorMat(i,i)/trace(EigenValuesCorMat);
   % Calculating Cumulative Explained Variance
   if(i==n)
        CumulativePercVarExplainedVarCov(i) = PercVarExplainedVarCov(i);
        CumulativePercVarExplainedCor(i) = PercVarExplainedCor(i);
   else
        CumulativePercVarExplainedVarCov(i) = CumulativePercVarExplainedVarCov(i+1) + PercVarExplainedVarCov(i);
        CumulativePercVarExplainedCor(i) = CumulativePercVarExplainedCor(i+1) + PercVarExplainedCor(i);
   end
end

% Formatting Output

% Graphs
FigHandle = figure('Position', [100, 100, 1049, 895]);

plotrows = ceil(n/2);
plotcols = 2;

for i=1:n
   componentindex = abs(i-(n+1));
   % Variance Covariance Plots 
   subplot(plotrows,plotcols,i);
   barVarCov = bar(EigenVectorsVarCovMat(:,componentindex));
   set(barVarCov,'BarWidth',0.6);
   hold on;
   barCor = bar(EigenVectorsCorMat(:,componentindex),'FaceColor', 'r', 'EdgeColor', 'r');
   set (barCor,'BarWidth',0.3);
   hold off;
   titlestring = {'Loadings component '};
   titlenum = num2str(i);
   titlecomplete = strcat(titlestring,titlenum);
   title(titlecomplete);
end

leg = legend ('VarCov','Corr');
rect = [0.49, 0.94, 0.05, 0.05];
set(leg,'Position',rect);

figure;
flipPercVarExplainedVarCov = flipud(PercVarExplainedVarCov);
line1 = plot(flipPercVarExplainedVarCov);
set(line1,'LineWidth',2)
hold on;
flipPercVarExplainedCor = flipud(PercVarExplainedCor);
line2 = plot(flipPercVarExplainedCor);
set(line2,'Color','red','LineWidth',2)
hold off;
set(gca,'XTick',0:1:(n+1));
axis([0 (n+1) 0 1]);
title('Portion of Variance explained by each component');
leg = legend ('VarCov','Corr');
rect = [0.70, 0.80, 0.05, 0.05];
set(leg,'Position',rect);

% OutPut Tables - Variance Covariance 

fprintf(1,'**Principal Components Analysis** - Variance Covariance Matrix\n');
fprintf('\n');
fprintf(1,'The number of variables is %d and the number of observations is %d\n',n, nobs);
fprintf('\n');
fprintf(1,'Importance of Components:');

for i=1:n
    fprintf('  Comp %d\t', i);
end

fprintf('\n');
fprintf(1,'Proportion of Variance\t');

for i=n:-1:1
   fprintf('%6.8f\t',PercVarExplainedVarCov(i));
end

fprintf('\n');
fprintf(1,'Cumulative Proportion\t');

for i=n:-1:1
    fprintf('%6.8f\t',CumulativePercVarExplainedVarCov(i));
end

fprintf('\n');
fprintf('\n');
fprintf(1,'Loadings:\t \t');

for i=1:n
    fprintf('  Comp %d\t', i);
end

fprintf('\n');

for i=1:n
    fprintf('Variable %d\t',i);
    fprintf('\t');
    for j=n:-1:1
      fprintf('%6.8f\t', EigenVectorsVarCovMat(i,j));
    end
    fprintf('\n');
end

% OutPut Tables - Correlation 

fprintf(1,'\n \n **Principal Components Analysis** - Correlation Matrix\n');
fprintf('\n');
fprintf(1,'The number of variables is %d and the number of observations is %d\n',n, nobs);
fprintf('\n');
fprintf(1,'Importance of Components:');

for i=1:n
    fprintf('  Comp %d\t', i);
end

fprintf('\n');
fprintf(1,'Proportion of Variance\t');

for i=n:-1:1
   fprintf('%6.8f\t',PercVarExplainedCor(i));
end

fprintf('\n');
fprintf(1,'Cumulative Proportion\t');

for i=n:-1:1
    fprintf('%6.8f\t',CumulativePercVarExplainedCor(i));
end

fprintf('\n');
fprintf('\n');
fprintf(1,'Loadings:\t \t');

for i=1:n
    fprintf('  Comp %d\t', i);
end

fprintf('\n');

for i=1:n
    fprintf('Variable %d\t',i);
    fprintf('\t');
    for j=n:-1:1
      fprintf('%6.8f\t', EigenVectorsCorMat(i,j));
    end
    fprintf('\n');
end

end


function [VarCovMat, CorMat] = VarianceCovarianceCorrelationMatrix (M)

[n,m] = size (M);

VarCovMat = zeros(m,m);
CorMat = zeros(m,m);
format long;

meanvar = mean(M);
stdvar = std(M);

%VarCovMAt cells loop
    for i=1:m
        for j=1:m
        %Variance/Covariance calculation
        varcovar = 0;    
            for k=1:n
            rowcalc = (M(k,i) - meanvar(i))*(M(k,j) - meanvar(j));
            varcovar = varcovar + rowcalc; 
            end
        varcovar = varcovar/(n-1);
        VarCovMat(i,j) = varcovar;
        CorMat (i,j) = VarCovMat(i,j)/(stdvar(i)*stdvar(j));
        end
    end
    
%Using Matlab functions
%VarCovMatML = cov(M);
%CorMatML = corrcov(VarCovMatML);

end

