% -------------------------------------------------------------------------
%
% Project: Uneven Growth
%
% Purpose: Create Figure 9
%
% -------------------------------------------------------------------------

clear variables
close all

cd('C:/Users/Pascual Restrepo/Dropbox/Research/wealth_inequality/AutomationInequality/replication_mrr_january_2021/additional_evidence/psz_irs_data/')
addpath 'C:/Users/Pascual Restrepo/Documents/MATLAB/scripts'
addpath 'output'

% -------------------------------------------------------------------------
%% - Preamble -
% -------------------------------------------------------------------------

% Color choices %
lightBlue = [40, 100, 150] / 255 ;
alpha=0.6;
scorp_color=([0.4660, 0.6740, 0.1880])*alpha+(1-alpha)*[1,1,1];
labor_color=lightBlue*alpha+(1-alpha)*[1,1,1];
cap_color=[0.9500, 0.4250, 0.40]*alpha+(1-alpha)*[1,1,1];

% Set latex as default interpreter
set(groot, 'defaulttextinterpreter','latex'); 
set(groot, 'defaultAxesTickLabelInterpreter','latex'); 
set(groot, 'defaultLegendInterpreter','latex');

% grid for extended x axis (need to fix this so that this represents
% average over population)
int_grid = [5 10:10:90 95 97.5 99 99.5 99.9 99.99 100];

% -------------------------------------------------------------------------
%% - Figure 10 and versions thereof -
% -------------------------------------------------------------------------

% - Read data -------------------------------------------------------------

irs_1980_2012 = csvread("irs_clean_1980_2012.csv",2);
irs_1980_2007 = csvread("irs_clean_1980_2007.csv",2);
psz_1980_2012 = csvread("psz_clean_1980_2012.csv",2);
psz_1980_2007 = csvread("psz_clean_1980_2007.csv",2);

% - Figure 10 -------------------------------------------------------------

% - Panel (a): IRS figure -------------------------------------------------
perc = irs_1980_2012(:,1);
growth = irs_1980_2012(:,2);
growth_lab = irs_1980_2012(:,3);
growth_scor = irs_1980_2012(:,4);
growth_cap = irs_1980_2012(:,5);

% - compute income growth for percentiles on the grid - %
growth=interp1(perc, growth , int_grid ,'spline');
growth_lab=interp1(perc, growth_lab , int_grid ,'spline');
growth_scor=interp1(perc, growth_scor , int_grid ,'spline');
growth_cap=interp1(perc, growth_cap , int_grid ,'spline');

% - approximate region by endpoint averages - %
growth=(growth(2:end)+growth(1:end-1))/2;
growth_lab=(growth_lab(2:end)+growth_lab(1:end-1))/2;
growth_cap=(growth_cap(2:end)+growth_cap(1:end-1))/2;
growth_scor=(growth_scor(2:end)+growth_scor(1:end-1))/2;

% - compute areas to be plotted - %
growth_scor2 = max(growth_lab,0)+growth_scor;
ind_dif = (growth_lab.*growth_scor)<0;
growth_scor2(ind_dif) = growth_scor(ind_dif);

% - define regions to be plotted - %
growth_cap2 = growth;
ind_dif = (growth.*growth_cap)<0;
growth_cap2(ind_dif) = growth_cap(ind_dif);

figure_irs = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

xaxis = 1:16;

axes('Position', [0.06 0.22 .9 .7]);
    hold on
    box on
   h1=area(xaxis,growth_cap2,'Facecolor',cap_color, 'Edgecolor', cap_color);
   h2=area(xaxis,growth_scor2,'Facecolor',scorp_color,'Edgecolor',scorp_color);
   h3=area(xaxis,growth_lab,'Facecolor',labor_color,'Edgecolor',labor_color);
   h4=plot(xaxis,growth,'-.', 'linewidth',1,'color',[0,0,0]);

   
    ylim([-1 6])
    xlim([1 16])

    set(gca,'XTick',1:length(xaxis),'XTickLabel',{'P5-P10','P10-P20','P20-P30','P30-P40','P40-P50','P50-P60','P60-P70','P70-P80' ...
        ,'P80-P90','P90-P95','P95-P97.5','P97.5-P99','P99-P99.5','P99.5-P99.9','P99.9-P99.99','P99.99-P100'}, ...
        'XTickLabelRotation',45)

