function [b,se,db,dse] = clsprobit(y,d,X,z)
%CLSPROBIT Estimates coefficients, marginal effects and 
%   s.e.s of the best probit least squares approximation to 
%   the LARF from Abadie (2003). Marginal effects are evaluated
%   at the mean of the covariates for compliers. CLSPROBIT uses 
%   a probit first step.
%
%   Usage:   [b,se,db,dse] = clsprobit(y,d,X,z) 
%
%   Input:   y   (nx1)  binary outcome variable     
%            d   (nx1)  binary treatment indicator 
%            X   (nxk)  covariates           
%            z   (nx1)  binary instrument 
%
%   Output:  b   ((k+1)x1)  estimated parameters
%            se  ((k+1)x1)  standard errors for b
%            db  ((k+1)x1)  estimated marginal effects 
%            dse ((k+1)x1)  standard errors for db
%
%            The first elements of b, se, db and dse 
%            correspond to the treatment indicator variable (d).

%Code last modified: Nov. 13, 1999. Written by Alberto Abadie. 

n = length(y);

gamma = probit(z,X); Phi = normcdf(X*gamma);

nc_ratio = zeros(n,1);
nc_ratio(d==0 & z==1) = 1./Phi(d==0 & z==1);
nc_ratio(d==1 & z==0) = 1./(1-Phi(d==1 & z==0));

kappa = 1 - nc_ratio;

b1 = probit(y,[d X]);

step = 1;
gtol = 0.0001;
dg = 1;


while dg > gtol
 b0 = b1;

 u = y - normcdf([d X]*b0);
 g = [d X]'*spdiags(kappa.*normpdf([d X]*b0),0,n,n);
 
 delta =  step*inv(g*g')*g*u;
 b1 = b0 + delta;
 dg = delta'*g*g'*delta;

end

b = b1;

u = y - normcdf([d X]*b);

dM_theta = ([d X]*b).*normpdf([d X]*b).*u + (normpdf([d X]*b).^2);
M_theta = [d X]'*spdiags(kappa.*dM_theta,0,n,n)*[d X];

derkappa = zeros(n,1);
derkappa(z==1 & d==0) =   1./(Phi(z==1 & d==0).^2);
derkappa(z==0 & d==1) = - 1./((1-Phi(z==0 & d==1)).^2);

M_gamma = - [d X]'*spdiags(normpdf([[d X]*b]).*u.*derkappa.*normpdf(X*gamma),0,n,n)*X;

lambda = zeros(n,1);
lambda(z==0) = - normpdf(X(z==0,:)*gamma)./(1-normcdf(X(z==0,:)*gamma));
lambda(z==1) =   normpdf(X(z==1,:)*gamma)./normcdf(X(z==1,:)*gamma);

H_gamma = X'*spdiags(lambda.^2,0,n,n)*X;

S_1 = [d X]'*spdiags((normpdf([d X]*b).*u.*kappa).^2,0,n,n)*[d X];
S_2 = M_gamma*inv(H_gamma)*M_gamma';
S_3 = M_gamma*inv(H_gamma)*X'*spdiags(lambda.*kappa.*(-u).*normpdf([d X]*b),0,n,n)*[d X];

S = S_1 + S_2 + S_3 + S_3';

V = inv(M_theta)*S*inv(M_theta);

se = sqrt(diag(V));

wbar = mean(spdiags(kappa,0,n,n)*[d X])/mean(kappa);

db =normpdf(wbar*b)*b;
G = normpdf(wbar*b)*(eye(size([d X],2)) - (wbar*b)*b*wbar);

dV = G*V*G';

dse = sqrt(diag(dV));

binary = (sum([d X]==0 | [d X]==1)==n)';

for i=1:length(binary)
   if binary(i)==1
      wbar1 = wbar; wbar1(i) = 1;
      wbar0 = wbar; wbar0(i) = 0;
      db(i) = normcdf(wbar1*b) - normcdf(wbar0*b);
      g = normpdf(wbar1*b)*wbar1 - normpdf(wbar0*b)*wbar0;
      dse(i) = sqrt(g*V*g');
   end
end


function b1 = probit(y,x)

n = length(y);

b1 = zeros(size(x,2),1);
step = 1;
gtol = 0.0001;
dg = 1;

lambda = zeros(n,1);

while dg > gtol
 b0 = b1;
 lambda(y==0) = - normpdf(x(y==0,:)*b0)./(1-normcdf(x(y==0,:)*b0));
 lambda(y==1) =   normpdf(x(y==1,:)*b0)./normcdf(x(y==1,:)*b0);
 
 g = x'*lambda;
 h = x'*spdiags(((-x*b0).*lambda) - lambda.^2,0,n,n)*x;
 
 delta = - step*inv(h)*g;
 b1 = b0 + delta;
 dg = g'*delta;
end

return;