%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Part of the replication package for the paper
%   "Marginal Effects for Probit and Tobit with Endogeneity"
%   by Kirill S. Evdokimov, Ilze Kalnina, and Andrei Zeleneev.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function score = fn_ivprobit_grad_new(Y,X,W,Z, theta,gam,piZW,sV)
  Y = Y(:); theta = theta(:); gam = gam(:); piZW = piZW(:);
  n = length(Y);
  dim_W = size(W,2);
  assert(all(size(X)==[n 1]));
  assert(all(size(Z)==[n 1]));
  assert(isscalar(gam) && isscalar(sV))
  assert(sV>0);
  assert(length(theta)==dim_W+1);
  assert(length(piZW)==dim_W+1);
  
  X_minus_ZW_pi = X - [Z W]*piZW;
  m_i = [X W]*theta + gam*X_minus_ZW_pi;
  Y21 = 2*Y-1;
  IMR = inv_mills_ratio(-Y21.*m_i);
  grad_m_i = Y21.*IMR;
  grad_sV = -1/sV+(X_minus_ZW_pi).^2/(sV^3);

  score = [ grad_m_i.*[X W X_minus_ZW_pi] (-grad_m_i*gam + X_minus_ZW_pi/sV^2).*[Z W] grad_sV];
end

function s = inv_mills_ratio(s)
  % computes  normpdf(s)./(1-normcdf(s))
  s = sqrt(2/pi)./erfcx(s/sqrt(2)); % works well for all s, including very large
end
