%% Name of file: "Incognito_code_W_v1.m"
 

% This Matlab code generates the socially optimal hiding probability, lambda^W, 
% in "Surfing Incognito: Welfare Effects of Anonymous Shopping" by Johan N. M.
% Lagerlöf. The latest version of the paper, the code and other material
% are available at www.johanlagerlof.com.

% The following code creates data used in Fig 6a in the paper (so for relatiely high values of delta).


clear all; % Clean up the memory.


Lbest=[0 0 0 0 0 0 0 0 0 0]; % This defines a vector with zeroes. The entries of the 
% vector will later be replaced with the lambda^W expressions that we solve for. 

for nd=1:9 % The start of the "outside" loop. This loop runs through the different delta values.
    
d=nd/10;  % Gives us the delta values (0.1 through 0.9).  

for nc=1:10 % The start of the "inside" loop. This loop runs through ten 
    % different values of c (the parameter in the quadratic cost function).
    
c=2*nc/5; % This defines the c values (from 2/5 through 20/5).

L = 0.0001:0.0001:0.9999; % This creates a vector with 9999 values of 
% lambda. The code will identify the lanbda value among these that 
% maxmizes the welfare function.

% Beolw is the welfare function W, as derived in the Supplementary Material.

W = (2 - L.*(3+L).*d + L.*(1+L).*d.^2) .* (6 + (2+L).*(1-5.*L).*d + L.*(1+L).*(1+2.*L).*d.^2) .* 0.5 .*(4 + (1-6.*L -3.*L.^2).*d + L.*(1+L).^2.*d.^2).^-2    +   d.* L.*(1-L).*(2 + (1-2.*L-L.^2).*d) .* (8 + 2.*(1-5.*L-4.*L.^2).*d - (1-3.*L-3.*L.^2-3.*L.^3).*d.^2) .* (1/8) .*(4 + (1-6.*L -3.*L.^2).*d + L.*(1+L).^2.*d.^2).^-2   +    d.*(2.*(3-L) + (1 - 9.*L - L.^2 + L.^3).*d  + 2.*L.*(1+L).*d.^2) .* (2.*(5+L) + (3 - 15.*L - 11.*L.^2 - L.^3).*d + 2.*L.*(L+1).*(2.*L+1).*d.^2) .* (1/8) .*(4 + (1-6.*L -3.*L.^2).*d + L.*(1+L).^2.*d.^2).^-2       -c.*L.^2;

W;

[Wb(nc), Lb(nc)]=max(W); % Identifies, for each value of the cost parameter,
% which entry in W that gives it its highest value, as well as the obtained
% value at that entry.

Lb(nc) = Lb(nc) ./ 9999; % Translates the entry number to the associated 
% value of lambda.

end

Lbest=[Lbest; Lb]; 

end


Lbest = Lbest([2:end],:) % Erases the first row of zeros in the Lbest matrix.




%%%%%% The following code creates data used in Fig 6b in the paper (so for
%%%%%% relatiely low values of delta).

clear all;

Lbest=[0 0 0 0 0 0 0 0 0 0];

for nd=1:11
    
d=(nd+4)/100; % Gives us the delta values (0.05 through 0.15). 

for nc=1:10
    
c=2*nc/5;

L = 0.00001:0.00001:0.99999;

W = (2 - L.*(3+L).*d + L.*(1+L).*d.^2) .* (6 + (2+L).*(1-5.*L).*d + L.*(1+L).*(1+2.*L).*d.^2) .* 0.5 .*(4 + (1-6.*L -3.*L.^2).*d + L.*(1+L).^2.*d.^2).^-2    +   d.* L.*(1-L).*(2 + (1-2.*L-L.^2).*d) .* (8 + 2.*(1-5.*L-4.*L.^2).*d - (1-3.*L-3.*L.^2-3.*L.^3).*d.^2) .* (1/8) .*(4 + (1-6.*L -3.*L.^2).*d + L.*(1+L).^2.*d.^2).^-2   +    d.*(2.*(3-L) + (1 - 9.*L - L.^2 + L.^3).*d  + 2.*L.*(1+L).*d.^2) .* (2.*(5+L) + (3 - 15.*L - 11.*L.^2 - L.^3).*d + 2.*L.*(L+1).*(2.*L+1).*d.^2) .* (1/8) .*(4 + (1-6.*L -3.*L.^2).*d + L.*(1+L).^2.*d.^2).^-2       -c.*L.^2;

W;

[Wb(nc), Lb(nc)]=max(W);

Lb(nc) = Lb(nc) ./ 99999;

end

hold on
T = [2/5 4/5 6/5 8/5 10/5 12/5 14/5 16/5 18/5 20/5];
plot(T,Lb)

Lbest=[Lbest; Lb];

end
hold off

Lbest = Lbest([2:end],:) ;

Lbest