function y = Pca(usf, ukf, gmf)
% Principal component analysis of factors
[M, N] = size(usf);
myscore = nan(M,N);
mypca = nan(M,3);
myr2 = nan(M,3);
for i = 1:M
%     if mod(i,2000)==0
%         disp(i);
%     end
    myy = [usf(i,:)' ukf(i,:)' gmf(i,:)'];
    [~,score,~,~,expl] = pca(myy);
    mypca(i,:) = expl;
    myscore(i,:) = score(:,1);
    ybar = mean(myy,1);  % 1x3
    myx = [ones(N,1) score(:,1)];
    sst = sum((myy - ybar).^2, 1); % 1x3
    e1 = myy(:,1) - myx*((myx'*myx)\myx'*myy(:,1));  % Nx1
    e2 = myy(:,2) - myx*((myx'*myx)\myx'*myy(:,2));
    e3 = myy(:,3) - myx*((myx'*myx)\myx'*myy(:,3));
    ssr1 = e1'*e1;
    ssr2 = e2'*e2;
    ssr3 = e3'*e3;
    rsq1 = 1 - ssr1/sst(1) ;
    rsq2 = 1 - ssr2/sst(2) ;
    rsq3 = 1 - ssr3/sst(3) ;
    myr2(i,:) = [rsq1 rsq2 rsq3];
end

y.pca = mean(mypca,1);
y.r2 = mean(myr2,1);

