%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_ivtobit_grad_new(Y,X,W,Z, theta,gam,piZW,sUgV,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(sUgV) && isscalar(sV))
  assert(sUgV>0 && 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;
  IYp = Y>0;
  neg_IMR = -inv_mills_ratio(m_i/sUgV);
  grad_m_i = (1-IYp).*neg_IMR/sUgV - IYp.*(m_i-Y)/(sUgV^2);
  grad_sUgV = -(1-IYp).*m_i.*neg_IMR/(sUgV^2) - IYp.*(1/sUgV-(m_i-Y).^2/(sUgV^3));
  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_sUgV 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
