%% Welfare in the multinomial logit model
% By Tobias Klein, May 2026
% We consider consumers with utilities u_{ij}=\beta_0 + \beta_1 p_{ij} + \varepsilon_{it}. Prices p_{ij} are individual-specific. They are drawn from a joint log-normal distribution with parameters \mu_p and \Sigma_p. \varepsilon_{ij} is type 1 extreme value.

% clean up and set seed
clear all
rng(17);

% parameters for data generating process
N=100000; % number consumers
beta=[-2;-0.1]; % parameters of the utility function
mu_p = [2.9 3 3]; % means of prices
sig_p = [.1 0  0;
         0 .2 .1
         0 .1 .2]; % variance-covariance for prices

% We generate data.

p = exp(mvnrnd(mu_p,sig_p,N)); %prices
X1=[ones(N,1) p(:,1)]; %matrix with explanatory variables for first alternative
X2=[ones(N,1) p(:,2)];
X3=[ones(N,1) p(:,3)];
epsilon1=-evrnd(0,1,N,1); %taste shock for first alternative
epsilon2=-evrnd(0,1,N,1);
epsilon3=-evrnd(0,1,N,1);
u1=X1*beta+epsilon1; %utility for first alternative
u2=X2*beta+epsilon2;
u3=X3*beta+epsilon3;
[~,y123] = max([u1 u2 u3],[],2); %choice
tabulate(y123)

% Now we investigate the effect of dropping alternative 2 from the choice set.

y13 = ones(size(u1));
y13(u3 > u1)  = 3;
tabulate(y13)

% cross tabulation
crosstab(y123,y13)

% maximal utilities
maxU123 = max([u1,u2,u3],[],2);
maxU13 = max([u1,u3],[],2);
[y123,maxU123,y13,maxU13,maxU13-maxU123]

% This table shows the optimal choice from the full choice set, the associated maximal utility, the choice from the smaller choice set without alternative 2, the associated maximal utility for that smaller choice set, and the utility loss associated with not being able to choose alternative 2 anymore.
% We can also translate this into monetary units and calculate the compensating variation. For this, we use the price coefficient.

moneyValueOneUtil = -1/beta(2)
utilityLoss = maxU13-maxU123;
compensatingVariation = -moneyValueOneUtil*utilityLoss
hist(compensatingVariation)
mean(compensatingVariation)

% From an ex ante perspective, realizations of \varepsilon_{ij} are not known. Based on the result \mathbb{E}[\max u_{ij}]=\gamma + \log \left(\sum_{j \in C} \exp (u_{ij})\right) we can calculate the expected utility loss for each consumer:

expMeanUtility = exp([X1*beta,X2*beta,X3*beta]);
eulerGamma = -psi(1);  % Euler–Mascheroni constant
expMaxUtility123 = eulerGamma + log(sum(expMeanUtility,2));
expMaxUtility13 = eulerGamma + log(expMeanUtility(:,1)+expMeanUtility(:,3));
[maxU123,expMaxUtility123,maxU13,expMaxUtility13]

% We see here realizations of utility maxima (column 1 and 3) and expected maximal utilities (column 2 and 4). These are of course not the same. However, they will be the same on average across consumers:

mean([maxU123,expMaxUtility123,maxU13,expMaxUtility13],1)

% From this, we can again calculate the expected compensating variation.

expCompensatingVariation = -moneyValueOneUtil*mean(expMaxUtility13-expMaxUtility123)