%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 [res_SE, res_PE] = fn_ivtobit_Naive_PE_SE(theta,sig_U,h,Sigma_th_sig_U,bProb)
  % bProb = 1 => PE^{Pr}, otherwise PE^{E}
  if nargin<5; bProb=1; end
  
  theta = theta(:); h = h(:);
  dim_theta = length(theta);
  assert(all(size(Sigma_th_sig_U)==(dim_theta+1)*[1 1])); 
  assert(length(h)==dim_theta);
  
  theta_h = theta'*h;
  phi = normpdf(theta_h/sig_U);
  common_mul = phi/sig_U;

  res_SE = zeros(dim_theta, 1);
  res_PE = zeros(dim_theta, 1);
  
  for j=1:dim_theta
    e_j = zeros(dim_theta,1); 
    e_j(j) = 1;
    
    if bProb
      PE_j = common_mul*theta(j);
      grad_j_th = (-(theta_h*theta(j)/sig_U^2)*h + e_j)*common_mul;
      grad_j_sigU = (theta_h^2/sig_U^2-1)*phi*theta(j)/sig_U^2;
    else
      ncdf = normcdf(theta_h/sig_U);
      PE_j = ncdf*theta(j);
      grad_j_th = common_mul*theta(j)*h + e_j*ncdf;
      grad_j_sigU = phi*theta(j)*(-theta_h/sig_U^2);
    end
    grad_j = [grad_j_th; grad_j_sigU];
    res_SE(j) = sqrt(grad_j'*Sigma_th_sig_U*grad_j);
    res_PE(j) = PE_j;
  end
  
end

