%TAXFUNC    Tax functions for wage and business income and average 
%           weighted marginal tax on dividends for 2007 used in:
%
%             Sweat Equity in US Private Business
%             by Anmol Bhandari and Ellen McGrattan
%
%           Notes: Run this code on a machine with excel.

%           Ellen McGrattan, 4-3-19
%           Revised, ERM, 4-3-19


%=======================================================================
% Read in raw data:
%
%   IRS Returns of Individuals 
%     Table 4.  Sources of incomes
%     Table 23. Exemptions by type
%
%   IRS Returns of S-corporations
%     Table 7S. Balance sheet and income statements
%
%   IRS tax schedule by marital status, 2007
%
%   BEA NIPA 
%     Table 7.14 Relation of nonfarm proprietors' income in NIPA and IRS
%
%   Fed FOF
%     Table B.101e. Balance sheet of Households, equity detail
%     Table L.118b. Private pension funds: DB plans
%     Table L.119b. Federal retirement funds: DB plans
%     Table L.120b. S&L retirement funds: DB plans
%     Table L.122.  Mutual funds
%     Table L.227.  Pension entitlements 
%
%=======================================================================

load ../RawData/IRS/fedrates07.dat

tab14        = xlsread('../RawData/IRS/in14ar/2007in14ar');
tab23        = xlsread('../RawData/IRS/in23ar/2007in23ar');
tab1120s     = xlsread('../RawData/IRS/co07s/2007co07s');
beamisrep    = xlsread('../RawData/BEA/NIPA/Section7all_xls','T71400-A');
btab_101eQ   = xlsread('../RawData/FOF/xlsx/b101e');
ltab_122Q    = xlsread('../RawData/FOF/xlsx/l122');
ltab_227Q    = xlsread('../RawData/FOF/xlsx/l227');
ltab_118bQ   = xlsread('../RawData/FOF/xlsx/l118b');
ltab_119bQ   = xlsread('../RawData/FOF/xlsx/l119b');
ltab_120bQ   = xlsread('../RawData/FOF/xlsx/l120b');

%=======================================================================
%
% Read off relevant rows/columns 
%
%=======================================================================

tab14        = [tab14(2,:);tab14(3,:)+tab14(4,:);tab14(5:21,:)];
tab23        = [tab23(2:8,:);tab23(9,:)+tab23(10,:);tab23(11,:)+ ...
                tab23(12,:);tab23(13,:)+tab23(14,:)+tab23(15,:); ...
                tab23(16:24,:)];

nragi        = tab14(:,1);
agi          = tab14(:,2);
nrsw         = tab14(:,3);
sw           = tab14(:,4);
nrbusi       = tab14(:,17);
busi         = tab14(:,18);
nrbusl       = tab14(:,19);
busl         = tab14(:,20);
nrspi        = tab14(:,91);
spi          = tab14(:,92);
nrspl        = tab14(:,93);
spl          = tab14(:,94);
exchk        = tab14(:,159);
nrtaxi       = tab14(:,161);
taxi         = tab14(:,162);
exem         = tab23(:,2);
kids         = tab23(:,10);
adults       = exem-kids;
sccomp1120s  = tab1120s(41,1);
scneti1120s  = tab1120s(57,1);

%=======================================================================
%
% Use taxable income to look up federal tax rates (done manually here)
%
%=======================================================================

taxichk      = fedrates07(:,4);
[n,m]        = size(fedrates07);
iwarn        = 1;
if (length(taxi)==n);
  if (sum(abs(taxi-taxichk))<.001); iwarn=0; end;
end;
if (iwarn==1); 
  disp('WARNING: Check taxable income ...')
end;

marwgt       = fedrates07(1,5:8);
fedrates     = [NaN*ones(1,4);fedrates07(2:n,5:8)];
mtr_fed      = fedrates*marwgt'/100;
mtr_stl      = 4*ones(n,1);
fica         = fedrates07(:,9);
mtr_nofica   = mtr_fed+mtr_stl;
mtr_fica     = mtr_fed+mtr_stl+fica;
mtr_qdiv     = (5.*(fedrates<=15)+15.*(fedrates>15))*marwgt'/100;

%=======================================================================
%
% Aggregate business returns with income and losses
%
%=======================================================================
bus          = busi-busl;
sp           = spi-spl;
nrbus        = nrbusi+nrbusl;
nrsp         = nrspi+nrspl;


