function [CS] = ConsumerSurplus_choice(Xa,bhat,ihess, EstimOpt,MP)

% Function calculates Consumer surplus after estimating MNL model
% Returns matrix CS
    % For each of individual we calculate
        % Estimate based on simulation (first row)
        % Standard error (second row)
        % Confidence interval (thrid and fourth row)
        % Estimate calculated directly (without simulation)

rng(100001);
tic

NAlt = EstimOpt.NAlt;
NP = EstimOpt.NP;
n = 10000; % no of. draws for se

% MP is a multiplier (e.g. if cost is divided in the model by 100, then MP should be 100)
if isempty(MP)
    MP = 1;
end

MP = -MP*sign(Xa(2,end)); % If cost enters with negative sign then this term is positive
bx = mvnrnd(bhat', ihess,n); 
CS = zeros(NP,4); % Estimate, Std, 95% CI, Estimate, no Sim

WTP = zeros(NP,n);
for i = 1:n
    fit = exp(Xa*bx(i,:)');
    fit = reshape(fit, NAlt, NP);
    WTP(:,i) = log(sum(fit,1)')/bx(i,end);
end
WTP = WTP*MP;
CS(:,1:4) = [mean(WTP,2), std(WTP,[],2), prctile(WTP, 2.5, 2), prctile(WTP, 97.5, 2)];
fit = exp(Xa*bhat);
fit = reshape(fit, NAlt, NP);
WTP2 = log(sum(fit,1)')/bhat(end);
CS(:,5) = MP*WTP2;

toc
        