%PLOTGINI  compare income distributions and Gini indices for model labor income
%          and US compensation.  For the model, compute both before- and after-
%          tax income.
%

%          Ellen McGrattan, 7-24-15
%          Revised, ERM, 7-24-15


%
% Data from CPS 2005 for 12 AGI brackets
%

uscomp  = [ 52.7, 167.7, 359.3, 483.3, 602.1, 596.9, 1098.5, 939.0, 1116.7, ...
           685.3, 748.7, 452.9]';
uspop   = [ 44.0,  18.5,  25.2,  23.5,  22.4,  18.4,   27.3,  18.3,   17.1, ...
             8.2,   5.8,   1.8]';
xusplot = [0;cumsum(uspop)]/sum(uspop)*100;
yusplot = [0;cumsum(uscomp)]/sum(uscomp)*100;
usgini  = gini(xusplot,yusplot);

%
% Draw results for several years (can be changed)
%
ny      = 4;
yr      = [1,10,50,240];

for k=1:ny;
  t                  = yr(k);
  labinput           = zeros(101,ne);
  popfrac            = zeros(101,ne);
  taxsched           = zeros(3,12);
  strt               = int2str(t);
  eval(['labinput(:) = l(',strt,',:);'])
  eval(['popfrac(:)  = popt(',strt,',:)/sum(popt(',strt,',:));'])
  eval(['taxsched(:) = wtaxt(',strt,',:);'])
  eval(['labinc      = w(',strt,')*labinput(1:45,:).*epsmat(1:45,:);'])
  taxsched           = taxsched';
  popfrac            = popfrac(1:45,:);
  for j=1:ne
    for i=1:45;
      tax(i,j)       = comptax(labinc(i,j),taxsched);
    end;
  end;
  mcompbt            = labinc(:).*popfrac(:); 
  mpopbt             = popfrac(:);
  [ym,ii]            = sort(mcompbt./mpopbt);
  mcompbt            = mcompbt(ii);
  mpopbt             = mpopbt(ii);

  mcompat            = (labinc(:)-tax(:)).*popfrac(:); 
  mpopat             = popfrac(:);
  [ymat,ii]          = sort(mcompat./mpopat);
  mcompat            = mcompat(ii);
  mpopat             = mpopat(ii);

  xmplot1(:,k)       = [0;cumsum(mpopbt)]/sum(mpopbt)*100;
  ymplot1(:,k)       = [0;cumsum(mcompbt)]/sum(mcompbt)*100;
  xmplot2(:,k)       = [0;cumsum(mpopat)]/sum(mpopat)*100;
  ymplot2(:,k)       = [0;cumsum(mcompat)]/sum(mcompat)*100;
  mgini1(k,1)        = gini(xmplot1(:,k),ymplot1(:,k));
  mgini2(k,1)        = gini(xmplot2(:,k),ymplot2(:,k));
   
end;

figure(1)
plot(xmplot1(:,1),ymplot1(:,1),'Linewidth',1.5);
for k=1:ny;
  line(xmplot1(:,k),ymplot1(:,k),'Linewidth',1.5);
end;
line([0:100],[0:100],'Linestyle',':')
%line(xusplot,yusplot,'Linewidth',2,'Color','k')


figure(2)
plot(xmplot2(:,1),ymplot2(:,1),'Linewidth',1.5);
for k=2:ny;
  line(xmplot2(:,k),ymplot2(:,k),'Linewidth',1.5);
end;
line([0:100],[0:100],'Linestyle',':')