%=======================================================================
%
% Include S-corporate wages (expensed on 1120S) with S-corp income
%
%=======================================================================

sppos        = sum(max(sp,0))-sp(1);
sccomp       = round(sccomp1120s/sppos*max(sp,0));
sccomp(1)    = sccomp1120s;
scshr        = scneti1120s/sp(1);
parshr       = 1-scshr;
scneti       = scshr*sp;
parneti      = parshr*sp;

%=======================================================================
%
% Use information on business income underreporting (Slemrod, GAO)
%   Sole proprietors: 57% of income not reported on Schedule C
%   S-corporations:   18% of income not reported on 1120s
%
%=======================================================================

busrate      = .57;
scrate       = .18;
bustru       = bus/(1-busrate);
sccomptru    = sccomp/(1-scrate);
scnetitru    = scneti/(1-scrate);

%=======================================================================
%
% Infer partners' underreporting from BEA Table 7.14
%   beamisrep  = misreporting by proprietors and partners according to BEA 
%   schcmisrep = misreporting by proprietors according to Slemrod
%   partmisrep = residual
%   
%   If $371 bil missing from Schedule C ==> 144 for partners
%   Add 144 to reported 1040 Schedule E (415) less S-corp 1120s income (287)
%   ==> true income is 272=144+128
%   ==> 53% is underreported (ie., 144/272)
%
%=======================================================================

beamisrep    = beamisrep(3,49)*1000;
busmisrep    = bustru(1)-bus(1);
parmisrep    = beamisrep-busmisrep;
parrate      = parmisrep/(parneti(1)+parmisrep);
partru       = parneti/(1-parrate);

%=======================================================================
%
% Compute effective tax rate for each AGI bin j across income types i
%   = Sum_i Compliant tax rate_{ij} *Reported income_{ij}
%    /Sum_i True income_{ij}
%
%=======================================================================

tauw         = mtr_fica;
taub         = (mtr_fica.*(bus+parneti+sccomp)+mtr_nofica.*scneti)./ ...
               (bustru+partru+sccomptru+scnetitru);


%=======================================================================
%
% Compute average marginal tax rates
%
%=======================================================================
[n,m]        = size(tauw);
btru         = bustru+partru+sccomptru+scnetitru;
amtrw        = sum(tauw(2:n).*sw(2:n))/sum(sw(2:n));
amtrb        = sum(taub(2:n).*btru(2:n))/sum(btru(2:n));

%=======================================================================
%
% Compute true income per return and per adult
%
%=======================================================================

yprw         = sw./nrsw;
yprb         = (bustru+partru+sccomptru+scnetitru)./(nrbus+nrsp);
ret2adult    = (nrsw+nrsp+nrbus)./adults;
yppw         = yprw.*ret2adult;
yppb         = yprb.*ret2adult;

%=======================================================================
%
% Normalize income per adult by per capita GDP (adjusted for intangibles)
%
%=======================================================================

tab1_1_5A    = xlsread('../RawData/BEA/NIPA/Section1all_xls','T10105-A');
[n,m]        = size(tab1_1_5A);
t            = 1929:1928+m;
i            = find(t==2007);
tab1_1_5A    = tab1_1_5A(:,i)/1000;
gdp          = tab1_1_5A(1);
cpu          = xlsread('../RawData/BLS/CPU_Data_construction','BLS-CPS');
cpu          = cpu(3:72,:);
t            = 1948:2017;
i            = find(t==2007);
pop1664      = cpu(i,20)/10^6;

nyppw        = yppw/(gdp*1.12/pop1664);
nyppb        = yppb/(gdp*1.12/pop1664);


%=======================================================================
%
% Compute marginal tax rate on dividends 
%
%=======================================================================

[n,m]        = size(tab14);
odiv         = tab14(:,10)-tab14(:,12);
qdiv         = tab14(:,12);
taudo        = (mtr_fed(2:n)+5).*odiv(2:n); 
taudq        = (mtr_qdiv(2:n)+5).*qdiv(2:n); %using Taxsim state mtr
taud         = taudo+taudq;
taudo        = sum(taudo)/odiv(1);
taudq        = sum(taudq)/qdiv(1);
taud         = sum(taud)/(odiv(1)+qdiv(1));

