%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 [G, Om, Sig, SigId, w] = gmm_get_G_Om_Sig(mom_fn, b0, d_eps)
  if nargin<3; d_eps = 1e-4; end
  w = mom_fn(b0);
  [obs,m] = size(w); 
  Om = gather(cov(w)); %iOm = inv(Om);

  dim_b = length(b0);
  if isscalar(d_eps)
    d_eps = d_eps*ones(dim_b,1);
  else
    d_eps = d_eps(:);
    assert(length(d_eps)==dim_b);
  end
  
  G = zeros(m, dim_b);
  for j=1:dim_b
    e_d = zeros(size(b0));
    e_d(j) = d_eps(j);
    G(:,j) = gather(mean(mom_fn(b0+e_d)-mom_fn(b0-e_d)))'/(2*d_eps(j));
  end
  
  if nargout>=3
    if dim_b==m % handle just-identified case separately
      % iG = inv(G);
      % Sig = iG*Om*(iG');
      Sig = G\Om/(G');
    else
      Sig = inv(G'*(Om\G)); % the inverse could be a problem, so don't compute Sig, unless the user asks to do it
    end
  end
end
