%% Program EPIDEMIC_COVID_CALIBRATIONS.M

% This program carries out the Covid-19 calibrations reported in Section 8
% of Goodkin-Gold et al. (2022) Review of Economics and Statistics. 


%% Parameters
theta = .8;

% Toggle for ancestral Covid parameters as of Oct. 2020
R0 = 2.8;
Ihat = .0019;
Rhat = .062;

% Toggle for alternative parameters for Delta
%R0 = 5.1;

% Display parameter values
fprintf('Parameters');
fprintf('----------');
fprintf('efficacy (theta) = %G\n',theta);
fprintf('basic reproductive number (R0) = %G\n',R0);
fprintf('initial infected proportion (Ihat) = %G\n',Ihat);
fprintf('initial recovered proportion (Rhat) = %G\n',Rhat);
Shat = 1 - Ihat - Rhat;
fprintf('initial susceptible proportion (Shat) = %G\n',Shat);

%% Define profit function and optimization routine for maximizing it

Slim = @(x)(1/R0)*abs(lambertw(-(Shat-theta*x)*R0*exp(-R0*(Ihat+Shat-theta*x))));
Phi = @(x)1-Slim(x)/(Shat-theta*x);
Profitneg = @(x)-x*Phi(x);
Qm = fmincon(Profitneg,.2,[],[],[],[],0,Shat);
Profit = @(x)-theta*Profitneg(x);

%% Display monopoly outcomes

fprintf('\n');
fprintf('Monopoly outcomes\n');
fprintf('-----------------\n');
fprintf('Monopoly quantity (Qm) = %G\n',Qm);
fprintf('Monopoly quantity as fraction of susceptibles = %G\n',Qm/Shat);
fprintf('Monopoly profit = %G\n',Profit(Qm));

Pm = theta * Phi(Qm);
fprintf('Monopoly price (Pm) = %G\n',Pm);
Welfarem = theta*Qm+Slim(Qm);
fprintf('Monopoly welfare (Wm) = %G\n',Welfarem);

GSm = Pm * ((theta/(1-theta)) * (R0 * Slim(Shat)/(1-R0*Slim(Shat))) - 1);
GSm = max(0,GSm); 
fprintf('Optimal subsidy under monopoly (GSm) = %G\n',GSm);

%% Increasing social returns results

% Determine threshold below which increasing social returns to vaccination
syms Q
eqn1 = R0*(Slim(Q) + Shat - theta * Q)/2;
IRSthresh = vpasolve(eqn1==1,Q,.6);
fprintf('\n');
fprintf('Increasing returns results\n');
fprintf('--------------------------\n');
fprintf('Threshold for increasing social returns to vaccination = %G\n',IRSthresh);

% Determine threshold stockpile size below which concentrate. 
eqn2 = Slim(Q) + theta * Q + Slim(0) - 2 *(Slim(Q/2) + theta * Q/2);
conc = vpasolve(eqn2==0,Q,.8);
fprintf('Threshold below which concentrate supplies in one of two regions = %G\n',conc);

%% Drug versus vaccine results 

Delta_profit = theta * (Ihat + Shat - Slim(0)) - Profit(Qm);
fprintf('\n');
fprintf('Drug vs. vaccine results\n');
fprintf('------------------------\n');
fprintf('Delta profit = %G\n',Delta_profit);

% Drug versus vaccine welfare
Delta_welfare = (1 - theta)*Slim(0) + theta * (Ihat + Shat) - Welfarem;
fprintf('Delta welfare = %G\n',Delta_welfare);