legend([h4 h3 h2 h1],{'total income','part due to labor income','part due to 75\% of s-corporations income','part due to capital income'},'Location','NorthWest')
legend('boxoff')
xlabel('income percentile')
ylabel('annual growth 1980-2012 (in \%)')

set(gcf,'color','w');

export_fig 'output/Figure_paper_consistent_irs_share_1980_2012_onepanel'.pdf 

% - Panel (b): PSZ figure -------------------------------------------------
perc = psz_1980_2012(:,1);
growth = psz_1980_2012(:,2);
growth_lab = psz_1980_2012(:,5);
growth_scor = psz_1980_2012(:,4);
growth_cap = psz_1980_2012(:,3);

% - compute income growth for percentiles on the grid - %
growth=interp1(perc, growth , int_grid ,'spline');
growth_lab=interp1(perc, growth_lab , int_grid ,'spline');
growth_scor=interp1(perc, growth_scor , int_grid ,'spline');
growth_cap=interp1(perc, growth_cap , int_grid ,'spline');

% - approximate region by endpoint averages - %
growth=(growth(2:end)+growth(1:end-1))/2;
growth_lab=(growth_lab(2:end)+growth_lab(1:end-1))/2;
growth_cap=(growth_cap(2:end)+growth_cap(1:end-1))/2;
growth_scor=(growth_scor(2:end)+growth_scor(1:end-1))/2;

% - compute areas to be plotted - %
growth_scor2 = max(growth_lab,0)+growth_scor;
ind_dif = (growth_lab.*growth_scor)<0;
growth_scor2(ind_dif) = growth_scor(ind_dif);

growth_cap2 = growth;
ind_dif = (growth.*growth_cap)<0;
growth_cap2(ind_dif) = growth_cap(ind_dif);

figure_psz = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

xaxis = 1:16;

axes('Position', [0.06 0.22 .9 .7]);
    hold on
    box on
   h1=area(xaxis,growth_cap2,'Facecolor',cap_color, 'Edgecolor', cap_color);
   h2=area(xaxis,growth_scor2,'Facecolor',scorp_color,'Edgecolor',scorp_color);
   h3=area(xaxis,growth_lab,'Facecolor',labor_color,'Edgecolor',labor_color);
   h4=plot(xaxis,growth,'-.', 'linewidth',1,'color',[0,0,0]);

   
    ylim([-1 6])
    xlim([1 16])

    set(gca,'XTick',1:length(xaxis),'XTickLabel',{'P5-P10','P10-P20','P20-P30','P30-P40','P40-P50','P50-P60','P60-P70','P70-P80' ...
        ,'P80-P90','P90-P95','P95-P97.5','P97.5-P99','P99-P99.5','P99.5-P99.9','P99.9-P99.99','P99.99-P100'}, ...
        'XTickLabelRotation',45)

legend([h4 h3 h2 h1],{'total income','part due to labor income','part due to 75\% of s-corporations income','part due to capital income'},'Location','NorthWest')
legend('boxoff')
xlabel('income percentile')
ylabel('annual growth 1980-2012 (in \%)')

set(gcf,'color','w');

export_fig 'output/Figure_paper_consistent_psz_share_1980_2012_onepanel'.pdf 

% - Figure 10 alternative: 1980 - 2007 ------------------------------------

% - Panel (a): IRS figure -------------------------------------------------
perc = irs_1980_2007(:,1);
growth = irs_1980_2007(:,2);
growth_lab = irs_1980_2007(:,3);
growth_scor = irs_1980_2007(:,4);
growth_cap = irs_1980_2007(:,5);

% - compute income growth for percentiles on the grid - %
growth=interp1(perc, growth , int_grid ,'spline');
growth_lab=interp1(perc, growth_lab , int_grid ,'spline');
growth_scor=interp1(perc, growth_scor , int_grid ,'spline');
growth_cap=interp1(perc, growth_cap , int_grid ,'spline');

% - approximate region by endpoint averages - %
growth=(growth(2:end)+growth(1:end-1))/2;
growth_lab=(growth_lab(2:end)+growth_lab(1:end-1))/2;
growth_cap=(growth_cap(2:end)+growth_cap(1:end-1))/2;
growth_scor=(growth_scor(2:end)+growth_scor(1:end-1))/2;

% - compute areas to be plotted - %
growth_scor2 = max(growth_lab,0)+growth_scor;
ind_dif = (growth_lab.*growth_scor)<0;
growth_scor2(ind_dif) = growth_scor(ind_dif);

