% This code provides information criteria selection based on Bai and Ng
% (2002). It is modified from the version of codes provided by Serena Ng on
% http://www.columbia.edu/~sn2294/research.html for Bai and Ng (2006).

% Inputs 
% X is observed data matrix
% kmax is the maximum number of factors considered
% jj specifies the information criteria that is imposed according to Bai
% and Ng (2002)
% DEMEAN=2 standardize data, DEMEAN=1, demean data, DEMEAN=0 raw data


% Outputs
% ic1 is the selected number of factors
% Fhat is T by r matrix of true factors
% chat is the estimated common component
% eigval is the eigenvalues associated with the principal component
% decomposition

function [ic1, chat,Fhat,eigval] = ICP(x,kmax,jj,DEMEAN)
    T=size(x,1);
    N=size(x,2);
    NT=N*T;
    NT1=N+T;
    CT=zeros(1,kmax);
    ii=1:1:kmax;
    if jj ==1;CT(1,:)=log(NT/NT1)*ii*NT1/NT;end;
    if jj==2; CT(1,:)=(NT1/NT)*log(min([N;T]))*ii;end;
    GCT=min([N;T]);
    if jj==3; CT(1,:)=ii*log(GCT)/GCT; end;
    if jj==4; CT(1,:)=2*ii/T; end;
    if jj==5; CT(1,:)=log(T)*ii/T;end;
    if jj==6; CT(1,:)=2*ii*NT1/NT; end;
    if jj==7; CT(1,:)=log(NT)*ii*NT1/NT;end;

    if DEMEAN == 0;
        X=x;
    elseif DEMEAN ==1;
        X=demean(x);
    elseif DEMEAN ==2;
        X=standard(x);
    end

    IC1=zeros(size(CT,1),kmax+1);
    Sigma=zeros(1,kmax+1);
    XX=X*X';
    [Fhat0,eigval,Fhat1]=svd(XX');
    sumeigval=cumsum(diag(eigval))/sum(diag(eigval));

    if jj < 8
        for i=kmax:-1:1
            Fhat=Fhat0(:,1:i);
            lambda=Fhat'*X;
            chat=Fhat*lambda;
            ehat=X-chat;
            Sigma(i)=mean(sum(ehat.*ehat/T));
            IC1(:,i)=log(Sigma(i))+CT(:,i);
        end
        
        Sigma(kmax+1)=mean(sum(X.*X/T));
        IC1(:,kmax+1)=log(Sigma(kmax+1));
        ic1=minindc(IC1')';
        ic1=ic1 .*(ic1 <= kmax);
    end
    
    if jj==8
      for j=1:size(sumeigval,1);
        if sumeigval(j) >= .6; ic1=j; break; end;% 0.6 modified by Tatevik
      end    
    end
    
    Fhat=[];
    Fhat=Fhat0(:,1:kmax);
    lambda=Fhat'*X;
    chat=Fhat*lambda;
    eigval=diag(eigval);
end

function x=seqa(a,b,c)
    x=linspace(a,(a+(c-1)*b),c)';
end

function pos=minindc(x)
    ncols=size(x,2);
    nrows=size(x,1);
    pos=zeros(ncols,1);
    seq=seqa(1,1,nrows);
    for i=1:ncols;
        dum=min(x(:,i));
        dum1= seq .* ( (x(:,i)-dum) ==0);
        pos(i)=sum(dum1);
    end
end