function [betamat,Emat]= gs_select(X,k1,k,m,n,r1,r2,mu0,V0,...
betadraw,v0,S0,Edraw,f,g,nf,ng,l0,r0,lpos,rpos,ydraw)
R = r1+r2;
betamat=zeros(k,r2);
Emat=zeros((m*(m+1)/2)-1,r2); %collects unique elements of E as a vector
% recall first element is one, so we can skip a row;
ff=(1:2:(2*n-1))'; % indexes odd rows only, i.e. X1 stuff, needed for draws of ydraw below
gg=(2:2:2*n)'; % indexes even rows only, i.e. X2 stuff, needed for draws of ydraw below
thoucount=1000;
% start main loop
i=1;
for i=1:R
%%%%%%%%%%%%%%%%%%%%%%
% draw beta
%%%%%%%%%%%%%%%%%%%%%
int1=0;
int2=0;
j=1;
for j=1:n
Xi=X{j}; % note: much faster than cell2mat!
int1=int1+(Xi'/Edraw)*Xi;
int2=int2+(Xi'/Edraw)*ydraw(j,:)';
end
V1=inv(inv(V0)+int1);
mu1=V1*(V0\mu0+int2);
betadraw=mvnrnd(mu1,V1)';
if i>r1
betamat(:,i-r1)=betadraw;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v1=v0+n;
yr= reshape(ydraw',2*n,1);%stacking of yi pairs - 2xn by 1
Xr=cell2mat(X); %2xn by k
int3=yr-Xr*betadraw;
int3=reshape(int3,m,n)';%n x 2
S1=S0+(int3'*int3);
Edraw=iw_sig11(v1,inv(S1),1);
if i>r1
Emat(:,i-r1)=[Edraw(1,2);Edraw(2,2)];
end
%%%%%%%%%%%%%%%%%%%
% draw ystar
%%%%%%%%%%%%%%%%%%%
sig11=Edraw(1,1); %"1", of course
sig12=Edraw(2,1);
sig22=Edraw(2,2);
rho=sig12/(sqrt(sig11*sig22));
beta1=betadraw(1:k1);
beta2=betadraw(k1+1:end);
X1=Xr;
X2=Xr;
X1(gg,:)=[]; % keep odd rows only
X1(:,k1+1:end)=[]; % get rid of trailing zeros to get n by k1
X2(ff,:)=[]; % keep even rows only
X2(:,1:k1)=[]; % get rid of leading zeros to get n by k2
y2=ydraw(:,2); % define y1 later, when it's updated
%Case 1: y1i=1
%%%%%%%%%%%%%%%%%%%
% STEP 1: draw y1star for nonzero amounts, conditional on y1star>0, y2star
% (observed)
mu=X1*beta1+(sig12/sig22)*(y2-X2*beta2);%conditional mean
mu(f)=[]; % eliminate zero amounts, so now mu is length g
yint=tnormrnd(mu,sqrt(sig11*(1-rho^2)),lpos,rpos);
ydraw(g,1)=yint;
% STEP 2: draw both y1star and y2star for zero amounts
% step 2a: draw y1star from its truncated marginal distribution
mu=X1*beta1;
mu(g)=[]; %eliminate nonzero amounts
yint=tnormrnd(mu,sqrt(sig11),l0,r0);
ydraw(f,1)=yint;
% step 2b: draw y2star|y1star
y1=ydraw(:,1); %KEY: don't forget to up date y1!!
mu=X2*beta2+(sig12/sig11)*(y1-X1*beta1);
mu(g)=[]; % eliminate nonzero amounts
yint=normrnd(mu,sqrt(sig22*(1-rho^2))); % no truncation here -
%it's like the unobserved wage of a nonworking woman could be anything...
ydraw(f,2)=yint;
% so the only element of ydraw that's never updated is ydraw(g,2), i.e.
% the positive amounts
if i== thoucount
i
thoucount=thoucount+1000;
end
end