growth_cap2 = growth;
ind_dif = (growth.*growth_cap)<0;
growth_cap2(ind_dif) = growth_cap(ind_dif);

figure_irs07 = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

xaxis = 1:16;

axes('Position', [0.06 0.22 .9 .7]);
    hold on
    box on
   h1=area(xaxis,growth_cap2,'Facecolor',cap_color, 'Edgecolor', cap_color);
   h2=area(xaxis,growth_scor2,'Facecolor',scorp_color,'Edgecolor',scorp_color);
   h3=area(xaxis,growth_lab,'Facecolor',labor_color,'Edgecolor',labor_color);
   h4=plot(xaxis,growth,'-.', 'linewidth',1,'color',[0,0,0]);

   
    ylim([-1 6])
    xlim([1 16])

    set(gca,'XTick',1:length(xaxis),'XTickLabel',{'P5-P10','P10-P20','P20-P30','P30-P40','P40-P50','P50-P60','P60-P70','P70-P80' ...
        ,'P80-P90','P90-P95','P95-P97.5','P97.5-P99','P99-P99.5','P99.5-P99.9','P99.9-P99.99','P99.99-P100'}, ...
        'XTickLabelRotation',45)

legend([h4 h3 h2 h1],{'total income','part due to labor income','part due to 75\% of s-corporations income','part due to capital income'},'Location','NorthWest')
legend('boxoff')
xlabel('income percentile')
ylabel('annual growth 1980-2012 (in \%)')

set(gcf,'color','w');

export_fig 'output/Figure_paper_consistent_irs_share_1980_2007_onepanel'.pdf 


% - Panel (b): PSZ figure -------------------------------------------------
perc = psz_1980_2007(:,1);
growth = psz_1980_2007(:,2);
growth_lab = psz_1980_2007(:,5);
growth_scor = psz_1980_2007(:,4);
growth_cap = psz_1980_2007(:,3);

% - compute income growth for percentiles on the grid - %
growth=interp1(perc, growth , int_grid ,'spline');
growth_lab=interp1(perc, growth_lab , int_grid ,'spline');
growth_scor=interp1(perc, growth_scor , int_grid ,'spline');
growth_cap=interp1(perc, growth_cap , int_grid ,'spline');

% - approximate region by endpoint averages - %
growth=(growth(2:end)+growth(1:end-1))/2;
growth_lab=(growth_lab(2:end)+growth_lab(1:end-1))/2;
growth_cap=(growth_cap(2:end)+growth_cap(1:end-1))/2;
growth_scor=(growth_scor(2:end)+growth_scor(1:end-1))/2;

% - compute areas to be plotted - %
growth_scor2 = max(growth_lab,0)+growth_scor;
ind_dif = (growth_lab.*growth_scor)<0;
growth_scor2(ind_dif) = growth_scor(ind_dif);

growth_cap2 = growth;
ind_dif = (growth.*growth_cap)<0;
growth_cap2(ind_dif) = growth_cap(ind_dif);

figure_psz07 = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

xaxis = 1:16;

axes('Position', [0.06 0.22 .9 .7]);
    hold on
    box on
   h1=area(xaxis,growth_cap2,'Facecolor',cap_color, 'Edgecolor', cap_color);
   h2=area(xaxis,growth_scor2,'Facecolor',scorp_color,'Edgecolor',scorp_color);
   h3=area(xaxis,growth_lab,'Facecolor',labor_color,'Edgecolor',labor_color);
   h4=plot(xaxis,growth,'-.', 'linewidth',1,'color',[0,0,0]);

   
    ylim([-1 6])
    xlim([1 16])

    set(gca,'XTick',1:length(xaxis),'XTickLabel',{'P5-P10','P10-P20','P20-P30','P30-P40','P40-P50','P50-P60','P60-P70','P70-P80' ...
        ,'P80-P90','P90-P95','P95-P97.5','P97.5-P99','P99-P99.5','P99.5-P99.9','P99.9-P99.99','P99.99-P100'}, ...
        'XTickLabelRotation',45)

legend([h4 h3 h2 h1],{'total income','part due to labor income','part due to 75\% of s-corporations income','part due to capital income'},'Location','NorthWest')
legend('boxoff')
xlabel('income percentile')
ylabel('annual growth 1980-2012 (in \%)')

set(gcf,'color','w');

