%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure replication
% Purpose: we summarize codes generating the data for figures
% and the code for
% Figures are drawn based on the data generated from these codes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% choose figure ID
%
data_id=5;
% 1: power curves
% 2: size curves (for various N (across N))
% 3: size curves (for various lambda_f (across lambda_f))
% 4: p-value (subset FAR) curves (LL data with three factors)
% 5: p_value curves (kroenker data (consumption growth factor))
%    single weak factor
% 6: sFAR distribution comparison
%
%
fig_id=5;
% 1: power curves
% 2: size curves (for various N (across N))
% 3: size curves (for various lambda_f (across lambda_f))
% 4: p-value (subset FAR) curves (LL data with three factors)
% 5: p_value curves (kroenker data (consumption growth factor))
%    single weak factor
% 6: sFAR distribution comparison
% 7: qq-plot for returns
%%


%% Data Generation %%%%%%%%%%%%%%%
%% Power curve data:
% generate the data for power curves of FAR, sCLR, JGLS and GLSLM for a
% single factor model
% and we compare the finite sample and asymptotic distributions
if data_id==1
    clear;
    rng(1);
    data=xlsread('data_KZ2017.xls');    % load data
    %     load GLSLM_qunatile_T100N52.mat  % pre-stored 95% critical values
    %%%%%%%%  Tuning parameters
    lambda_f0 =0;
    T=100;
    T1=50;
    rdd_num=50000;
    power_disturb=(-0.10:0.01:0.10);
    asset_num_mius2 =43;
    % if r0=0, then we donot fix r (the identification strength)
    r0=15;
    %%%%%%%%
    power_length=length(power_disturb);
    reject_frequency_GLSLM=NaN(asset_num_mius2,2,power_length);
    reject_frequency_FAR=NaN(asset_num_mius2,2,power_length);
    reject_frequency_JGLS=NaN(asset_num_mius2,2,power_length);
    reject_frequency_CLR=NaN(asset_num_mius2,2,power_length);
    G=data(:,4);                        % consumption growth
    K=size(G,2);
    Vff=cov(G);
    F0=mvnrnd(zeros(1,K),Vff, T);
    for n=3:20:asset_num_mius2
        fprintf('%2.0f out of %2.0f\n', n, asset_num_mius2);
        tic;
        % step one: caliberating from the empirical data
        R=[data(:,9:9+n) data(:,7)];         % test assets
        N=size(R,2);
        T0=size(R,1);
        fprintf('N=%i', N);
        K=size(G,2);
        V=cov([G R]);
        Vc=V(1:1,1:1);                      % V_f
        Vr=V(2:end,2:end);                  % V_r
        Vcr=V(1:1,2:end);                   % V_fr
        Vrc=Vcr';                           % V_rf
        beta=Vrc/Vc;
        Omega=Vr-beta*Vc*beta';
        C_FAR=(T-K-N+1)/( (N-1))/(T-K-1);
        C_JGLS=(T-N+1)/(N-K-1)/(T-K-1);
        %%% renormalize the factor loading to fix r
        if r0>0
            Jn=[eye(N-1), -ones(N-1,1)];
            R=R*Jn';
            R_demean=(R-ones(T0,1)*mean(R));
            F_demean=(G-ones(T0,1)*mean(G));
            F_pluslambda= F_demean+ ones(T0,1)*lambda_f0;
            F_pluslambda_2=(F_pluslambda'* F_pluslambda);
            B_tilde=R'*F_pluslambda/(F_pluslambda'* F_pluslambda);
            B_hat=  R_demean'*G/(G'* G);
            Sigma_hat=(R_demean-F_demean*B_hat')'*(R_demean-F_demean*B_hat')/(T0-K);
            r=F_pluslambda_2/T0 *B_tilde'/Sigma_hat*B_tilde;
            beta= beta*sqrt(r0/r);
        end
        for j=1:power_length
            reject=0;
            rejectF=0;
            rejectJ=0;
            rejectC=0;
            reject2=0;
            rejectF2=0;
            rejectJ2=0;
            rejectC2=0;
            lambda_f= lambda_f0+power_disturb(j);
            parfor i=1:rdd_num
                F=F0;
                v=mvnrnd(zeros(1,N),Omega,T);
                R=(F+kron(ones(T,1),lambda_f'))*beta'+v;
                [~, p_val_F, p_val_chi]=FAR_CLR_JGLS_GLSLM_test_joint_rp(R,F,lambda_f,T1,2000, 1);
                reject=reject + (p_val_chi.GLSLM < 0.05);
                rejectF=rejectF + (p_val_chi.FAR < 0.05);
                rejectJ=rejectJ + (p_val_chi.JGLS < 0.05);
                rejectC=rejectC + (p_val_chi.sCLR < 0.05);
                reject2=reject2 + (p_val_F.GLSLM < 0.05);
                rejectF2=rejectF2 + (p_val_F.FAR < 0.05);
                rejectJ2=rejectJ2 + (p_val_F.JGLS < 0.05);
                rejectC2=rejectC2 + (p_val_F.sCLR < 0.05);
            end
            %           asymptotic
            reject_frequency_GLSLM(n,1,j)= (reject)/rdd_num;
            reject_frequency_FAR(n,1,j)=(rejectF)/rdd_num;
            reject_frequency_JGLS(n,1,j)=   (rejectJ)/rdd_num;
            reject_frequency_CLR(n,1,j)=   (rejectC)/rdd_num;
            %           finite sample
            reject_frequency_GLSLM(n,2,j)= (reject2)/rdd_num;
            reject_frequency_FAR(n,2,j)=(rejectF2)/rdd_num;
            reject_frequency_JGLS(n,2,j)=   (rejectJ2)/rdd_num;
            reject_frequency_CLR(n,2,j)=   (rejectC2)/rdd_num;
        end
        toc;
    end
    %     % figure for checking
    %     figure()
    %     hold on
    %     plot(reshape(reject_frequency_CLR(3,1,:), 1, []))
    %     plot(reshape(reject_frequency_CLR(3,2,:), 1, []))
    %     hold off
    % %   we save the data for drawing the plot
    % save Power_acrossNT100_N_3_step20_43_fix_r_15.mat reject_frequency_FAR  reject_frequency_GLSLM  reject_frequency_JGLS reject_frequency_CLR  power_disturb lambda_f0
end

%% Size curve data (for various N):
if data_id==2
    clear;
    data=xlsread('data_KZ2017.xls');    % load data
    load GLSLM_qunatile_T50.mat         % pre-stored
    G=data(:,4);                        % consumption growth factor
    K=size(G,2);
    Vff=cov(G);
    %%%%%%%%%%%   Tuning paramsters
    lambda_f =0;
    rdd_num=10000;
    T=100;
    T1=15;
    asset_num_mimus2 =43;
    % if r0=0, then we donot fix r (the identification strength)
    %     r0=15;
    r0=0;
    %%%%%%%%%%%
    reject_frequency_GLSLM=NaN(asset_num_mimus2,2);
    reject_frequency_FAR=NaN(asset_num_mimus2,2);
    reject_frequency_JGLS=NaN(asset_num_mimus2,2);
    reject_frequency_CLR=NaN(asset_num_mimus2,4);
    F0=mvnrnd(zeros(1,K),Vff,T);
    F0=F0-ones(T,1)*mean(F0);
    tic;
    for n=3:asset_num_mimus2
        fprintf('%2.0f out of %2.0f\n', n, asset_num_mimus2);
        % step one: caliberating from the empirical data
        R=[data(:,9:9+n) data(:,7)];         % test assets
        N=size(R,2);
        fprintf('N=%i', N);
        K=size(G,2);
        V=cov([G R]);
        Vc=V(1:1,1:1);                      % V_f
        Vr=V(2:end,2:end);                  % V_r
        Vcr=V(1:1,2:end);                   % V_fr
        Vrc=Vcr';                           % V_rf
        Sigma=V(2:end,2:end);                  % V_r
        beta=Vrc/Vc;
        Omega=Vr-beta*Vc*beta';
        % stage three: random dra from the distribution
        rejectG=0;
        rejectF=0;
        rejectJ=0;
        rejectC=0;
        rejectG2=0;
        rejectF2=0;
        rejectJ2=0;
        rejectC2=0;
        C_FAR=(T-K-N+1)/( (N-1))/(T-K-1);
        C_JGLS=(T-N+1)/(N-K-1)/(T-K-1);
        %%% renormalize the factor loading to fix r
        if r0>0
            Jn=[eye(N-1), -ones(N-1,1)];
            R=R*Jn';
            R_demean=(R-ones(T0,1)*mean(R));
            F_demean=(G-ones(T0,1)*mean(G));
            F_pluslambda= F_demean+ ones(T0,1)*lambda_f0;
            F_pluslambda_2=(F_pluslambda'* F_pluslambda);
            B_tilde=R'*F_pluslambda/(F_pluslambda'* F_pluslambda);
            B_hat=  R_demean'*G/(G'* G);
            Sigma_hat=(R_demean-F_demean*B_hat')'*(R_demean-F_demean*B_hat')/(T0-K);
            r=F_pluslambda_2/T0 *B_tilde'/Sigma_hat*B_tilde;
            beta= beta*sqrt(r0/r);
        end
        parfor i=1:rdd_num
            Fi=F0;
            v=mvnrnd(zeros(1,N),Omega,T);
            R=(Fi+kron(ones(T,1),lambda_f'))*beta'+v;
            [~, p_val_F, p_val_chi]=FAR_CLR_JGLS_GLSLM_test_joint_rp(R,Fi,lambda_f,T1,2000, 1);
            reject=reject + (p_val_chi.GLSLM < 0.05);
            rejectF=rejectF + (p_val_chi.FAR < 0.05);
            rejectJ=rejectJ + (p_val_chi.JGLS < 0.05);
            rejectC=rejectC + (p_val_chi.sCLR < 0.05);
            reject2=reject2 + (p_val_F.GLSLM < 0.05);
            rejectF2=rejectF2 + (p_val_F.FAR < 0.05);
            rejectJ2=rejectJ2 + (p_val_F.JGLS < 0.05);
            rejectC2=rejectC2 + (p_val_F.sCLR < 0.05);
        end
        reject_frequency_GLSLM(n,1)= (rejectG)/rdd_num;
        reject_frequency_FAR(n,1)=(rejectF)/rdd_num;
        reject_frequency_JGLS(n,1)=   (rejectJ)/rdd_num;
        reject_frequency_CLR(n,1)=   (rejectC)/rdd_num;
        reject_frequency_GLSLM(n,2)= (rejectG2)/rdd_num;
        reject_frequency_FAR(n,2)=(rejectF2)/rdd_num;
        reject_frequency_JGLS(n,2)=   (rejectJ2)/rdd_num;
        reject_frequency_CLR(n,2)=   (rejectC2)/rdd_num;
        %         fprintf('\n reject_frequency_CLR_1=%f \n reject_frequency_CLR_2=%f \n reject_frequency_CLR_2=%f \n',reject_frequency_CLR(n,:));
        toc;
        %% Recording Memo
    end
    %     save size_acrossNT100_new.mat reject_frequency_FAR  reject_frequency_GLSLM  reject_frequency_JGLS reject_frequency_CLR T asset_num_mimus2
    %     save size_acrossNT100_fizxedr_new2 reject_frequency_FAR  reject_frequency_GLSLM  reject_frequency_JGLS reject_frequency_CLR T asset_num_mimus2
    
end

%% Size curve data (for various lambda_f)
if data_id==3
    clear;
    rng(1);
    data=xlsread('data_KZ2017.xls');
    load GLSLM_qunatile_T100N52.mat  % per-stored GLSLM critical values
    G=data(:,4);                        % consumption growth
    K=size(G,2);
    Vff=cov(G);
    %%%%%   Tuning parameters
    power_disturb=(-0.15:0.01:0.15);
    asset_num_mius2 =43;
    lambda_f0 =0;
    T=100;
    T1=50;
    rdd_num=50;
    F0=mvnrnd(zeros(1,K),Vff, T);
    F0=F0-ones(T,1)*mean(F0);
    r0=15;
    %%%%%
    power_length=length(power_disturb);
    reject_frequency_GLSLM=NaN(asset_num_mius2,2,power_length);
    reject_frequency_FAR=NaN(asset_num_mius2,2,power_length);
    reject_frequency_JGLS=NaN(asset_num_mius2,2,power_length);
    reject_frequency_CLR=NaN(asset_num_mius2,3,power_length);
    for n=3:20:asset_num_mius2
        tic;
        fprintf('%2.0f out of %2.0f\n', n, asset_num_mius2);
        % step one: caliberating from the empirical data
        R=[data(:,9:9+n) data(:,7)];         % test assets
        N=size(R,2);
        T0=size(R,1);
        fprintf('N=%i', N);
        K=size(G,2);
        V=cov([G R]);
        Vc=V(1:1,1:1);                      % V_f
        Vr=V(2:end,2:end);                  % V_r
        Vcr=V(1:1,2:end);                   % V_fr
        Vrc=Vcr';                           % V_rf
        beta=Vrc/Vc;
        Omega=Vr-beta*Vc*beta';
        C_FAR=(T-K-N+1)/( (N-1))/(T-K-1);
        C_JGLS=(T-N+1)/(N-K-1)/(T-K-1);
        %%% renormalize the factor loading to fix r
        if r0>0
            Jn=[eye(N-1), -ones(N-1,1)];
            R=R*Jn';
            R_demean=(R-ones(T0,1)*mean(R));
            F_demean=(G-ones(T0,1)*mean(G));
            F_pluslambda= F_demean+ ones(T0,1)*lambda_f0;
            F_pluslambda_2=(F_pluslambda'* F_pluslambda);
            B_tilde=R'*F_pluslambda/(F_pluslambda'* F_pluslambda);
            B_hat=  R_demean'*G/(G'* G);
            Sigma_hat=(R_demean-F_demean*B_hat')'*(R_demean-F_demean*B_hat')/(T0-K);
            r=F_pluslambda_2/T0 *B_tilde'/Sigma_hat*B_tilde;
            beta= beta*sqrt(r0/r);
        end
        for j=1:power_length
            reject=0;
            rejectF=0;
            rejectJ=0;
            rejectC=0;
            reject2=0;
            rejectF2=0;
            rejectJ2=0;
            rejectC2=0;
            lambda_f= lambda_f0+power_disturb(j);
            parfor i=1:rdd_num
                Fi=F0;
                v=mvnrnd(zeros(1,N),Omega,T);
                R=(Fi+kron(ones(T,1),lambda_f'))*beta'+v;
                [~, p_val_F, p_val_chi]=FAR_CLR_JGLS_GLSLM_test_joint_rp(R,Fi,lambda_f,T1,2000, 1);
                reject=reject + (p_val_chi.GLSLM < 0.05);
                rejectF=rejectF + (p_val_chi.FAR < 0.05);
                rejectJ=rejectJ + (p_val_chi.JGLS < 0.05);
                rejectC=rejectC + (p_val_chi.sCLR < 0.05);
                reject2=reject2 + (p_val_F.GLSLM < 0.05);
                rejectF2=rejectF2 + (p_val_F.FAR < 0.05);
                rejectJ2=rejectJ2 + (p_val_F.JGLS < 0.05);
                rejectC2=rejectC2 + (p_val_F.sCLR < 0.05);
            end
            reject_frequency_GLSLM(n,1,j)= (reject)/rdd_num;
            reject_frequency_FAR(n,1,j)=(rejectF)/rdd_num;
            reject_frequency_JGLS(n,1,j)=   (rejectJ)/rdd_num;
            reject_frequency_CLR(n,1,j)=   (rejectC)/rdd_num;
            reject_frequency_GLSLM(n,2,j)= (reject2)/rdd_num;
            reject_frequency_FAR(n,2,j)=(rejectF2)/rdd_num;
            reject_frequency_JGLS(n,2,j)=   (rejectJ2)/rdd_num;
            reject_frequency_CLR(n,2,j)=   (rejectC2)/rdd_num;
            toc;
        end
    end
    %     figure()
    %     hold on
    %     plot(reshape(reject_frequency_CLR(3,1,:), 1, []))
    %     plot(reshape(reject_frequency_CLR(3,2,:), 1, []))
    %     hold off
    %     save Size_acrossNT100_N_3_step20_43_fix_r_15.mat reject_frequency_FAR  reject_frequency_GLSLM  reject_frequency_JGLS reject_frequency_CLR  power_disturb lambda_f0
end

%% sFAR P-value curve data:
%  using Lettau Ludvigson (2001) data
if data_id==4
    clear;
    rng(1);
    load LL_data.mat
    
    
    xxxx
    %%%%%%%% Tuning parameters
    lambda_f=(-20:0.1:20);
    %%%%%%%%
    length_lambda=length(lambda_f);
    pvalue_sFAR=NaN(3,2,length_lambda);
    for n=1:3
        switch(n)
            case 1
                F =[cay, Delta_c, Delta_c_cay];
            case 2
                F =[Delta_c, cay, Delta_c_cay];
            case 3
                F =[Delta_c_cay, cay, Delta_c];
        end
        R=FF25;         % test assets
        for i=1:length_lambda
            [sFAR,  pvalue_sFAR(n,2,i), pvalue_sFAR(n,1,i)]=sFAR_test_single_rp(R, F, lambda_f(i));
        end
    end
    %     save P_value_LL_new.mat pvalue_sFAR lambda_f
    %     figure()
    %     hold on
    %     plot(reshape(pvalue_sFAR(1,1,:), 1, []))
    %     plot(reshape(pvalue_sFAR(1,2,:), 1, []))
    %     hold off
end

%% single weak factor P-value curve data
if data_id==5
    clear;
    tic;
    rng(1);
    data=xlsread('data_KZ2017.xls');
    % load GLSLM_qunatile_T142.mat
    R=data(:,9:28); % test assets
    F=data(:,4);
    
    lambda_f=(-10:0.1:10);
    T1=25;
    m=50000;
    GLSLM_id=1;
    %%%%%%%%%%%%%
    
    
    
    length_lambda=length(lambda_f);
    pvalue_FAR=NaN(1,length_lambda);
    pvalue_sCLR=NaN(1,length_lambda);
    pvalue_JGLS=NaN(1,length_lambda);
    pvalue_GLSLM=NaN(1,length_lambda);
    pvalueF_FAR=NaN(1,length_lambda);
    pvalueF_sCLR=NaN(1,length_lambda);
    pvalueF_JGLS=NaN(1,length_lambda);
    pvalueF_GLSLM=NaN(1,length_lambda);
    parfor i=1:length_lambda
        [~, p_val_F, p_val_chi]= FAR_CLR_JGLS_GLSLM_test_joint_rp( R, F, lambda_f(i), T1, m,GLSLM_id);              % test assets, , lambda_f(i));
        pvalueF_FAR(i)=p_val_F.FAR;
        pvalue_FAR(i)=p_val_chi.FAR;
        pvalueF_sCLR(i)=p_val_F.sCLR;
        pvalue_sCLR(i)=p_val_chi.sCLR;
        pvalueF_JGLS(i)=p_val_F.JGLS;
        pvalue_JGLS(i)=p_val_chi.JGLS;
        pvalueF_GLSLM(i)=p_val_F.GLSLM;
        pvalue_GLSLM(i)=p_val_chi.GLSLM;
    end
    % save P_value_k.mat pvalue_FAR pvalue_sCLR pvalue_JGLS pvalue_GLSLM pvalueF_FAR pvalueF_sCLR pvalueF_JGLS pvalueF_GLSLM  lambda_f
end

%% sFAR distribution simulation data
if data_id==6
    clear;
    rng(1);
    data=xlsread('data_KZ2018.xls');    % load data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This code is used to simulate the distribution of the subset FAR
    % statistic and to compare with the the JGLS statistic, and compare with the
    % two bounded distributions
    %   (1) sFAR
    %   (2) JGLS        the normal JGLS statisic
    %   (3) bound_f     the bounds using the finite sample distribution of the JGLS
    %   (4) bound_a     the bounds using the asymptotic chi2 distribution
    
    % we copnsider a multi-factor simulation: using three FF factors and
    % a consumption series, where we only test the risk premium with respect to
    % the consumption series (weak factor)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    G=data(:,4);                        % consumption growth
    FF=data(:,120:122);               % FF three factors
    G=[G, FF];   % the four factors we gonna use
    K=size(G,2);
    Vff=cov(G);
    %  Tuning parameters
    lambda_f =2*ones(K,1);
    T=500;
    rdd_num=50000;
    asset_num_minus2 =43;
    
    lambda_a=2;
    a=4;
    b=(1:3);
    %   Simulated factors
    F1=mvnrnd(zeros(1,K),Vff, T);
    F=F1-mean(F1,1);
    F_pluslambda= F+ ones(T,1)*lambda_f';
    F_pluslambda_2=(F_pluslambda'* F_pluslambda);
    F_tilde=F_pluslambda/F_pluslambda_2;
    F_hat= F/(F'* F);
    Q=inv((F'* F)/T);
    
    n=asset_num_minus2;
    % step one: caliberating from the empirical data
    R=[data(:,9:9+n) data(:,7)];         % test assets
    N=size(R,2);
    K=size(G,2);
    V=cov([G R]);
    Vc=V(1:K,1:K);                      % V_f
    Vr=V(K+1:end,K+1:end);                  % V_r
    Vcr=V(1:K,K+1:end);                   % V_fr
    Vrc=Vcr';                           % V_rf
    Qff=V(1:K,1:K);                      % V_f
    Sigma=V(K+1:end,K+1:end);                  % V_r
    beta=Vrc/Vc;
    Omega=100*(Vr-beta*Vc*beta');
    
    sFAR = NaN(rdd_num,1);
    JGLS = NaN(rdd_num,1);
    FAR = NaN(rdd_num,1);
    F_dist=NaN(rdd_num,1);
    F2_dist=NaN(rdd_num,1);
    Chi2_dist=NaN(rdd_num,1);
    Chi22_dist=NaN(rdd_num,1);
    C_JGLS=(T-N+1)/(N-K-1)/(T-K-1);
    C_FAR=(T-K-N+1)/(N-1)/(T-K-1);
    F_quantile(n)= finv(0.95, N-K-1,T-N+1)/C_JGLS;
    parfor i=1:rdd_num
        v=mvnrnd(zeros(1,N),Omega,T);
        R=ones(T,N)*lambda0+(F+kron(ones(T,1),lambda_f'))*beta'+v;
        Jn=[eye(N-1), -ones(N-1,1)];
        R=R*Jn'; % now  N-1 cross section
        B_tilde=R'*F_tilde;
        B_hat=  R'*F_hat;
        Sigma_hat=(R-F*B_hat')'*(R-F*B_hat')/(T-K); % this might be the poblem??
        Rmean_betalambda=mean(R)'-B_tilde*lambda_f;
        GLSLM=Rmean_betalambda'/ Sigma_hat* B_tilde/ (B_tilde'/Sigma_hat*B_tilde) *B_tilde'/Sigma_hat*Rmean_betalambda;
        GLS_LM= T*GLSLM/(1-lambda_f'/(F_pluslambda_2/T)*lambda_f);
        FAR(i) = T* (Rmean_betalambda'/ Sigma_hat*Rmean_betalambda) /(1-lambda_f'/(F_pluslambda_2/T)*lambda_f);
        JGLS(i,:)=FAR(i)- GLS_LM ;
        % sFAR statistic
        A1=[mean(R)'-B_hat(:,a)*lambda_a  ,B_hat(:,b)]'/Sigma_hat*[mean(R)'-B_hat(:,a)*lambda_a  ,B_hat(:,b)];
        A2= ([1+lambda_a'*Q(a,a)*lambda_a, -lambda_a'* Q(a,b); Q(b,a)*lambda_a, Q(b,b)]);
        A= A1/A2;
        [V,D] = eig(A);
        D=real(D); % may cantain small complex part due to calculation issues
        [d,ind] = sort(diag(D));
        sFAR(i,:) = T*d(1);
        
        F_dist(i,:)=frnd(N-K-1,T-N+1);
        F2_dist(i,:)=frnd(N-1,T-N+1-K);
        Chi2_dist(i,:)=chi2rnd(N-K-1);
        Chi22_dist(i,:)=chi2rnd(N-1);
    end
    %     save sFAR_dist_comparison_T500.mat sFAR FAR F_dist F2_dist JGLS Chi2_dist Chi22_dist C_JGLS C_FAR
end






%% Draw Figures %%%%%%%%%%%%%%%%%
%% Draw power curves:
if fig_id==1
    %% 3*1 subplots for fixed T,r, variopus N
    % load Power_acrossNT100_N_3_step20_43
    load  Power_acrossNT500_N_3_step20_43_fix_r_15_new2
    x=power_disturb;
    a=(1:11);
    figure()
    subplot(3,1,1)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(3,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(3,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(3,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(3,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(3,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(3,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(3,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(3,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=5,T=100')
    subplot(3,1,2)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(23,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(23,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(23,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(23,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(23,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(23,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(23,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(23,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=25,T=100')
    subplot(3,1,3)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(43,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(43,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(43,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(43,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(43,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(43,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(43,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(43,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=45,T=100')
    %% 3*3 subplots for fixed T, 3 various Ns and rs
    clear
    load Power_acrossNT100_N_3_step20_43_fix_r
    x=power_disturb;
    length_x=length(x);
    a=(21:length_x-20);
    x=x(a);
    figure
    subplot(3,3,1)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(3,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(3,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(3,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(3,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(3,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(3,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(3,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(3,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=5,T=100, r=1.5')
    subplot(3,3,4)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(23,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(23,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(23,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(23,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(23,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(23,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(23,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(23,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=25,T=100, r=1.5')
    subplot(3,3,7)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(43,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(43,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(43,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(43,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(43,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(43,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(43,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(43,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=45,T=100, r=1.5')
    load Power_acrossNT500_N_3_step20_43_fix_r_15_new2
    x=power_disturb;
    length_x=length(x);
    a=(11:length_x-10);
    x=x(a);
    subplot(3,3,2)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(3,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(3,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(3,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(3,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(3,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(3,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(3,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(3,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=5,T=100, r=15')
    subplot(3,3,5)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(23,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(23,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(23,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(23,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(23,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(23,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(23,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(23,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=25,T=100, r=15')
    subplot(3,3,8)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(43,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(43,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(43,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(43,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(43,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(43,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(43,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(43,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=45,T=100, r=15')
    load Power_acrossNT500_N_3_step20_43_fix_r_150_new2
    x=power_disturb;
    length_x=length(x);
    a=(11:length_x-10);
    x=x(a);
    subplot(3,3,3)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(3,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(3,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(3,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(3,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(3,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(3,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(3,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(3,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=5,T=100, r=150')
    subplot(3,3,6)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(23,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(23,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(23,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(23,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(23,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(23,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(23,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(23,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=25,T=100, r=150')
    subplot(3,3,9)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(43,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(43,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(43,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(43,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(43,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(43,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(43,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(43,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=45,T=100, r=150')
end

%% Draw size curves (across N)
if fig_id==2
    %% 1*1 size curves across N for fixed r fixed T
    load size_acrossNT100_fizxedr_new2
    %   load size_acrossNT100_new
    figure()
    hold on
    plot(reject_frequency_CLR(:,2),'k-.','LineWidth',1.5)
    plot(reject_frequency_CLR(:,1),'r-.','LineWidth',1.5)
    plot(reject_frequency_FAR(:,1),'k-','LineWidth',1.5)
    plot(reject_frequency_FAR(:,2),'r-','LineWidth',1.5)
    plot(reject_frequency_GLSLM(:,1),'k:','LineWidth',1.5)
    plot(reject_frequency_GLSLM(:,2),'r:','LineWidth',1.5)
    plot(reject_frequency_JGLS(:,1),'k--','LineWidth',1.5)
    plot(reject_frequency_JGLS(:,2),'r--','LineWidth',1.5)
    hold off
end

%% Draw size curves (across lambda_f)
if fig_id==3
    clear
    load Size_acrossNT100_N_3_step20_43_fix_r_15
    x=power_disturb;
    a=(1:41);
    figure
    subplot(3,1,1)   % add first plot in 2 x 1 grid
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(3,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(3,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(3,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(3,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(3,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(3,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(3,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(3,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=5,T=100')
    subplot(3,1,2)   % add first plot in 2 x 1 grid
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(23,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(23,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(23,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(23,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(23,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(23,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(23,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(23,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=25,T=100')
    subplot(3,1,3)   % add first plot in 2 x 1 grid
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(reject_frequency_FAR(43,1,a), 1, [])','k-')
    plot(x, reshape(reject_frequency_FAR(43,2,a), 1, [])','r-')
    plot(x, reshape(reject_frequency_JGLS(43,1,a), 1, [])','k--')
    plot(x, reshape(reject_frequency_JGLS(43,2,a), 1, [])','r--')
    plot(x, reshape(reject_frequency_GLSLM(43,1,a), 1, [])','k:')
    plot(x, reshape(reject_frequency_GLSLM(43,2,a), 1, [])','r:')
    plot(x, reshape(reject_frequency_CLR(43,1,a), 1, [])','r-.')
    plot(x, reshape(reject_frequency_CLR(43,2,a), 1, [])','k-.')
    hold off
    ylim([0 1])
    title('N=45,T=100')
end

%% Draw sFAR P-value curve:
if fig_id==4
    clear
    load P_value_LL_new.mat
    x=lambda_f;
    figure();
    subplot(1,3,1)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(pvalue_sFAR(1,1,:), 1, [])','k-')
    plot(x, reshape(pvalue_sFAR(1,2,:), 1, [])','r-')
    hold off
    ylim([0 0.4])
    title('cay')
    subplot(1,3,2)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(pvalue_sFAR(2,1,:), 1, [])','k-')
    plot(x, reshape(pvalue_sFAR(2,2,:), 1, [])','r-')
    hold off
    ylim([0 0.4])
    title('\Delta Concumption')
    subplot(1,3,3)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, reshape(pvalue_sFAR(3,1,:), 1, [])','k-')
    plot(x, reshape(pvalue_sFAR(3,2,:), 1, [])','r-')
    hold off
    ylim([0 0.4])
    title('\Delta Concumption \times cay')
end

%% Draw single weak factor P-value curve:
if fig_id==5
    
    clear;
    load P_value_k.mat;
    
    figure();
    subplot(2,2,1)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, pvalueF_FAR ,'k-','LineWidth',1.5)
    plot(x, pvalue_FAR,'r-','LineWidth',1.5)
    hold off
    ylim([0 1])
    title('FAR')
    subplot(2,2,2)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, pvalueF_sCLR ,'k-.','LineWidth',1.5)
    plot(x, pvalue_sCLR,'r-.','LineWidth',1.5)
    hold off
    ylim([0 1])
    title('sCLR')
    subplot(2,2,3)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, pvalueF_JGLS ,'k--','LineWidth',1.5)
    plot(x, pvalue_JGLS,'r--','LineWidth',1.5)
    hold off
    ylim([0 1])
    title('JGLS')
    subplot(2,2,4)
    hold on
    plot(x, 0.05*(1+0*x),'k')
    plot(x, pvalueF_GLSLM,'k:','LineWidth',1.5)
    plot(x, pvalue_GLSLM,'r:','LineWidth',1.5)
    hold off
    ylim([0 1])
    title('GLSLM')
end

%% Draw sFAR distribution
if fig_id==6
    
    %% PDF for test statistics
    figure();
    histogram(sFAR,'FaceColor','r','Normalization','probability','FaceAlpha',1);
    hold on;
    histogram(F_dist/C_JGLS,'FaceColor','b','Normalization','probability','FaceAlpha',0.5);
    % histogram(F2_dist/C_FAR,'FaceColor','b','Normalization','probability','FaceAlpha',0.5);
    histogram(JGLS,'FaceColor','y','Normalization','probability');
    histogram(Chi2_dist,'FaceColor','g','Normalization','probability','FaceAlpha',0.5);
    % histogram(Chi22_dist,'FaceColor','g','Normalization','probability','FaceAlpha',0.5);
    hold off;
    xlabel(''); ylabel('Densities');
    % title(['Densities of simulated \min FAR, JGLS and randoms draws from \chi^2(N-K-1), \mathcal{F}(N-K-1,T-N+1)/C_{JGLS} given'  '; N=' num2str(N) '; T= ' num2str(T)])
    % clear
    % load minFAR_dist_comparison_T500.mat
    figure();
    [f,x] = ecdf(sFAR);
    plot(x,f,'r-');
    xlim([0,200]);
    hold on;
    [f,x] = ecdf(F_dist/C_JGLS);
    plot(x,f,'b--');
    % [f,x] = ecdf(F2_dist/C_FAR);
    % plot(x,f,'b--');
    [f,x] = ecdf(JGLS);
    plot(x,f,'y--');
    % [f,x] = ecdf(FAR);
    % plot(x,f,'y--');
    [f,x] = ecdf(Chi2_dist);
    plot(x,f,'g--');
    % [f,x] = ecdf(Chi22_dist);
    % plot(x,f,'g--');
    hold off;
    xlabel(''); ylabel('Empirical CDF');
end

%% Draw qq-plots for the return data
if fig_id==7
    % Purpose: This programs draw the Q-Q plots of returns versus standard normal distributions.
    % and do a one-sample Kolmogorov-Smirnov (normality) test for these
    % returns.
    % Data: daily returns (6 Portfolios Formed on Size and Book-to-Market (2 x 3) [Daily]),
    % rest returns are the FF benchmark portfolio returns from
    % 1986jan02-2008Jan02
    clear
    % daily
    load daily_return.mat    % load data
    data_d= table2array(daily_return);
    figure()
    qqplot(zscore(data_d(:,3)))
    kstest(zscore(data_d(:,3)))
    figure()
    qqplot(zscore(data_d(:,4)))
    kstest(zscore(data_d(:,4)))
    figure()
    qqplot(zscore(data_d(:,5)))
    kstest(zscore(data_d(:,5)))
    
    % monthly
    load monthly_return.mat    % load data
    data_m= table2array(Monthly_return);
    
    
    figure()
    qqplot(zscore(data_m(:,6)))
    kstest(zscore(data_m(:,6)))
    figure()
    qqplot(zscore(data_m(:,4)))
    kstest(zscore(data_m(:,4)))
    figure()
    qqplot(zscore(data_m(:,5)))
    kstest(zscore(data_m(:,5)))
    
    
    % quarterly
    load quarterly_return.mat    % load data
    data_q= table2array(quarterly_return);
    
    figure()
    qqplot(zscore(data_q(:,6)))
    kstest(zscore(data_q(:,6)))
    figure()
    qqplot(zscore(data_q(:,4)))
    kstest(zscore(data_m(:,4)))
    figure()
    qqplot(zscore(data_q(:,6)))
    kstest(zscore(data_m(:,6)))
    
    % annual
    load  annual_return.mat    % load data
    data_a= table2array(annual_return);
    
    figure()
    qqplot(zscore(data_a(:,3)))
    figure()
    qqplot(zscore(data_a(:,4)))
    figure()
    qqplot(zscore(data_a(:,6)))
 
end



