clear
tic;

N=50;
T=50;
r=5;
theta=2;

Strongfactor=0; %Strongfactor=0 if there is no dominant factor; Strongfactor=1 if there is one dominant factor; Strongfactor=2 if there is one weak factor.
Case=4; %Case=1 for not correlated errors (rho=beta=J=0); Case=2 for serially correlated errors (rho=0.5, beta=J=0);
%Case=3 for cross-sectionally correlated errors (rho=0, beta=0.5,
%J=max(10,N/20); Case 4 for both serially and cross-sectionally
%corrected errors (rho=0.5, beta=0.5, J=max(10,N/20).

R=10; %R is the number of replications

m1=min(N,T);
c=log(m1)/10/m1;

r_est=zeros(1,R);
for ii=1:R
    if Strongfactor==1
        F1=normrnd(0,2,T,1);
        F2=normrnd(0,1,T,r-1);
        F0=[F1,F2];
    elseif Strongfactor==2
        F1=normrnd(0,1,T,r-1);
        F2=normrnd(0,0.5,T,1);
        F0=[F1,F2];
    else
        F0=normrnd(0,1,T,r);
    end
    if Case==1
        rho=0;
        beta=0;
        J=0;
        u=normrnd(0,1,N,T);
    elseif Case==2
        rho=0.5;
        beta=0;
        J=0;
        uu=normrnd(0,1,N,T);
        u=zeros(N,T);
        u(:,1)=uu(:,1);
        for t=1:T-1
            u(:,t+1)=rho*u(:,t)+uu(:,t+1);
        end
    elseif Case==3
        rho=0;
        beta=0.5;
        J=max(10,N/20);
        uu=normrnd(0,1,N,T);
        u=zeros(N,T);
        for i=1:N
            u(i,:)=uu(i,:)+sum(beta*uu(max(i-J,1):i-1,:))+sum(beta*uu(i+1:min(i+J,N),:));
        end
    else
        rho=0.5;
        beta=0.5;
        J=max(10,N/20);
        uu=normrnd(0,1,N,T);
        u=zeros(N,T);
        u(:,1)=uu(:,1);
        for i=1:N
            for t=1:T-1
                u(i,t+1)=rho*u(i,t)+uu(i,t+1)+sum(beta*uu(max(i-J,1):i-1,t+1))+sum(beta*uu(i+1:min(i+J,N),t+1));
            end
        end
    end
    
    Lam0=normrnd(0,1,N,r);
    e=sqrt((1-(rho)^2)/(1+2*J*(beta)^2))*u;
    y=Lam0*F0'+sqrt(theta)*e;
    y = y - kron(mean(y,2), ones(T,1)');
    
    
    approx=(4:min(floor(3.3*sqrt(N))+3,floor(N/2-0.1)))';
    Xf=(1/sqrt(N))*y'*exp(-sqrt(-1)*2*pi*(1:N)'*approx'/N);
    m=length(approx);
    rmax=8;
    kmax=rmax;
    
    
    %eigv1是原始特特征值，eigv0是r=0时候的第0个特征值，eigv是两个特征值合并，此时它的维数是(N+1)*1,gamma与delta也是一样的
    covy=(y'*y)/(N*T);
    [V,D]=eig(covy);
    eigv1=diag(D);
    eigv1=(fliplr(eigv1'))';
    eigv0=min(1/sqrt(log(m1)),m1/sqrt(log(m1))*eigv1(1));%mock eigv
    eigv=[eigv0;eigv1];
    
    
    delta=zeros(m1,1);
    covy1=Xf*Xf'/m/T;
    gamma=eig(covy1);
    gamma1=fliplr((gamma'))';
    gamma0=min(1/sqrt(log(m)),m1/sqrt(log(m))*gamma1(1));%mock gamma
    gamma=[gamma0;gamma1];
    
    [f1,~]= eigs(covy1, m1);
    for k=1:m1
        v=f1(:,k);
        delta1(k,1)=real(v'*(y')*y*v)/N/T;
    end
    delta0=min(1/sqrt(log(m)),m1/sqrt(log(m))*delta1(1));%mock delta
    delta=[delta0;delta1];
    
     % ER
      rr1=(eigv(1:rmax-1))./(eigv(2:rmax));    
      rER(ii)=find(rr1==max(rr1))-1;
    %FNER
    per=gamma(1:rmax-1)./gamma(2:rmax);
    rFER(ii)=find(per==max(per))-1;
    
    %PER
    per=delta(1:rmax-1)./delta(2:rmax);
    rPER(ii)=find(per==max(per))-1;
    
    %GR
     VV=sum(eigv(1:m1))*ones(m1,1)-cumsum(eigv(1:m1));
    grr0=zeros(rmax,1);
    grr0(1)=log(sum(eigv)/VV(1))/log(VV(1)/VV(2));
    grr0(2:rmax)=log(VV(1:rmax-1)./VV(2:rmax))./log(VV(2:rmax)./VV(3:rmax+1));
    rGR(ii)=find(grr0==max(grr0))-1;
    
    %FNGR
    VV=sum(gamma(1:m1))*ones(m1,1)-cumsum(gamma(1:m1));
    grr0=zeros(rmax,1);
    grr0(1)=log(sum(gamma)/VV(1))/log(VV(1)/VV(2));
    grr0(2:rmax)=log(VV(1:rmax-1)./VV(2:rmax))./log(VV(2:rmax)./VV(3:rmax+1));
    rFGR(ii)=find(grr0==max(grr0))-1;
 
   
    % PGR
    VV=sum(delta(1:m1))*ones(m1,1)-cumsum(delta(1:m1));
    grr0=zeros(rmax,1);
    grr0(1)=log(sum(delta)/VV(1))/log(VV(1)/VV(2));
    grr0(2:rmax)=log(VV(1:rmax-1)./VV(2:rmax))./log(VV(2:rmax)./VV(3:rmax+1));
    rPGR(ii)=find(grr0==max(grr0))-1;
    
    %ON
    threshold=(1+N^(-1/3))*(2.7*eigv1(kmax+1)-1.7*eigv1(2*kmax+1));
    rr7=max(find(eigv1>threshold,1,'last'));
    if length(rr7)==0;
        rON(ii)=0;
    else
         rON(ii)= rr7-1;
    end

    %ABC
    [rhat1, rhat2]=ABC_crit(y,kmax,floor(N/10),3, 0);
    rr10=rhat1;
    if length(rr10)==0;
        rABC(ii)=0;
    else
         rABC(ii)= rr10;
    end
end
%parameters

ABC=[mean(rABC),round(sum(rABC<r)*100/R),round(sum(rABC==r)*100/R),round(sum(rABC>r)*100/R)]

ER=[mean(rER),round(sum(rER<r)*100/R),round(sum(rER==r)*100/R),round(sum(rER>r)*100/R)]
FNER=[mean(rFER),round(sum(rFER<r)*100/R),round(sum(rFER==r)*100/R),round(sum(rFER>r)*100/R)]
PER=[mean(rPER),round(sum(rPER<r)*100/R),round(sum(rPER==r)*100/R),round(sum(rPER>r)*100/R)]

GR=[mean(rGR),round(sum(rGR<r)*100/R),round(sum(rGR==r)*100/R),round(sum(rGR>r)*100/R)]
FNGR=[mean(rFGR),round(sum(rFGR<r)*100/R),round(sum(rFGR==r)*100/R),round(sum(rFGR>r)*100/R)]
PGR=[mean(rPGR),round(sum(rPGR<r)*100/R),round(sum(rPGR==r)*100/R),round(sum(rPGR>r)*100/R)]

ON=[mean(rON),round(sum(rON<r)*100/R),round(sum(rON==r)*100/R),round(sum(rON>r)*100/R)]
[N, T, r]
toc;