clear
close all
clc;

%% figure formatting

set(0,'DefaultTextInterpreter','latex')
set(0,'DefaultAxesTickLabelInterpreter','latex');
set(0,'DefaultLegendInterpreter', 'latex')

set(0,'DefaultTextFontSize', 14)
set(0,'DefaultAxesFontSize', 14)
set(0,'DefaultLineLineWidth',2)

colors = get(gca,'ColorOrder');

c1 = colors(1,:);
c2 = colors(2,:);
c3 = colors(3,:);
close all

country = 1; % 1: US, 2: France
if country == 1 % US
    data = xlsread('TopShareUS.csv','A6:G106');
    TopShares = fliplr(data(:,2:end))/100;
elseif country == 2 % France
    data = xlsread('TopShareFrance.csv','A2:G103');
    TopShares = data(:,2:end);
end

year = data(:,1);
T = length(year);

percent = [0.01 0.1 0.5 1 5 10];

%% estimate Pareto exponent for each year using CUMDE
Kvec = [3 5];
Alpha = zeros(T,numel(Kvec));
sigma_Alpha = zeros(T,numel(Kvec));
Gval = zeros(T,numel(Kvec));
for k=1:numel(Kvec)
    K = Kvec(k);
    critVal = chi2inv(0.95,K-1);
    for t=1:T
        p = percent(1:K+1)/100;
        S = TopShares(t,1:K+1);
        [alpha,G,xi,sigma_alpha] = ParetoCUMDE(p,S);
        Alpha(t,k) = alpha;
        sigma_Alpha(t,k) = sigma_alpha;
        Gval(t,k) = G;
    end
end

%% estimate Pareto exponent using simple estimator
p = percent(2);
q = percent(4);
Sp = TopShares(:,2);
Sq = TopShares(:,4);
AlphaSimple = 1./(1-(log(Sp./Sq))/(log(p/q)));

%% optimal tax rate example
e = 0.3; % income elasticity
tau_CUMDE = 100/(1+Alpha(end,1)*e)
tau_Simple = 100/(1+AlphaSimple(end)*e)

%{
%% save output
if country == 1
    save ParetoUS
elseif country == 2
    save ParetoFrance
end
%}

%% plot figures

figure(1)
plot(year,100*TopShares(:,Kvec(1)),'-','Color',c1); hold on
plot(year,100*TopShares(:,Kvec(2)),':','Color',c2); hold off
xlabel('Year')
ylabel('Top income shares (\%)')
xlim([min(year) max(year)])
legend('Top 1\%','Top 10\%','Location','Best')

%save figure in pdf format
fig = gcf;
fig.PaperPositionMode = 'auto';
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
if country == 1 % US
    print(fig,'Figure_1a','-dpdf')
elseif country == 2
    print(fig,'Figure_2a','-dpdf')
end

figure(2)
plot(year,Alpha(:,1),'-','Color',c1); hold on
plot(year,Alpha(:,2),':','Color',c2);
plot(year,AlphaSimple,'-.','Color',c3);
xlabel('Year')
ylabel('Pareto exponent')
xlim([min(year) max(year)])
legend('Top 1\%','Top 10\%','Simple 0.1/1','Location','NW')

%save figure in pdf format
fig = gcf;
fig.PaperPositionMode = 'auto';
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
if country == 1 % US
    print(fig,'Figure_1b','-dpdf')
elseif country == 2
    print(fig,'Figure_2b','-dpdf')
end

% conservative confidence interval
n = 1e6; % conservative population

figure(3)
plot(year,2*1.96*sigma_Alpha(:,1)/sqrt(n),'-','Color',c1); hold on
plot(year,2*1.96*sigma_Alpha(:,2)/sqrt(n),':','Color',c2);
xlabel('Year')
ylabel('Length of 95\% CI')
xlim([min(year) max(year)])
legend('Top 1\%','Top 10\%','Location','Best')

%save figure in pdf format
fig = gcf;
fig.PaperPositionMode = 'auto';
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
if country == 1 % US
    print(fig,'Figure_3a','-dpdf')
elseif country == 2
    print(fig,'Figure_3b','-dpdf')
end
