function [lfst,lfst_pvalue] = lfst_pv(X,g)
% compute LFST statistic and pvalue
% X  = cosine transforms
% g = value of g under alternative
q = size(X,1);
tmp = (1:1:q)';
tmp = (tmp*pi).^2;
cov_rw = 1./tmp;
cov_ha = ones(q,1)+(g*g*cov_rw);
cov_ha_inv = 1./cov_ha;
chol_ha = sqrt(cov_ha_inv);
Xa = X.*chol_ha;
lfst = sum(X.^2)/sum(Xa.^2);

% Generate 1000 draws under Ho for pvalue calculation
nsim = 50000;
Xsim = randn(q,nsim);
Xasim = Xsim.*repmat(chol_ha,1,nsim);
nx = sum(Xsim.^2);
dx = sum(Xasim.^2);
lfst_sim = (nx./dx)';
lfst_pvalue = mean(lfst_sim >= lfst);

end

