function [L] = logl_ss(b)
% The following is the correct loglikelihood function for a probit model
% with sample selectivity.

global x w y z N

b0 = b(1);
b1 = b(2);
sigmau = b(3);
g0 = b(4);
g1 = b(5);
rho = b(6);

XB = b0*ones(size(x,1),1) + b1*x ;
WG = g0*ones(size(w,1),1) + g1*w ;

L=0 ;

for i = 1:N
  if z(i) == 1
    L = L + log((1/sigmau)*normpdf((y(i) - XB(i))/sigmau))... 
            + log(normcdf((1/sqrt(abs(1-rho^2)))*(WG(i) + (rho/sigmau)*(y(i) - XB(i)))));
  elseif z(i) == 0
    L = L + log(normcdf(-WG(i)));
  end   
end

L = -L;

end

