%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code is used to generate critical values for the GLS-LM statistic \
% for various N with fixed T = 50 100 200 2000
% the emprirical distribution of the test statistic under the null
% with various N and T
% Q1: can I find the ounds of such test statistics??

% for the simulation, we only need to know Q_ff and lmabda_F and Sigma
% and we start with a single factor setting
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
rng(1); 
% Input values N varies from 3 to asset_num
% T (sample zie),K (factor number)  is fixed
T=50;
K=1;
asset_num =48;
rdd_num=500000;



GLSLM_qunatile=NaN(asset_num,1);
for n=3:asset_num
    fprintf('%2.0f out of %2.0f\n', n, asset_num);
    N=n;
    GLS_LM=NaN(rdd_num,1);
    
    parfor i=1:rdd_num
        psi=mvnrnd(zeros(N-1,1),eye(N-1),1);
        WZ_i =mvnrnd(zeros(1,N-1), eye(N-1), T-K-1);
        W=WZ_i'*WZ_i/(T-K-1);
        PSI=reshape(psi',N-1,[]);
        C=eye((N-1),(N-1-K));
        GLS_LM(i)= PSI'/W*PSI -  PSI'*C/(C'*W* C)*C'*PSI;
    end
    GLSLM_qunatile(n)= quantile(GLS_LM,0.95);

end
    num=2+(1:asset_num)';
    save GLSLM_qunatile_T50.mat num GLSLM_qunatile
    GLSLM_qunatile=GLSLM_qunatile(1:asset_num);
    figure()
    plot(GLSLM_qunatile(1:asset_num))
