%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 [b_hat, asyV_textbook, loglik] = probit_MLE(Y, X)
  n = size(Y,1);  assert(all(size(Y)==[n,1]));
  assert(size(X,1)==n);
  b0 = 2.5*(X\Y); %b_ols ~ 0.4 b_probit
  ix_intercept = var(X)==0;
  if sum(ix_intercept)==1 % found intercept
    b0(ix_intercept) = b0(ix_intercept)-1.25; %b_ols ~ 0.4 b_probit + 0.5 for the intercept
  end
  options = optimoptions(@fminunc,'Algorithm', 'quasi-newton' ... 'trust-region' ... 'quasi-newton' ... 
      , 'SpecifyObjectiveGradient', true ...  %, 'HessianFcn','objective' ... %,'GradObj','on' ...
      , 'OptimalityTolerance', 1e-9, 'StepTolerance', 1e-9, 'MaxFunctionEvaluations', 1e3*(length(b0)+1) ...
      , 'Display', 'off'); % ... 'final-detailed' ...
    
  Z = X.*(2*Y-1); %bsxfun(@times, X, (2*Y-1));
  [b_hat,loglik,exitflag,output] = fminunc(@negative_probit_loglik, b0, options); %#ok<ASGLU>
    
  if nargout>=2
    [~, I] = negative_probit_loglik(b_hat,1);
    asyV_textbook = inv(I);
  end

function [loglik, grad, H] = negative_probit_loglik(b,bReturnScoreVariance)
  if nargin<2; bReturnScoreVariance=0; end %if true, returns covariance of the score in 'grad' output
  Zb = full(Z*b(:));
  loglik = -full(sum(log_normcdf(Zb)));
  
  if nargout>=2
    phi_Phi_Z = Z.*inv_mills_ratio(-Zb);
    
    if bReturnScoreVariance
      grad = full(cov(phi_Phi_Z));
    else
      grad = -full(sum(phi_Phi_Z));
    end
    if nargout>=3
      H = +full( phi_Phi_Z'*phi_Phi_Z + phi_Phi_Z'*(bsxfun(@times, Z, Zb)) );
    end
  end
end

end %probit_MLE()

function s = log_normcdf(s)
  % computes log(normcdf(s))
  todo = s<0; %for these we will replace 1-normcdf(s) with erfcx...
  s(~todo) = log(erfc(-s(~todo)/sqrt(2)))-log(2);
  s(todo) = log(erfcx(-s(todo)/sqrt(2)))-(s(todo).^2/2)-log(2);
end

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