
% X is observed
% r is the true number of true factors
% F is T by r matrix of true factors
% Lambda N by r is the true loading matrix
% C=F*Lambda' T by N is the true common component
% chat is the estimated common component

function [pc1,ic1]=mypcp(X,kmax);
T=size(X,1);
N=size(X,2);

m1=ones(T,1)*mean(X);
m2=ones(T,1)*std(X);
X=(X-m1)./m2;


NT=N*T;
NT1=N+T;
CT=zeros(3,kmax+1);
ii=[(1:1:kmax) 0];
CT(1,:)=log(NT/NT1)*ii*NT1/NT;
CT(2,:)=(NT1/NT)*log(min([N;T]))*ii;
GCT=min([N;T]);
CT(3,:)=ii*log(GCT)/GCT;
PC1=zeros(size(CT,1),kmax+1);
IC1=zeros(size(CT,1),kmax+1);
rr=zeros(1,kmax);
Sigma=zeros(1,kmax+1);
logSig=zeros(1,kmax+1);
chat=zeros(T,N);
XX=X*X';
[Fhat0,eigval,Fhat1]=svd(XX');
for i=kmax:-1:1;
Fhat=Fhat0(:,1:i);
%lambda=Fhat\X(:,j);
lambda=Fhat'*X;
chat=Fhat*lambda;
ehat=X-chat;
Sigma(i)=mean(sum(ehat.*ehat/T));
PC1(:,i)=Sigma(i)+CT(:,i)*Sigma(kmax);
IC1(:,i)=log(Sigma(i))+CT(:,i);
end;
Sigma(kmax+1)=mean(sum(X.*X/T));
PC1(:,kmax+1)=Sigma(kmax+1);
IC1(:,kmax+1)=log(Sigma(kmax+1));
pc1=minindc(PC1')';
pc1=pc1 .*(pc1 <= kmax);
ic1=minindc(IC1')';
ic1=ic1 .*(ic1 <= kmax);
