clear

n=25;
T=200;
r00=3;
m1=min(n,T);
R=1000;
for ii=1:R
    
    r_est=zeros(1,R);
    c=log(m1)/(10*m1);
    
    
    
    %%%one strong factor
%    F1=normrnd(0,2,T,1);
%    F2=normrnd(0,1,T,r00-1);
%    F0=[F1,F2];
    
    %%%normal
    F0=normrnd(0,1,T,r00);
    Lam0=randn(n,r00);
    % %%%%%%%% i.i.d. case1
    %u=normrnd(0,1,n,T);
    %beta=0;
    %J=0;
    %rho=0;
    
     %%%% serial correlation case2
     %beta=0;
     %  J=0;
     % rho=0.7;
     %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);
    
    
     %%% cross-sectional case3
         rho=0;
          beta=0.2;
         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
    
    %%%% cross-sectional case4
    %rho=0.7;
    %beta=0.2;
    %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
    
    e=sqrt((1-(rho)^2)/(1+2*J*(beta)^2))*u;
    
    y=Lam0*F0'+e;
    %y=y1-ones(n,1)*mean(y1,1)-mean(y1,2)*ones(T,1)'+mean(mean(y1))*ones(n,T) ;
    y = y- kron(mean(y,2), ones(T,1)');
    %y = y ./ kron(std(y')', ones(T,1)');%产生随机数据，这个过程没有问题
    
    
    %计算omega和dft
    %对时间T做dft
    approxt=(1:min(floor(3.3*sqrt(T))+3,floor(T/2-0.1)))';
    mt=length(approxt);
    Xft=(1/sqrt(T))*y*exp(-sqrt(-1)*2*pi*(1:T)'*approxt'/T);
    %对个体n做dft
    approx=(1:min(floor(3.3*sqrt(n))+3,floor(n/2-0.1)))';
    m=length(approx);
    Xf=(1/sqrt(n))*y'*exp(-sqrt(-1)*2*pi*(1:n)'*approx'/n);
    %dft之前的特征值
    covy=(y*y')/(n*T);
    [F,D]=eig(covy);
    eigv1=diag(D);
    eigv(:,ii)=(fliplr(eigv1'))';
    
    rmax=2*r00;
    
    %dft之后的特征值
    covyt=Xft*Xft'/mt/T;
    gammat1=eig(covyt);
    gammat(:,ii)=fliplr((gammat1'))';
    %fin_gammat=mean(gammat,2);
    
    %covy1=Xf*Xf'/m/T;
    covy1=Xf*Xf'/m/T;
    gamma2=eig(covy1);
    gamma(:,ii)=fliplr((gamma2'))';
    %fin_gamma=mean(gamma,2);
    
    
    [ft,~]= eigs(covyt, m1);
    for k=1:m1
        v=ft(:,k);
        delta_t(k)=v'*y*(y')*v/n/T;
    end
    %deltaT(:,ii)=flipud(real(delta_t'));%%%这是我转置后的delta
    deltaT(:,ii)=real(delta_t');
    %fin_deltaT=mean(deltaT,2);
    
    
    [f1,~]= eigs(covy1, m1);
    for k=1:m1
        v=f1(:,k);
        delta1(k)=v'*(y')*y*v/n/T;
    end
    %delta(:,ii)=flipud(real(delta1'));%%%这是我转置后的delta
    delta(:,ii)=real(delta1');
    %fin_delta=mean(delta,2);
end
fin_eigv=mean(eigv,2);
fin_gamma=mean(gamma,2);
fin_gammat=mean(gammat,2);
fin_delta=mean(delta,2);
fin_deltaT=mean(deltaT,2);

fin_eigv(1:10)'
fin_gamma(1:10)'
%fin_deltaT(1:10)'
%figure(1)

X1=(1:10);
Y0=fin_eigv(1:10);

%P1=fin_delta(1:10);
%P2=fin_deltaT(1:10);
%Y3=gamma(1:10)';
%plot(X1,Y0,'r-o',X1,P1,'b-*')
%title('Eigenvalues vs. pseudo-eigenvalues')
%xlabel('k')
%ylabel('values')
%legend('eigenvalues','pseudo-eigenvalues')

figure(2)

Y1=fin_gamma(1:10);
%Y2=fin_gammat(1:10);
%Y3=gamma(1:10)';
%plot(X1,Y0,'r-o',X1,Y1,'b-*',X1,Y2,'k-x')
plot(X1,Y0,'r-o',X1,Y1,'b-*')
title('Eigenvalues vs. F-eigenvalues')
xlabel('k')
ylabel('values')
legend('eigenvalues','F-eigenvalues')

fin_eigv(r00)/fin_eigv(r00+1)
fin_gamma(r00)/fin_gamma(r00+1)
