function psi = psi_sim(p,Eps)
%--------------------------------------------------------------------------
% psi_sim.m: Calculate psi(p) from  Fosgerau, Melo, Shum, Sørensen (2021),
%            Equation (16), p. 4, up to simulation error.
%--------------------------------------------------------------------------
%
% DESCRIPTION: Calculate a solution psi_S(p) to the problem
%           minimize{exp(W_S(v)) - v'*p subject to v in R^J}
% where W_S is the empirical approximation 
%           W_S(v):=(1/S)sum_{s=S}^{S} max_j{v_j+eps_{sj}}
% to the surplus W(v):=E[max_j{v_j+eps_j}]
% based on S independent draws from the joint distribution of
% eps=(eps_1,...,eps_J).
%
% INPUT ARGUMENTS:
% p:        J x 1 vector of choice probabilities (p_j>0 all j)
% Eps:      J x S matrix of simulation draws
%
% OUTPUT ARGUMENT: 
% psi: The scalar psi_S(p)
%
% AUTHOR: The MATLAB code was written by Jesper Riis-Vestergaard Sørensen, 
% Department of Economics, University of Copenhagen, Denmark.
%   e:  jrvs@econ.ku.dk
%   w:  https://sites.google.com/site/jesperrvs
%
% DATE: 26 May 2022
%
% REFERENCE: When using this software, please cite
%   Fosgerau, Melo, Shum, Sørensen (2021) "Some remarks on CCP-based
%   estimators of dynamic models"
%   Economics Letters, Vol. 204, July 2021.
%   https://doi.org/10.1016/j.econlet.2021.109911
if nargin < 2
   error('More input arguments needed.');
end
J=size(p,1);
if nargout==1
    W_S=@(v) mean(max(v+Eps));                      % approximate surplus
    f=@(v) exp(W_S(v))-v'*p;                        % criterion
    psi0=zeros(J,1);                                % starting value
    opts=optimoptions(@fminunc,'Display','off');    % turn off display
    psi=fminunc(f,psi0,opts);                       % optimize
else
    error('Too many output arguments');
end
end
 
% EXAMPLES: For each of the examples below let
%   J=5;                                % #alternatives
%   p=(1:J)'/sum(1:J);                  % choice probabilities
%   S=1000000;                          % #simulation draws

% Independent Gumbel(0,1)/Type 1 Extreme Value
%   Eps=-evrnd(0,1,J,S);                % J x S indep. Gumbel(0,1) draws
%   psi=psi_sim(p,Eps);
%
% Independent N(0,1)
%   Eps=randn(J,S);                     % J x S indep. N(0,1) draws
%   psi=psi_sim(p,Eps);
%
% Correlated Normal
%   rho=.5;
%   Sigma=nan(J,J);
%   for j=1:J
%     for k=1:J       
%         Sigma(j,k)=rho^(abs(j-k));    % Toeplitz covariance matrix
%     end
%   end
%   Eps=mvnrnd(zeros(J,1),Sigma,S)';    % N(0_Jx1,Sigma) draws
%   psi=psi_sim(p,Eps);