function pv = pvalue_LStat(L_Stat,k)
    % Can be sped up 
    T = 1000;
    nrep = 5000;
    tstat = NaN(nrep,1);
    for i = 1:nrep
      X = randn(T,k);
      Xm = X - mean(X);
      cXm = cumsum(flipud(Xm));
      Lm_mat = sum(cXm.^2)/T^2;
      tstat(i) = sum(Lm_mat);
    end
    pv = mean(tstat>L_Stat);
end