clear
tic;

R=1000;
theta=2;

A=zeros(30,8);
for Strongfactor=0:1
    for Case=1:4
        rER=zeros(1,R);
        rFNER=zeros(1,R);
        rFTER=zeros(1,R);
        rFNTER=zeros(1,R);
        rPER=zeros(1,R);
        rGR=zeros(1,R);
        rFNGR=zeros(1,R);
        rFTGR=zeros(1,R);
        rFNTGR=zeros(1,R);
        rPGR=zeros(1,R);
        rRRE=zeros(1,R);
        rFNRRE=zeros(1,R);
        rFTRRE=zeros(1,R);        
        rFNTRRE=zeros(1,R);
        rPRRE=zeros(1,R);
        rTCR=zeros(1,R);
        rFNTCR=zeros(1,R);
        rFTTCR=zeros(1,R);
        rFNTTCR=zeros(1,R);
        rPTCR=zeros(1,R);
        rTER=zeros(1,R);
        rNTER=zeros(1,R);        
        rTTER=zeros(1,R);        
        rNTTER=zeros(1,R);
        rPTER=zeros(1,R);
        rED=zeros(1,R);
        rFNED=zeros(1,R);
        rFTED=zeros(1,R);
        rFNTED=zeros(1,R);
        rPED=zeros(1,R);
        
        for ii=1:R
            N=randi([25,200],1);
            T=randi([25,200],1);
            r=randi([1,10],1);
            
            m1=min(N,T);
            M1=max(N,T);
            k1=m1;
                        
            if Strongfactor==1
                F1=normrnd(0,2,T,1);
                F2=normrnd(0,1,T,r-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);
                u(:,2:T)=rho*u(:,1:T-1)+uu(:,2:T);
            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)');
            
            
                        % PhiN
                        omegaN=(4:min(floor(3.3*sqrt(N))+3,floor(N/2-0.1)))';
                        omegaN=2*pi*omegaN/N;
                        mN=length(omegaN);
                        PhiN=zeros(mN,N);
                        for j=1:mN
                            for k=1:N
                                PhiN(j,k)=exp(-sqrt(-1)*omegaN(j)*k)/sqrt(N);
                            end
                        end
                        
                        %PhiT
                        omegaT=(4:min(floor(3.3*sqrt(T))+3,floor(T/2-0.1)))';
                        omegaT=2*pi*omegaT/T;
                        mT=length(omegaT);
                        PhiT=zeros(T,mT);
                        for j=1:mT
                            for k=1:T
                                PhiT(k,j)=exp(-sqrt(-1)*omegaT(j)*k)/sqrt(T);
                            end
                        end
                        
                        yN=PhiN*y;
                        yT=y*PhiT;
                        yNT=PhiN*y*PhiT;
                        
                        covy=(y'*y)/N/T;
                        eigv=eig(covy);
                        eigv=(fliplr(eigv'))';
                        
                        covyN=yN'*yN/mN/T;
                        gammaN=eig(covyN);
                        gammaN=fliplr((gammaN'))';
                        
                        covyT=yT'*yT/mT/N;
                        gammaT=eig(covyT);
                        gammaT=fliplr((gammaT'))';
                        
                        covyNT=yNT'*yNT/mT/mN;
                        gammaNT=eig(covyNT);
                        gammaNT=fliplr((gammaNT'))';
                        
            [f1,~]= eigs(covyN, m1);
            for k=1:m1
                v=f1(:,k);
                delta1(k)=v'*(y')*y*v/N/T;
            end
            delta=real(delta1');
            
            m1hat=min(mN,T);
            m1tilde=min(N,mT);
            m1bar=min(mN,mT);
            
            c=log(m1)/(10*m1); 
            chat=log(m1hat)/10/m1hat;
            ctilde=log(m1tilde)/10/m1tilde;
            cbar=log(m1bar)/10/m1bar;
            
            rmax=min(mT,min(mN,r*2))-1;
            
            % %ER
            er=eigv(1:rmax-1)./eigv(2:rmax);
            r_er=find(er==max(er));
            if r_er==r
                rER(ii)=1;
            end
            
            % %FNER
            fner=gammaN(1:rmax-1)./gammaN(2:rmax);
            r_fner=find(fner==max(fner));
            if r_fner==r
                rFNER(ii)=1;
            end
            
            % %FTER
            fter=gammaT(1:rmax-1)./gammaT(2:rmax);
            r_fter=find(fter==max(fter));
            if r_fter==r
                rFTER(ii)=1;
            end
            
            % %FNTER
            fnter=gammaNT(1:rmax-1)./gammaNT(2:rmax);
            r_fnter=find(fnter==max(fnter));
            if r_fnter==r
                rFNTER(ii)=1;
            end
            
            % %PER
            per=delta(1:rmax-1)./delta(2:rmax);
            r_per=find(per==max(per));
            if r_per==r
                rPER(ii)=1;
            end
                        
            % % 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));
            r_gr=find(grr0==max(grr0));
            if r_gr==r
                rGR(ii)=1;
            end
            
            % % FNGR
            VV=sum(gammaN(1:m1hat))*ones(m1hat,1)-cumsum(gammaN(1:m1hat));
            grr0=zeros(rmax,1);
            grr0(1)=log(sum(gammaN)/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));
            r_fngr=find(grr0==max(grr0));
            if r_fngr==r
                rFNGR(ii)=1;
            end
            
            % % FTGR
            VV=sum(gammaT(1:m1tilde))*ones(m1tilde,1)-cumsum(gammaT(1:m1tilde));
            grr0=zeros(rmax,1);
            grr0(1)=log(sum(gammaT)/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));
            r_ftgr=find(grr0==max(grr0));
            if r_ftgr==r
                rFTGR(ii)=1;
            end
            
            % % FNTGR
            VV=sum(gammaNT(1:m1bar))*ones(m1bar,1)-cumsum(gammaNT(1:m1bar));
            grr0=zeros(rmax,1);
            grr0(1)=log(sum(gammaNT)/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));
            r_fntgr=find(grr0==max(grr0));
            if r_fntgr==r
                rFNTGR(ii)=1;
            end
           
            % %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));
            r_pgr=find(grr0==max(grr0));
            if r_pgr==r
                rPGR(ii)=1;
            end
             
            %  RRE
            rre=(eigv(1:k1-1)+c)./(eigv(2:k1)+c);
            r_rre=find(rre==max(rre));
            if r_rre==r
                rRRE(ii)=1;
            end
            
            % FNRRE
            frr3=(gammaN(1:m1hat-1)+chat)./(gammaN(2:m1hat)+chat);
            r_fnrre=find(frr3==max(frr3));
            if r_fnrre==r
                rFNRRE(ii)=1;
            end
            
            % FTRRE
            frr3=(gammaT(1:m1tilde-1)+ctilde)./(gammaT(2:m1tilde)+ctilde);
            r_ftrre=find(frr3==max(frr3));
            if r_ftrre==r
                rFTRRE(ii)=1;
            end
            
            
            % FNTRRE
            frr3=(gammaNT(1:m1bar-1)+cbar)./(gammaNT(2:m1bar)+cbar);
            r_fntrre=find(frr3==max(frr3));
            if r_fntrre==r
                rFNTRRE(ii)=1;
            end
            
            %  PRRE
            prre=(delta(1:k1-1)+c)./(delta(2:k1)+c);
            r_prre=find(prre==max(prre));
            if r_prre==r
                rPRRE(ii)=1;
            end
            
            
            % % TCR
            VV=sum(eigv(1:m1))*ones(m1,1)-cumsum(eigv(1:m1));
            tcr0=zeros(rmax,1);
            tcr0(1)=log(1+eigv(1)/sum(eigv))/log(1+eigv(2)/VV(1));
            tcr0(2:rmax)=log(1+eigv(2:rmax)./VV(1:rmax-1))./log(1+eigv(3:rmax+1)./VV(2:rmax));
            r_tcr=find(tcr0==max(tcr0));
            if r_tcr==r
                rTCR(ii)=1;
            end
            
            % FNTCR
            VV=sum(gammaN(1:m1hat))*ones(m1hat,1)-cumsum(gammaN(1:m1hat));
            tcr0=zeros(rmax,1);
            tcr0(1)=log(1+gammaN(1)/sum(gammaN))/log(1+gammaN(2)/VV(1));
            tcr0(2:rmax)=log(1+gammaN(2:rmax)./VV(1:rmax-1))./log(1+gammaN(3:rmax+1)./VV(2:rmax));
            r_fntcr=find(tcr0==max(tcr0));
            if r_fntcr==r
                rFNTCR(ii)=1;
            end
            
            % FTTCR
            VV=sum(gammaT(1:m1tilde))*ones(m1tilde,1)-cumsum(gammaT(1:m1tilde));
            tcr0=zeros(rmax,1);
            tcr0(1)=log(1+gammaT(1)/sum(gammaT))/log(1+gammaT(2)/VV(1));
            tcr0(2:rmax)=log(1+gammaT(2:rmax)./VV(1:rmax-1))./log(1+gammaT(3:rmax+1)./VV(2:rmax));
            r_fttcr=find(tcr0==max(tcr0));
            if r_fttcr==r
                rFTTCR(ii)=1;
            end
            
            
            % FNTTCR
            VV=sum(gammaNT(1:m1bar))*ones(m1bar,1)-cumsum(gammaNT(1:m1bar));
            tcr0=zeros(rmax,1);
            tcr0(1)=log(1+gammaNT(1)/sum(gammaNT))/log(1+gammaNT(2)/VV(1));
            tcr0(2:rmax)=log(1+gammaNT(2:rmax)./VV(1:rmax-1))./log(1+gammaNT(3:rmax+1)./VV(2:rmax));
            r_fnttcr=find(tcr0==max(tcr0));
            if r_fnttcr==r
                rFNTTCR(ii)=1;
            end
            
            % %PTCR
            VV=sum(delta(1:m1))*ones(m1,1)-cumsum(delta(1:m1));
            tcr0=zeros(rmax,1);
            tcr0(1)=log(1+delta(1)/sum(delta))/log(1+delta(2)/VV(1));
            tcr0(2:rmax)=log(1+delta(2:rmax)./VV(1:rmax-1))./log(1+delta(3:rmax+1)./VV(2:rmax));
            r_ptcr=find(tcr0==max(tcr0));
            if r_ptcr==r
                rPTCR(ii)=1;
            end
            
            
            %  BIC
            logei=log(eigv(1:k1,1)+1)-eigv(1:k1,1);
            flogei=T/2*cumsum(logei)/sum(logei)-2*log(T)*cumsum(ones(k1,1))/k1;
            r_bic=find(flogei==max(flogei));
            if r_bic==r
                rBIC(ii)=1;
            end
            
            %PBIC
            logei=log(delta(1:k1,1)+1)-delta(1:k1,1);
            flogei=T/2*cumsum(logei)/sum(logei)-2*log(T)*cumsum(ones(k1,1))/k1;
            r_pbic=find(flogei==max(flogei));
            if r_pbic==r
                rPBIC(ii)=1;
            end
            
            
            % %TER
            eigv08=chi2cdf((eigv+c).^2,1);
            ed=eigv08(1:m1-1)./eigv08(2:m1);
            r_ter=find(ed==max(ed));
            if r_ter==r
                rTER(ii)=1;
            end
            
            
            % NTER
            feigv07=chi2cdf((gammaN+chat).^2,1);
            frr7=feigv07(1:m1hat-1)./feigv07(2:m1hat);
            r_nter=find(frr7==max(frr7));
            if r_nter==r
                rNTER(ii)=1;
            end
            
            
            % TTER
            feigv07=chi2cdf((gammaT+ctilde).^2,1);
            frr7=feigv07(1:m1tilde-1)./feigv07(2:m1tilde);
            r_tter=find(frr7==max(frr7));
            if r_tter==r
                rTTER(ii)=1;
            end
            
            % NTTER
            feigv07=chi2cdf((gammaNT+cbar).^2,1);
            frr7=feigv07(1:m1bar-1)./feigv07(2:m1bar);
            r_ntter=find(frr7==max(frr7));
            if r_ntter==r
                rNTTER(ii)=1;
            end
            
            %%PTER
            delta08=chi2cdf((delta+c).^2,1);
            ped=delta08(1:m1-1)./delta08(2:m1);
            r_pter=find(ped==max(ped));
            if r_pter==r
                rPTER(ii)=1;
            end
            
            
            % %ED
            eigv08=chi2cdf((0.43*log(m1)*eigv).^2,1);
            ed=eigv08(1:m1-1)-eigv08(2:m1);
            r_ed=find(ed==max(ed));
            if r_ed==r
                rED(ii)=1;
            end
            
            
            %%FNED
            gamma08=chi2cdf((0.43*log(m1hat)*gammaN).^2,1);
            fed=gamma08(1:m1hat-1)-gamma08(2:m1hat);
            r_fned=find(fed==max(fed));
            if r_fned==r
                rFNED(ii)=1;
            end
            
            
            %%FTED
            gamma08=chi2cdf((0.43*log(m1tilde)*gammaT).^2,1);
            fed=gamma08(1:m1tilde-1)-gamma08(2:m1tilde);
            r_fted=find(fed==max(fed));
            if r_fted==r
                rFTED(ii)=1;
            end
            
            
            %%FNTED
            gamma08=chi2cdf((0.43*log(m1bar)*gammaNT).^2,1);
            fed=gamma08(1:m1bar-1)-gamma08(2:m1bar);
            r_fnted=find(fed==max(fed));
            if r_fnted==r
                rFNTED(ii)=1;
            end
            
            %%PED
            delta08=chi2cdf((0.43*log(m1)*delta).^2,1);
            ped=delta08(1:m1-1)-delta08(2:m1);
            r_ped=find(ped==max(ped));
            if r_ped==r
                rPED(ii)=1;
            end
            
            
        end
        A(1,4*Strongfactor+Case)=round(sum(rER)*100/R);
        A(2,4*Strongfactor+Case)=round(sum(rFNER)*100/R);
        A(3,4*Strongfactor+Case)=round(sum(rFTER)*100/R);
        A(4,4*Strongfactor+Case)=round(sum(rFNTER)*100/R);
        A(5,4*Strongfactor+Case)=round(sum(rPER)*100/R);
        A(6,4*Strongfactor+Case)=round(sum(rGR)*100/R);
        A(7,4*Strongfactor+Case)=round(sum(rFNGR)*100/R);
        A(8,4*Strongfactor+Case)=round(sum(rFTGR)*100/R);
        A(9,4*Strongfactor+Case)=round(sum(rFNTGR)*100/R);
        A(10,4*Strongfactor+Case)=round(sum(rPGR)*100/R);
        A(11,4*Strongfactor+Case)=round(sum(rRRE)*100/R);
        A(12,4*Strongfactor+Case)=round(sum(rFNRRE)*100/R);
        A(13,4*Strongfactor+Case)=round(sum(rFTRRE)*100/R);
        A(14,4*Strongfactor+Case)=round(sum(rFNTRRE)*100/R);
        A(15,4*Strongfactor+Case)=round(sum(rPRRE)*100/R);
        A(16,4*Strongfactor+Case)=round(sum(rTCR)*100/R);
        A(17,4*Strongfactor+Case)=round(sum(rFNTCR)*100/R);
        A(18,4*Strongfactor+Case)=round(sum(rFTTCR)*100/R);
        A(19,4*Strongfactor+Case)=round(sum(rFNTTCR)*100/R);
        A(20,4*Strongfactor+Case)=round(sum(rPTCR)*100/R);
        A(21,4*Strongfactor+Case)=round(sum(rTER)*100/R);
        A(22,4*Strongfactor+Case)=round(sum(rNTER)*100/R);
        A(23,4*Strongfactor+Case)=round(sum(rTTER)*100/R);
        A(24,4*Strongfactor+Case)=round(sum(rNTTER)*100/R);
        A(25,4*Strongfactor+Case)=round(sum(rPTER)*100/R);
        A(26,4*Strongfactor+Case)=round(sum(rED)*100/R);
        A(27,4*Strongfactor+Case)=round(sum(rFNED)*100/R);
        A(28,4*Strongfactor+Case)=round(sum(rFTED)*100/R);
        A(29,4*Strongfactor+Case)=round(sum(rFNTED)*100/R);
        A(30,4*Strongfactor+Case)=round(sum(rPED)*100/R);
    end
end
xlswrite('accucomp.xlsx',A,theta,'B2:I31');

toc;