export_fig 'output/Figure_paper_consistent_psz_share_1980_2007_onepanel'.pdf 


% - Figure 10 alternative: 1980 - 2012 with alternative s-corp -----------

perc = psz_1980_2007(:,1);
growth = psz_1980_2007(:,2);
growth_lab = psz_1980_2007(:,5);
growth_cap = psz_1980_2007(:,6);
growth_scor = psz_1980_2007(:,7);

% - compute income growth for percentiles on the grid - %
growth=interp1(perc, growth , int_grid ,'spline');
growth_lab=interp1(perc, growth_lab , int_grid ,'spline');
growth_scor=interp1(perc, growth_scor , int_grid ,'spline');
growth_cap=interp1(perc, growth_cap , int_grid ,'spline');

% - approximate region by endpoint averages - %
growth=(growth(2:end)+growth(1:end-1))/2;
growth_lab=(growth_lab(2:end)+growth_lab(1:end-1))/2;
growth_cap=(growth_cap(2:end)+growth_cap(1:end-1))/2;
growth_scor=(growth_scor(2:end)+growth_scor(1:end-1))/2;

% - compute areas to be plotted - %
growth_scor2 = max(growth_lab,0)+growth_scor;
ind_dif = (growth_lab.*growth_scor)<0;
growth_scor2(ind_dif) = growth_scor(ind_dif);

growth_cap2 = growth;
ind_dif = (growth.*growth_cap)<0;
growth_cap2(ind_dif) = growth_cap(ind_dif);

figure_psz_alt = figure('units','centimeters','outerposition',[10 10 16*1.2 12*1.2]);

xaxis = 1:16;

axes('Position', [0.06 0.22 .9 .7]);
    hold on
    box on
   h1=area(xaxis,growth_cap2,'Facecolor',cap_color, 'Edgecolor', cap_color);
   h2=area(xaxis,growth_scor2,'Facecolor',scorp_color,'Edgecolor',scorp_color);
   h3=area(xaxis,growth_lab,'Facecolor',labor_color,'Edgecolor',labor_color);
   h4=plot(xaxis,growth,'-.', 'linewidth',1,'color',[0,0,0]);

   
    ylim([-1 6])
    xlim([1 16])

    set(gca,'XTick',1:length(xaxis),'XTickLabel',{'P5-P10','P10-P20','P20-P30','P30-P40','P40-P50','P50-P60','P60-P70','P70-P80' ...
        ,'P80-P90','P90-P95','P95-P97.5','P97.5-P99','P99-P99.5','P99.5-P99.9','P99.9-P99.99','P99.99-P100'}, ...
        'XTickLabelRotation',45)

legend([h4 h3 h2 h1],{'total income','part due to labor income','part due to 75\% of s-corporations income','part due to capital income'},'Location','NorthWest')
legend('boxoff')
xlabel('income percentile')
ylabel('annual growth 1980-2012 (in \%)')

set(gcf,'color','w');

export_fig 'output/Figure_paper_consistent_psz_share_1980_2012_alternative_onepanel'.pdf 


%% Note: alternative code for integration step:
% - integrate over grid - %
%cum_growth = cumtrapz(perc,growth);
%cum_growth_int=interp1(perc, cum_growth , int_grid ,'spline');
%growth=(cum_growth_int(2:end)-cum_growth_int(1:end-1))./(int_grid(2:end)-int_grid(1:end-1));
%cum_growth_lab = cumtrapz(perc,growth_lab);
%cum_growth_lab_int=interp1(perc, cum_growth_lab , int_grid ,'spline');
%growth_lab=(cum_growth_lab_int(2:end)-cum_growth_lab_int(1:end-1))./(int_grid(2:end)-int_grid(1:end-1));
%cum_growth_cap = cumtrapz(perc,growth_cap);
%cum_growth_cap_int=interp1(perc, cum_growth_cap , int_grid ,'spline');
%growth_cap=(cum_growth_cap_int(2:end)-cum_growth_cap_int(1:end-1))./(int_grid(2:end)-int_grid(1:end-1));
%cum_growth_scor = cumtrapz(perc,growth_scor);
%cum_growth_scor_int=interp1(perc, cum_growth_scor , int_grid ,'spline');
%growth_scor=(cum_growth_scor_int(2:end)-cum_growth_scor_int(1:end-1))./(int_grid(2:end)-int_grid(1:end-1));
