%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 [CI_max, cv] = inference_on_max(X, Sigma, alpha)
  % constructs a 1-sided CI for the max(mu_1,mu_2)
  % Inputs: 
  %   X is either a 2-element vector or an Kx2 matrix.
  %   Sigma is a 2x2 cov matrix
  %
  % Output: 
  %   CI_max is Kx1 matrix of the CIs
  assert(all(size(Sigma)==[2 2]));
  if numel(X)==2
    X = X(:)'; %row vector 1x2
  else
    assert(size(X,2)==2); % Kx2 matrix
  end
  assert(isscalar(alpha));
  assert(0<alpha && alpha<1/2);
  arr_sig = sqrt(diag(Sigma));

  D = load("stored_grid_CV_for_max_1-sided.mat"); %disp(D);
  [d_q, i_q] = min(abs(D.arr_alpha-alpha));
  if d_q>1e-7
    warning('myFun:inference_on_max__alpha', 'inference_on_max: alpha %7.5f not available', alpha);
  end
  
  % handle numerically slightly asymmetric matrices
  rho = [Sigma(1,2) Sigma(2,1)]/prod(arr_sig);
  %assert(abs(rho(1)-rho(2))<1e-10); %issymmetric(Sigma));
  if abs(rho(1)-rho(2))>1e-10
    warning('myFun:inference_on_max:asymSig', 'inference_on_max: Sigma is not symmetric?');
  end
  rho = (rho(1)+rho(2))/2;
  
  i_rho = find(D.arr_rho<=min(rho,1),1,'last');
  
  cv = D.arr_Q(i_rho, i_q);
  % cv = interp1(D.arr_rho, D.arr_Q(:, i_q), rho, 'spline');

  CI_max = max(X-arr_sig'*cv,[],2);
end