%=======================================================================
%
% Multiply marginal rate by taxable share of equities held by HHs
%
%=======================================================================


[n,m]        = size(btab_101eQ);
btab_101eA   = btab_101eQ(1:4:n,:)';
[n,m]        = size(btab_101eA);
t            = 1945:1944+m;
i            = find(t==2007);
btab_101eA   = btab_101eA(:,i)/1000;

[n,m]        = size(ltab_118bQ);
ltab_118bA    = ltab_118bQ(1:4:n,:)';
[n,m]        = size(ltab_118bA);
ltab_118bA    = ltab_118bA(:,i)/1000;

[n,m]        = size(ltab_119bQ);
ltab_119bA    = ltab_119bQ(1:4:n,:)';
[n,m]        = size(ltab_119bA);
ltab_119bA    = ltab_119bA(:,i)/1000;

[n,m]        = size(ltab_120bQ);
ltab_120bA    = ltab_120bQ(1:4:n,:)';
[n,m]        = size(ltab_120bA);
ltab_120bA    = ltab_120bA(:,i)/1000;

[n,m]        = size(ltab_122Q);
ltab_122A    = ltab_122Q(1:4:n,:)';
[n,m]        = size(ltab_122A);
ltab_122A    = ltab_122A(:,i)/1000;

[n,m]        = size(ltab_227Q);
ltab_227A    = ltab_227Q(1:4:n,:)';
[n,m]        = size(ltab_227A);
ltab_227A    = ltab_227A(:,i)/1000;

prvdbce      = ltab_118bA(12);
prvdbmf      = ltab_118bA(13);
feddbce      = ltab_119bA(8);
stldbce      = ltab_120bA(13);
stldbmf      = ltab_120bA(14);
irasmf       = ltab_227A(16);
irassd       = ltab_227A(17);
mfce         = ltab_122A(10);
mftot        = ltab_122A(12);

eqiras       = (irasmf+irassd)*mfce/mftot;
eqdefben     = prvdbce+feddbce+stldbce+(prvdbmf+stldbmf)*mfce/mftot;
hhequity     = btab_101eA(7)+eqdefben+eqiras;
untaxed      = sum(btab_101eA(10:13))+eqdefben+eqiras;

hheqtaxed    = (hhequity-untaxed)/hhequity;
taudeff      = taud*hheqtaxed;


%=======================================================================
%
% Print results 
%
%=======================================================================
n            = length(tauw);
lbw          = [-inf;.5*nyppw(2:n-1)+.5*nyppw(3:n)];
ubw          = [.5*nyppw(2:n-1)+.5*nyppw(3:n);inf];
lbb          = [-inf;.5*nyppb(2:n-1)+.5*nyppb(3:n)];
ubb          = [.5*nyppb(2:n-1)+.5*nyppb(3:n);inf];

clc
disp('  Tax Function for Wages T^w(.)   ')
disp('----------------------------------')
disp('   Wage brackets       Tax rate   ')
disp('----------------------------------')
for i=1:n-1;
  disp(sprintf('%7.3f %7.3f       %7.1f',lbw(i),ubw(i),tauw(i+1)))
end;
disp('----------------------------------')
disp(' ')
disp(' ')
disp(' ')
disp(' Tax Function for Business T^b(.) ')
disp('----------------------------------')
disp('   Income brackets     Tax rate   ')
disp('----------------------------------')
for i=1:n-1;
  disp(sprintf('%7.3f %7.3f       %7.1f',lbb(i),ubb(i),taub(i+1)))
end;
disp('----------------------------------')
disp(' ')
disp(' ')
disp(' ')
disp('----------------------------------')
disp('      Avg Marginal Rates          ')
disp('----------------------------------')
disp(sprintf('Wages                 %7.1f',amtrw))
disp(sprintf('Business              %7.1f',amtrb))
disp('----------------------------------')
disp(' ')
disp(' ')
disp(' ')
disp('----------------------------------')
disp('      Tax Rates on Dividends      ')
disp('----------------------------------')
disp('Dividend type:         Tax rate   ')
disp(sprintf('  Ordinary            %7.1f',taudo))
disp(sprintf('  Qualified           %7.1f',taudq))
disp(sprintf('  Effective           %7.1f',taudeff))
disp(' ')
disp(sprintf('%% taxable             %7.1f',hheqtaxed*100))
disp('----------------------------------')
