function k = find_kalpha_20Dec2020(n,mbar,Ztilde1,Ztilde2,X1,D,W_order_B,B,Omega_hat,delta_1,delta_2,gamma,p_eta,alpha1)
 
%Create vectors for storing
mbar_star = zeros(size(mbar,2),B); %these would be m x j x n if Gamma_LRR(beta) was non-singleton. (See Sep 25 code)
Qstat = zeros(B,1);

for b = 1:B
%Redrawn covariates
X1b = X1(1,W_order_B(:,b));
Db = D(W_order_B(:,b));
     
%Indicator functions for X1, D.
X1_eq_1_D_eq_1 = (X1b'==1 & Db==1);
X1_eq_1_D_eq_0 = (X1b'==1 & Db==0);
X1_eq_0_D_eq_1 = (X1b'==0 & Db==1);
X1_eq_0_D_eq_0 = (X1b'==0 & Db==0);

%Redrawn Ztilde1, Ztilde2
Z1t = Ztilde1(W_order_B(:,b));
Z2t = Ztilde2(W_order_B(:,b));
  
 %Moment conditions for those parameters
m1 = (Z1t.*X1_eq_1_D_eq_1 - (p_eta*(delta_1+delta_2+gamma)+ (1-p_eta)*(log(delta_1+delta_2 + gamma))).*X1_eq_1_D_eq_1);
m2 = (Z1t.*X1_eq_1_D_eq_0 - (p_eta*(delta_1+delta_2) + (1-p_eta)*log(delta_1 + delta_2)).*X1_eq_1_D_eq_0);
m3 = (Z1t.*X1_eq_0_D_eq_1 - (p_eta*(delta_1+gamma) + (1-p_eta)*log(delta_1 + gamma)).*X1_eq_0_D_eq_1);
m4 = (Z1t.*X1_eq_0_D_eq_0 - (p_eta*(delta_1) + (1-p_eta)*log(delta_1)).*X1_eq_0_D_eq_0);

m5 = (p_eta*(delta_1+delta_2+gamma) + (1-p_eta)*log(delta_1+delta_2+gamma)).*X1_eq_1_D_eq_1-Z2t.*X1_eq_1_D_eq_1;
m6 = (p_eta*(delta_1+delta_2) + (1-p_eta)*log(delta_1+delta_2)).*X1_eq_1_D_eq_0-Z2t.*X1_eq_1_D_eq_0;
m7 = (p_eta*(delta_1+gamma) + (1-p_eta)*log(delta_1+gamma)).*X1_eq_0_D_eq_1-Z2t.*X1_eq_0_D_eq_1;
m8 = (p_eta*(delta_1) + (1-p_eta)*log(delta_1)).*X1_eq_0_D_eq_0-Z2t.*X1_eq_0_D_eq_0;

  
%mbar_star average bootstrap moments
mbar_star(:,b) = [mean(m1),mean(m2),mean(m3),mean(m4),mean(m5),mean(m6),mean(m7),mean(m8)];

%Test Statistic
Qstat(b) = min(sqrt(n)*(mbar_star(:,b)'-mbar)./sqrt(Omega_hat)); 
end

%sort stats
T_sorted = sort(Qstat); 
k = T_sorted(floor(alpha1*B)); %find the largest stat that is smaller than percentile.

end

