clear
tic;

R=1000;

N=100;
T=200;
r=5;
rho=(0:0.1:1);
Beta=(0:0.1:1);

for ii=1:length(rho)
    for kk=1:length(Beta)
        for jj=1:R
            
            m1=min(N,T);
            %%%one strong factor
            F1=normrnd(0,2,1,T);
            F2=normrnd(0,1,r-1,T);
            F0=[F1;F2];
            
            %%%normal
            %F0=normrnd(0,1,r,T);
            
            Lam0=normrnd(0,1,N,r);
            J=max(10,floor(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(ii)*u(i,t)+uu(i,t+1)+sum(Beta(kk)*uu(max(i-J,1):i-1,t))+sum(Beta(kk)*uu(i+1:min(i+J,N),t));
                end
            end        
            
            e=sqrt((1-(rho(ii))^2)/(1+2*J*(Beta(kk))^2))*u;
            y=Lam0*F0+sqrt(2)*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'))';
            
            rmax=min(mN,min(mT,2*r));
            
            [f1,~]= eigs(covyN, rmax);
            for k=1:rmax
                v=f1(:,k);
                delta(k)=real(v'*(y')*y*v)/N/T;
            end
            
            %%ER
            er=eigv(1:rmax-1)./eigv(2:rmax);
            rER1(jj)=find(er==max(er));
            
            %%FNER
            fner=gammaN(1:rmax-1)./gammaN(2:rmax);
            rFNER1(jj)=find(fner==max(fner));
            
            % %FTER
            fter=gammaT(1:rmax-1)./gammaT(2:rmax);
            rFTER1(jj)=find(fter==max(fter));
            
            
            % %FNTER
            fnter=gammaNT(1:rmax-1)./gammaNT(2:rmax);
            rFNTER1(jj)=find(fnter==max(fnter));
            
            %%PER
            per=delta(1:rmax-1)./delta(2:rmax);
            rPER1(jj)=find(per==max(per));
            
        end
        rER(ii,kk)=mean(rER1);
        rFNER(ii,kk)=mean(rFNER1);
        rFTER(ii,kk)=mean(rFTER1);
        rFNTER(ii,kk)=mean(rFNTER1);
        rPER(ii,kk)=mean(rPER1);
    end
end
figure(1)
r_ER=rER(:);
r_hat1=reshape(r_ER,length(rho),length(Beta));
peaks;colormap white; colorbar ;
mesh(rho,Beta,r_hat1);
surf(rho,Beta,r_hat1);
%plot(par,r_est)
xlabel('rho');
ylabel('Beta');
zlim([0 5]);
%zlabel('\hat{r}');
zlabel('r');
title('ER');

figure(2)
r_FNER=rFNER(:);
r_hat2=reshape(r_FNER,length(rho),length(Beta));
peaks;colormap white; colorbar ;
mesh(rho,Beta,r_hat2);
surf(rho,Beta,r_hat2);
%plot(par,r_est)
xlabel('rho');
ylabel('Beta');
zlim([0 5]);
%zlabel('\hat{r}');
zlabel('r');
title('FNER');

figure(3)
r_FTER=rFTER(:);
r_hat3=reshape(r_FTER,length(rho),length(Beta));
peaks;colormap white; colorbar ;
mesh(rho,Beta,r_hat3);
surf(rho,Beta,r_hat3);
%plot(par,r_est)
xlabel('rho');
ylabel('Beta');
zlim([0 5]);
%zlabel('\hat{r}');
zlabel('r');
title('FTER');

figure(4)
r_FNTER=rFNTER(:);
r_hat3=reshape(r_FNTER,length(rho),length(Beta));
peaks;colormap white; colorbar ;
mesh(rho,Beta,r_hat3);
surf(rho,Beta,r_hat3);
%plot(par,r_est)
xlabel('rho');
ylabel('Beta');
zlim([0 5]);
%zlabel('\hat{r}');
zlabel('r');
title('FNTER');

figure(5)
r_PER=rPER(:);
r_hat3=reshape(r_PER,length(rho),length(Beta));
peaks;colormap white; colorbar ;
mesh(rho,Beta,r_hat3);
surf(rho,Beta,r_hat3);
%plot(par,r_est)
xlabel('rho');
ylabel('Beta');
zlim([0 5]);
%zlabel('\hat{r}');
zlabel('r');
title('PER');
toc;