%% Program EPIDEMIC_FIGURES.M

% This program generates the five-panel Figure 1 in Goodkin-Gold et al. 
% (2022) Review of Economics and Statistics. 

clc; 

% Set up graph tiling
tiledlayout(5,1)
t.Padding = 'compact';
t.TileSpacing = 'compact';

% Set parameters
theta = .7;
c = 0.3;
Ihat = .1;
Shat = .8;
R0 = (0.01:.01:8)';

% Derived scalars
sizeR0 = size(R0,1);
tc = c / theta;
ltc = log(1-tc);
altc = abs(ltc);

% Thresholds for cases
% Program maximizes monopoly profit so border between interior and 
% boundary solution set endogenously
thresh1 = altc/(Ihat + tc*Shat);
thresh2 = altc/(Ihat + (1-theta)*tc*Shat);

% Asymptotic susceptible population
Slim0 = (1./R0) .* abs(lambertw(-(Shat-theta*0)*R0.*exp(-R0*(Ihat+Shat-theta*0))));
SlimS0 = (1./R0) .* abs(lambertw(-(Shat-theta*Shat)*R0.*exp(-R0*(Ihat+Shat-theta*Shat))));

%%%%%%%%%%
% QUANTITY
%%%%%%%%%%

% Competition
Qc2 = (1 + (Ihat + ltc./R0)/(tc*Shat))/theta;
Qc = 0 + (((R0 > thresh1).*(R0 <= thresh2)) .* Qc2) + (R0 > thresh2);
nexttile
plot(R0,Qc,'--k','LineWidth',2)
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
ylim([0 1.2])
hold on

% Monopoly
Qm2 = zeros(sizeR0,1);
options = optimoptions('fmincon','Display','none');
for i=1:sizeR0
    Slimx = @(x)(1/R0(i))*abs(lambertw(-(Shat-theta*x)*R0(i)*exp(-R0(i)*(Ihat+Shat-theta*x))));
    MPBx = @(x)1-Slimx(x)/(Shat-theta*x);
    Profitnegx = @(x)-(MPBx(x)-tc)*x;
    Qm2(i) = fmincon(Profitnegx,.2,[],[],[],[],0,Shat,[],options);
end
Qm = 0 + (R0 > thresh1).*Qm2;
SlimQm = (1./R0) .* abs(lambertw(-(Shat-theta*Qm) .* R0 .* exp(-R0 .* (Ihat+Shat-theta*Qm))));
tQm = Qm / Shat;
plot(R0,tQm,'-k','LineWidth',2)
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold off

%%%%%%%%%%%
% Recovered
%%%%%%%%%%%

% Competition
Recovered1c = 1 - Slim0;
Recovered2c = 1 - Shat - Ihat + altc ./ R0;
Recovered3c = 1 - SlimS0 - theta * Shat;
Recoveredc = (Recovered1c.*(R0 <= thresh1)) + (Recovered2c.*(R0>thresh1).*(R0<=thresh2)) + (Recovered3c.*(R0>thresh2));
nexttile
plot(R0,Recoveredc,'--k','LineWidth',2)
ylim([0 inf])
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold on

% Monopoly
Recovered1m = 1 - Slim0;
Recovered2m = 1 - SlimQm - theta * Qm;
Recoveredm = (Recovered1m.*(R0 <= thresh1)) + (Recovered2m.*(R0>thresh1)); 
plot(R0,Recoveredm,'-k','LineWidth',2)
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold off

%%%%%
% MPB 
%%%%%

% Competition
MPB1c = 1 - Slim0/Shat;
MPB2c = tc;
MPB3c = 1 - SlimS0/((1-theta)*Shat);
MPBc = (MPB1c.*(R0 <= thresh1)) + (MPB2c.*(R0>thresh1).*(R0<=thresh2)) + (MPB3c.*(R0>thresh2));
nexttile
plot(R0,MPBc,'--k','LineWidth',2);
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
ylim([0 1.2])
hold on

% Monopoly
MPB1m = 1 - Slim0/Shat;
MPB2m = 1 - SlimQm ./ (Shat - theta * Qm);
MPBm = (MPB1m.*(R0 <= thresh1)) + (MPB2m.*(R0>thresh1)); 
plot(R0,MPBm,'-k','LineWidth',2)
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold off

%%%%%
% MEX
%%%%%

% Cmpetition
MEX1c = (MPB1c.*R0.*Slim0) ./ (1 - R0 .* Slim0);
MEX2c = (tc*(1-tc)*(altc - Ihat .* R0)) ./ (tc + (1-tc)*(ltc + Ihat .* R0));
MEX3c = (MPB3c.*R0.*SlimS0) ./ (1 - R0 .* SlimS0);
MEX = (MEX1c.*(R0 <= thresh1)) + (MEX2c.*(R0>thresh1).*(R0<=thresh2)) + (MEX3c.*(R0>thresh2));
nexttile
plot(R0,MEX,'--k','LineWidth',2);
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold on

% Monopoly
MEX1m = (MPB1c.*R0.*Slim0) ./ (1 - R0 .* Slim0);
MEX2m = (MPB2m.*R0.*SlimQm) ./ (1 - R0 .* SlimQm);
MEXm = (MEX1m.*(R0 <= thresh1)) + (MEX2m.*(R0>thresh1)); 
plot(R0,MEXm,'-k','LineWidth',2)
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold off

%%%%%%%%%
% Welfare 
%%%%%%%%%

% Competition
W1c = Slim0 / (theta*Shat);
W2c = (1-tc)/theta;
W3c = 1 - tc + SlimS0 / (theta*Shat);
Wc = (W1c.*(R0 <= thresh1)) + (W2c.*(R0>thresh1).*(R0<=thresh2)) + (W3c.*(R0>thresh2));
nexttile
plot(R0,Wc,'--k','LineWidth',2);
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold on

% Monopoly
W1m = Slim0 / (theta*Shat);
W2m = (SlimQm / (theta * Shat)) + (1-tc)*tQm;
Wm = (W1m.*(R0 <= thresh1)) + (W2m.*(R0>thresh1)); 
plot(R0,Wm,'-k','LineWidth',2)
set(gca,'xtick',[],'xticklabel',[],'ytick',[],'yticklabel',[],'box','off')
hold off

