%ACCOUNTS   National accounting used for
%
%             Sweat Equity in US Private Business
%             by Anmol Bhandari and Ellen McGrattan
%
%           Notes: Run this code on a machine with excel.

%           Ellen McGrattan, 2-20-19
%           Revised, ERM, 5-16-19

%
% User-defined printing options (iprint) and samples (t0:t1)
%   iprint = 0,  print only final results for paper
%          = 1,  print all intermediate results
%
iprint       = 1;
t0           = 2003;
t1           = 2013;
t0           = 2007;
t1           = 2007;
nyr          = t1-t0+1;


%=======================================================================
%
% NIPA Table 1.1.5 Gross Domestic Product (Annual, Mil$)
%
%=======================================================================
tab1_1_5A    = xlsread('../RawData/BEA/NIPA/Section1all_xls','T10105-A');
[n,m]        = size(tab1_1_5A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab1_1_5A    = tab1_1_5A(:,i)/1000;
gdp          = tab1_1_5A(1,:)';
pce          = tab1_1_5A(2,:)';
pced         = tab1_1_5A(4,:)';
pcend        = tab1_1_5A(5,:)';
pces         = tab1_1_5A(6,:)';
gpdi         = tab1_1_5A(7,:)';
gpdinrstr    = tab1_1_5A(10,:)';
gpdinreq     = tab1_1_5A(11,:)';
gpdinripp    = tab1_1_5A(12,:)';
gpdires      = tab1_1_5A(13,:)';
gpdiinv      = tab1_1_5A(14,:)';
netxm        = tab1_1_5A(15,:)';
netx         = tab1_1_5A(16,:)';
netm         = tab1_1_5A(19,:)';
gov          = tab1_1_5A(22,:)';


%=======================================================================
%
% NIPA Table 1.1.9 Implicit price deflators (2012=100)
%
%=======================================================================

tab1_1_9A    = xlsread('../RawData/BEA/NIPA/Section1all_xls','T10109-A');
[n,m]        = size(tab1_1_9A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab1_1_9A    = tab1_1_9A(2:n,i);

dgdp0        = tab1_1_9A(1,:)';
dgdp         = dgdp0/dgdp0(1);


%=======================================================================
%
% NIPA Table 1.7.5 Relation of GDP, GNP, etc. (Annual, Bil$)
%
%=======================================================================

tab1_7_5A    = xlsread('../RawData/BEA/NIPA/Section1all_xls','T10705-A');
[n,m]        = size(tab1_7_5A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab1_7_5A    = tab1_7_5A(:,i);

gnp          = tab1_7_5A(4,:)';
gni          = tab1_7_5A(28,:)';

%=======================================================================
%
% NIPA Table 1.10 Gross domestic income (Annual, Mil$)
%
%=======================================================================

tab1_10A     = xlsread('../RawData/BEA/NIPA/Section1all_xls','T11000-A');
[n,m]        = size(tab1_10A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab1_10A     = tab1_10A(:,i)/1000;

gdi          = tab1_10A(1,:)';
comp         = tab1_10A(2,:)';
ibt          = tab1_10A(7,:)'-tab1_10A(8,:)'+tab1_10A(12,:)';
nint         = tab1_10A(11,:)';
pinc         = tab1_10A(13,:)';
rent         = tab1_10A(14,:)';
prof         = tab1_10A(15,:)';
ctax         = tab1_10A(16,:)';
gprof        = tab1_10A(20,:)';
dep          = tab1_10A(21,:)';
desc         = tab1_10A(25,:)';

%=======================================================================
%
% NIPA Table 1.13  National income by sector (Anuual, Mil$)
%
%=======================================================================

tab1_13A     = xlsread('../RawData/BEA/NIPA/Section1all_xls','T11300-A');
[n,m]        = size(tab1_13A);
t            = 1948:1947+m;
i            = find(t==t0):find(t==t1);
tab1_13A     = tab1_13A(:,i)/1000;

dbus         = tab1_13A(2,:)';
dcbus        = tab1_13A(3,:)';
dcbuscomp    = tab1_13A(4,:)';
dcbusprof    = tab1_13A(7,:)';
dcbusnint    = tab1_13A(8,:)';
dcbusibt     = tab1_13A(9,:)';
dncbus       = tab1_13A(10,:)';
dncbuscomp   = tab1_13A(11,:)';
dncbuspinc   = tab1_13A(14,:)';
dncbusrent   = tab1_13A(15,:)';
dncbusnint   = tab1_13A(16,:)';
dncbusibt    = tab1_13A(17,:)';
dncbusgprof  = tab1_13A(18,:)';
dppbus       = tab1_13A(19,:)';
dppbuscomp   = tab1_13A(20,:)';
dppbuspinc   = tab1_13A(23,:)';
dppbusfpinc  = tab1_13A(24,:)';
dppbusnfpinc = tab1_13A(25,:)';
dppbusnint   = tab1_13A(26,:)';
dppbusibt    = tab1_13A(27,:)';
dopbus       = tab1_13A(28,:)';
dopbuscomp   = tab1_13A(29,:)';
dopbuspinc   = tab1_13A(32,:)';
dopbusrent   = tab1_13A(33,:)';
dopbusnint   = tab1_13A(34,:)';
dopbusibt    = tab1_13A(35,:)';
dgebus       = tab1_13A(36,:)';
dgebuscomp   = tab1_13A(37,:)';
dgebusgprof  = tab1_13A(40,:)';
dhh          = tab1_13A(42,:)';
dhhcomp      = tab1_13A(43,:)';
dhhrent      = tab1_13A(46,:)';
dhhnint      = tab1_13A(47,:)';
dhhibt       = tab1_13A(48,:)';
dnp          = tab1_13A(49,:)';
dnpcomp      = tab1_13A(50,:)';
dnprent      = tab1_13A(53,:)';
dnpnint      = tab1_13A(54,:)';
dnpibt       = tab1_13A(55,:)';
dgg          = tab1_13A(56,:)';
dggcomp      = tab1_13A(57,:)';
row          = tab1_13A(60,:)';
rowcomp      = tab1_13A(61,:)';
rowprof      = tab1_13A(62,:)';
rownint      = tab1_13A(63,:)';
piadj        = sum(tab1_13A([74,77,78],:))';
nbuscomp     = dgebuscomp+dggcomp+dhhcomp+dnpcomp;

%=======================================================================
%
% NIPA Table 5.2.5 Gross and Net Domestic Investment (Annual, Mil$)
%
%=======================================================================

tab5_2_5A    = xlsread('../RawData/BEA/NIPA/Section5all_xls','T50205-A');
[n,m]        = size(tab5_2_5A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab5_2_5A    = tab5_2_5A(:,i)/1000;

depipp       = tab5_2_5A(20,:)';
depres       = tab5_2_5A(23,:)';
depgov       = tab5_2_5A(27,:)';

%=======================================================================
%
% NIPA Table 5.8.5A/B Private inventories and final sales (Quarterly, Bil$)
%
%=======================================================================

tab5_8_5AQ   = xlsread('../RawData/BEA/NIPA/Section5all_xls','T50805A-Q');
tab5_8_5BQ   = xlsread('../RawData/BEA/NIPA/Section5all_xls','T50805B-Q');
[n,m]        = size(tab5_8_5BQ);
tab5_8_5Q    = [tab5_8_5AQ(1:2,:),tab5_8_5BQ(1:2,6:m)];
[n,m]        = size(tab5_8_5Q);
tab5_8_5A    = tab5_8_5Q(:,4:4:m);
[n,m]        = size(tab5_8_5A);
t            = 1947:1946+m;
i            = find(t==t0):find(t==t1);
tab5_8_5A    = tab5_8_5A(:,i)/1000;

kinv         = tab5_8_5A(1,:)';
kcinv        = tab5_8_5A(1,:)'-.9*tab5_8_5A(2,:)';
kncinv       = .9*tab5_8_5A(2,:)';

%=======================================================================
%
% NIPA Table 3.1 Government current receipts and expenditures (Annual, Mil$)
%
%=======================================================================

tab3_1A      = xlsread('../RawData/BEA/NIPA/Section3all_xls','T30100-A');
[n,m]        = size(tab3_1A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab3_1A      = tab3_1A(:,i)/1000;

intrec       = tab3_1A(12,:)';
intpay       = tab3_1A(27,:)';

%=======================================================================
%
% NIPA Table 3.2 Federal govt current receipts and exp. (Annual, Mil$)
%
%=======================================================================

tab3_2A      = xlsread('../RawData/BEA/NIPA/Section3all_xls','T30200-A');
[n,m]        = size(tab3_2A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab3_2A      = tab3_2A(:,i)/1000;

ctaxfed      = tab3_2A(8,:)';

%=======================================================================
%
% NIPA Table 3.3 S&L govt current receipts and exp. (Annual, Mil$)
%
%=======================================================================

tab3_3A      = xlsread('../RawData/BEA/NIPA/Section3all_xls','T30300-A');
[n,m]        = size(tab3_3A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab3_3A      = tab3_3A(:,i)/1000;

ctaxsl       = tab3_3A(11,:)';

%=======================================================================
%
% NIPA Table 3.5 Taxes on production and imports (Annual, Mil$)
%
%=======================================================================

tab3_5A      = xlsread('../RawData/BEA/NIPA/Section3all_xls','T30500-A');
[n,m]        = size(tab3_5A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab3_5A      = tab3_5A(:,i)/1000;

ptax         = tab3_5A(38,:)';
stax         = tab3_5A(1,:)'-ptax;

%=======================================================================
%
% NIPA Table 3.9.5 Govt consumption and gross investment (Annual, Mil$)
%
%=======================================================================

tab3_9_5A    = xlsread('../RawData/BEA/NIPA/Section3all_xls','T30905-A');
[n,m]        = size(tab3_9_5A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab3_9_5A    = tab3_9_5A(:,i)/1000;

govc         = tab3_9_5A(2,:)';
govi         = tab3_9_5A(3,:)';
govist       = tab3_9_5A(4,:)';
govieq       = tab3_9_5A(5,:)';
goviipp      = tab3_9_5A(6,:)';

%=======================================================================
%
% NIPA Table 3.12 Govt social benefits (Annual, Mil$)
%
%=======================================================================

tab3_12A     = xlsread('../RawData/BEA/NIPA/Section3all_xls','T31200-A');
[n,m]        = size(tab3_12A);
t            = 1929:1928+m;
i            = find(t==t0):find(t==t1);
tab3_12A     = tab3_12A(:,i)/1000;

socsec       = tab3_12A(5,:)';
medicare     = tab3_12A(6,:)';
rrret        = tab3_12A(12,:)';
pbg          = tab3_12A(13,:)';
vetpen       = tab3_12A(18,:)';
retfinance   = socsec+medicare+rrret+pbg+vetpen;

%=======================================================================
%
% NIPA Table 6.5 Full-time Equivalent Employees (Annual, Thous.)
%
%=======================================================================

tab6_5CA     = xlsread('../RawData/BEA/NIPA/Section6all_xls','T60500C-A');
tab6_5DA     = xlsread('../RawData/BEA/NIPA/Section6all_xls','T60500D-A');
[n,m]        = size(tab6_5CA);
tab6_5A      = [tab6_5CA(1:3,1:m-3),tab6_5DA(1:3,:)];
[n,m]        = size(tab6_5A);
t            = 1987:1986+m;
i            = find(t==t0):find(t==t1);
tab6_5A      = tab6_5A(:,i)/1000;

ftempl       = tab6_5A(1,:)';
ftempld      = tab6_5A(2,:)';
ftemplp      = tab6_5A(3,:)';
ftemplg      = ftempld-ftemplp;


%=======================================================================
%
% NIPA Table 6.8 Persons Engaged in Production (Annual, Thous.)
%
%=======================================================================

tab6_8CA     = xlsread('../RawData/BEA/NIPA/Section6all_xls','T60800C-A');
tab6_8DA     = xlsread('../RawData/BEA/NIPA/Section6all_xls','T60800D-A');
[n,m]        = size(tab6_8CA);
tab6_8A      = [tab6_8CA(1:3,1:m-3),tab6_8DA(1:3,:)];
[n,m]        = size(tab6_8A);
t            = 1987:1986+m;
i            = find(t==t0):find(t==t1);
tab6_8A      = tab6_8A(:,i)/1000;

pereng       = tab6_8A(1,:)';
perengd      = tab6_8A(2,:)';
perengp      = tab6_8A(3,:)';
perengg      = perengd-perengp;
selfemp      = pereng-ftempl;

%=======================================================================
%
% NIPA Table 6.13 Noncorporate CCA (Annual, Mil$)
%
%=======================================================================

tab6_13CA    = xlsread('../RawData/BEA/NIPA/Section6all_xls','T61300C-A');
tab6_13DA    = xlsread('../RawData/BEA/NIPA/Section6all_xls','T61300D-A');
[n,m]        = size(tab6_13CA);
tab6_13A     = [tab6_13CA(1:3,1:m-3),tab6_13DA(1:3,:)];
[n,m]        = size(tab6_13A);
t            = 1987:1986+m;
i            = find(t==t0):find(t==t1);
tab6_13A     = tab6_13A(:,i)/1000;

nccca        = tab6_13A(1,:)';

%=======================================================================
%
% NIPA Table 6.22 Corporate CCA (Annual, Mil$)
%
%=======================================================================

tab6_22CA    = xlsread('../RawData/BEA/NIPA/Section6all_xls','T62200C-A');
tab6_22DA    = xlsread('../RawData/BEA/NIPA/Section6all_xls','T62200D-A');
[n,m]        = size(tab6_22CA);
tab6_22A     = [tab6_22CA(1:3,1:m-3),tab6_22DA(1:3,:)];
[n,m]        = size(tab6_22A);
t            = 1987:1986+m;
i            = find(t==t0):find(t==t1);
tab6_22A     = tab6_22A(:,i)/1000;

ccca         = tab6_22A(1,:)';

%=======================================================================
%
% FA Table 1.1 Current-cost net stocks (Annual, Bil$, year-end) 
%
%=======================================================================

tab1_1A      = xlsread('../RawData/BEA/FA/Section1all_xls','FAAt101-A');
[n,m]        = size(tab1_1A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab1_1A      = tab1_1A(:,i);

kfix         = tab1_1A(2,:)';
kg           = tab1_1A(9,:)';
kpfix        = kfix-kg;
kcd          = tab1_1A(15,:)';

%=======================================================================
%
% FA Table 1.5 Investments (Annual, Bil$) 
%
%=======================================================================

tab1_5A      = xlsread('../RawData/BEA/FA/Section1all_xls','FAAt105-A');
[n,m]        = size(tab1_5A);
t            = 1901:1900+m;
i            = find(t==t0):find(t==t1);
tab1_5A      = tab1_5A(:,i);

xtot         = tab1_5A(1,:)';
xfix         = tab1_5A(2,:)';
xg           = tab1_5A(9,:)';
xcd          = tab1_5A(15,:)';


%=======================================================================
%
% FA Table 4.1 Current-cost nonresidential net stocks (Annual, Bil$, year-end) 
%
%=======================================================================

tab4_1A      = xlsread('../RawData/BEA/FA/Section4all_xls','FAAt401-A');
[n,m]        = size(tab4_1A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab4_1A      = tab4_1A(:,i);

kpnr         = tab4_1A(1,:)';
kpnreq       = tab4_1A(2,:)';
kpnrstr      = tab4_1A(3,:)';
kpnripp      = tab4_1A(4,:)';
kpnrcs       = tab4_1A(19,:)';
kpnrcseq     = tab4_1A(20,:)';
kpnrcsstr    = tab4_1A(21,:)';
kpnrcsipp    = tab4_1A(22,:)';
kpnrsp       = tab4_1A(63,:)';
kpnrspeq     = tab4_1A(64,:)';
kpnrspstr    = tab4_1A(65,:)';
kpnrspipp    = tab4_1A(66,:)';
kpnrpa       = tab4_1A(67,:)';
kpnrpaeq     = tab4_1A(68,:)';
kpnrpastr    = tab4_1A(69,:)';
kpnrpaipp    = tab4_1A(70,:)';
kpnrnp       = tab4_1A(71,:)';
kpnrnpeq     = tab4_1A(72,:)';
kpnrnpstr    = tab4_1A(73,:)';
kpnrnpipp    = tab4_1A(74,:)';
kpnrhh       = tab4_1A(75,:)';
kpnrhheq     = tab4_1A(76,:)';
kpnrhhstr    = tab4_1A(77,:)';
kpnrhhipp    = tab4_1A(78,:)';
kpnrte       = tab4_1A(79,:)';
kpnrteeq     = tab4_1A(80,:)';
kpnrtestr    = tab4_1A(81,:)';
kpnrteipp    = tab4_1A(82,:)';

%=======================================================================
%
% FA Table 4.4 Current-cost depreciation (Annual, Bil$, year-end) 
%
%=======================================================================

tab4_4A      = xlsread('../RawData/BEA/FA/Section4all_xls','FAAt404-A');
[n,m]        = size(tab4_4A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab4_4A      = tab4_4A(:,i);

dpnr         = tab4_4A(1,:)';
dpnreq       = tab4_4A(2,:)';
dpnrstr      = tab4_4A(3,:)';
dpnripp      = tab4_4A(4,:)';
dpnrcs       = tab4_4A(19,:)';
dpnrcseq     = tab4_4A(20,:)';
dpnrcsstr    = tab4_4A(21,:)';
dpnrcsipp    = tab4_4A(22,:)';
dpnrsp       = tab4_4A(63,:)';
dpnrspeq     = tab4_4A(64,:)';
dpnrspstr    = tab4_4A(65,:)';
dpnrspipp    = tab4_4A(66,:)';
dpnrpa       = tab4_4A(67,:)';
dpnrpaeq     = tab4_4A(68,:)';
dpnrpastr    = tab4_4A(69,:)';
dpnrpaipp    = tab4_4A(70,:)';
dpnrnp       = tab4_4A(71,:)';
dpnrnpeq     = tab4_4A(72,:)';
dpnrnpstr    = tab4_4A(73,:)';
dpnrnpipp    = tab4_4A(74,:)';
dpnrhh       = tab4_4A(75,:)';
dpnrhheq     = tab4_4A(76,:)';
dpnrhhstr    = tab4_4A(77,:)';
dpnrhhipp    = tab4_4A(78,:)';
dpnrte       = tab4_4A(79,:)';
dpnrteeq     = tab4_4A(80,:)';
dpnrtestr    = tab4_4A(81,:)';
dpnrteipp    = tab4_4A(82,:)';


%=======================================================================
%
% FA Table 4.6 Historical-cost depreciation (Annual, Bil$, year-end) 
%
%=======================================================================

tab4_6A      = xlsread('../RawData/BEA/FA/Section4all_xls','FAAt406-A');
[n,m]        = size(tab4_6A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab4_6A      = tab4_6A(:,i);

hdpnr        = tab4_6A(1,:)';
hdpnreq      = tab4_6A(2,:)';
hdpnrstr     = tab4_6A(3,:)';
hdpnripp     = tab4_6A(4,:)';
hdpnrcs      = tab4_6A(19,:)';
hdpnrcseq    = tab4_6A(20,:)';
hdpnrcsstr   = tab4_6A(21,:)';
hdpnrcsipp   = tab4_6A(22,:)';
hdpnrsp      = tab4_6A(63,:)';
hdpnrspeq    = tab4_6A(64,:)';
hdpnrspstr   = tab4_6A(65,:)';
hdpnrspipp   = tab4_6A(66,:)';
hdpnrpa      = tab4_6A(67,:)';
hdpnrpaeq    = tab4_6A(68,:)';
hdpnrpastr   = tab4_6A(69,:)';
hdpnrpaipp   = tab4_6A(70,:)';
hdpnrnp      = tab4_6A(71,:)';
hdpnrnpeq    = tab4_6A(72,:)';
hdpnrnpstr   = tab4_6A(73,:)';
hdpnrnpipp   = tab4_6A(74,:)';
hdpnrhh      = tab4_6A(75,:)';
hdpnrhheq    = tab4_6A(76,:)';
hdpnrhhstr   = tab4_6A(77,:)';
hdpnrhhipp   = tab4_6A(78,:)';
hdpnrte      = tab4_6A(79,:)';
hdpnrteeq    = tab4_6A(80,:)';
hdpnrtestr   = tab4_6A(81,:)';
hdpnrteipp   = tab4_6A(82,:)';

%=======================================================================
%
% FA Table 4.7 Investment, private nonresidential (Annual, Bil$) 
%
%=======================================================================

tab4_7A      = xlsread('../RawData/BEA/FA/Section4all_xls','FAAt407-A');
[n,m]        = size(tab4_7A);
t            = 1901:1900+m;
i            = find(t==t0):find(t==t1);
tab4_7A      = tab4_7A(:,i);

xpnr         = tab4_7A(1,:)';
xpnreq       = tab4_7A(2,:)';
xpnrstr      = tab4_7A(3,:)';
xpnripp      = tab4_7A(4,:)';
xpnrcs       = tab4_7A(19,:)';
xpnrcseq     = tab4_7A(20,:)';
xpnrcsstr    = tab4_7A(21,:)';
xpnrcsipp    = tab4_7A(22,:)';
xpnrsp       = tab4_7A(63,:)';
xpnrspeq     = tab4_7A(64,:)';
xpnrspstr    = tab4_7A(65,:)';
xpnrspipp    = tab4_7A(66,:)';
xpnrpa       = tab4_7A(67,:)';
xpnrpaeq     = tab4_7A(68,:)';
xpnrpastr    = tab4_7A(69,:)';
xpnrpaipp    = tab4_7A(70,:)';
xpnrnp       = tab4_7A(71,:)';
xpnrnpeq     = tab4_7A(72,:)';
xpnrnpstr    = tab4_7A(73,:)';
xpnrnpipp    = tab4_7A(74,:)';
xpnrhh       = tab4_7A(75,:)';
xpnrhheq     = tab4_7A(76,:)';
xpnrhhstr    = tab4_7A(77,:)';
xpnrhhipp    = tab4_7A(78,:)';
xpnrte       = tab4_7A(79,:)';
xpnrteeq     = tab4_7A(80,:)';
xpnrtestr    = tab4_7A(81,:)';
xpnrteipp    = tab4_7A(82,:)';


%=======================================================================
%
% FA Table 5.1 Current-cost residential net stocks (Annual, Bil$, year-end) 
%
%=======================================================================

tab5_1A      = xlsread('../RawData/BEA/FA/Section5all_xls','FAAt501-A');
[n,m]        = size(tab5_1A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab5_1A      = tab5_1A(:,i);

kres         = tab5_1A(1,:)';
kpres        = tab5_1A(3,:)';
kprescs      = tab5_1A(4,:)';
kpressp      = tab5_1A(6,:)';
kpresnp      = tab5_1A(7,:)';
kpreshh      = tab5_1A(8,:)';
kgres        = tab5_1A(9,:)';

%=======================================================================
%
% FA Table 5.4 Current-cost depreciation of residential (Annual, Bil$, year-end) 
%
%=======================================================================

tab5_4A      = xlsread('../RawData/BEA/FA/Section5all_xls','FAAt504-A');
[n,m]        = size(tab5_4A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab5_4A      = tab5_4A(:,i);

dres         = tab5_4A(1,:)';
dpres        = tab5_4A(3,:)';
dprescs      = tab5_4A(4,:)';
dpressp      = tab5_4A(6,:)';
dpresnp      = tab5_4A(7,:)';
dpreshh      = tab5_4A(8,:)';
dgres        = tab5_4A(9,:)';

%=======================================================================
%
% FA Table 5.6 Historical-cost depreciation of residential (Annual, Bil$, year-end) 
%
%=======================================================================

tab5_6A      = xlsread('../RawData/BEA/FA/Section5all_xls','FAAt506-A');
[n,m]        = size(tab5_6A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab5_6A      = tab5_6A(:,i);

hdres        = tab5_6A(1,:)';
hdpres       = tab5_6A(3,:)';
hdprescs     = tab5_6A(4,:)';
hdpressp     = tab5_6A(6,:)';
hdpresnp     = tab5_6A(7,:)';
hdpreshh     = tab5_6A(8,:)';
hdgres       = tab5_6A(9,:)';

%=======================================================================
%
% FA Table 5.7 Investment, residential (Annual, Bil$) 
%
%=======================================================================

tab5_7A      = xlsread('../RawData/BEA/FA/Section5all_xls','FAAt507-A');
[n,m]        = size(tab5_7A);
t            = 1901:1900+m;
i            = find(t==t0):find(t==t1);
tab5_7A      = tab5_7A(:,i);

xres         = tab5_7A(1,:)';
xpres        = tab5_7A(3,:)';
xprescs      = tab5_7A(4,:)';
xpressp      = tab5_7A(6,:)';
xpresnp      = tab5_7A(7,:)';
xpreshh      = tab5_7A(8,:)';
xgres        = tab5_7A(9,:)';


%=======================================================================
%
% FA Table 8.4 Current-cost Durable Depreciation (Annual, $Bil, year-end) 
%
%=======================================================================

tab8_4A      = xlsread('../RawData/BEA/FA/Section8all_xls','FAAt804-A');
[n,m]        = size(tab8_4A);
t            = 1925:1924+m;
i            = find(t==t0):find(t==t1);
tab8_4A      = tab8_4A(:,i);

cddep        = tab8_4A(1,:)';


%=======================================================================
%
% FOF Table B.101. Balance sheet of Households (Annual, $Mil)
%
%=======================================================================
btab_101Q    = xlsread('../RawData/FOF/xlsx/b101');
[n,m]        = size(btab_101Q);
btab_101A    = btab_101Q(1:4:n,:)';
[n,m]        = size(btab_101A);
t            = 1945:1944+m;
i            = find(t==t0):find(t==t1);
btab_101A    = btab_101A(:,i)/1000;

hhland       = btab_101A(3,:)'-sum(btab_101A([45,48],:))';

%=======================================================================
%
% FOF Table B.103. Balance sheet of Nonfinancial Corporations (Annual, $Mil)
%
%=======================================================================

btab_103Q    = xlsread('../RawData/FOF/xlsx/b103');
[n,m]        = size(btab_103Q);
btab_103A    = btab_103Q(1:4:n,:)';
[n,m]        = size(btab_103A);
t            = 1945:1944+m;
i            = find(t==t0):find(t==t1);
btab_103A    = btab_103A(:,i)/1000;

cland        = btab_103A(3,:)'-sum(btab_103A([40:41],:))';

%=======================================================================
%
% FOF Table B.104. Balance sheet of Nonfinancial Noncorp Bus (Annual, $Mil)
%
%=======================================================================

btab_104Q    = xlsread('../RawData/FOF/xlsx/b104');
[n,m]        = size(btab_104Q);
btab_104A    = btab_104Q(1:4:n,:)';
[n,m]        = size(btab_104A);
t            = 1945:1944+m;
i            = find(t==t0):find(t==t1);
btab_104A    = btab_104A(:,i)/1000;

ncland       = btab_104A(3,:)'-sum(btab_104A([36:37],:))';


%=======================================================================
%
% IRS Returns of Active Corporations (All)
%     Table 2. Balance Sheets and Income Statements ($Thous)
%
%=======================================================================

load ../RawData/IRS/co02ccr/co02ccr.dat
t            = co02ccr(:,1);
i            = find(t==t0):find(t==t1);

asset        = co02ccr(i,2);
cogs         = co02ccr(i,3);
coff         = co02ccr(i,4);
salw         = co02ccr(i,5);
labc         = co02ccr(i,6);
pena         = co02ccr(i,7);
empb         = co02ccr(i,8);


%=======================================================================
%
% IRS Returns of Active S Corporations 
%     Table 7-S. Balance Sheets and Income Statements ($Thous)
%
%=======================================================================

load ../RawData/IRS/co07s/co07sall.dat
load ../RawData/IRS/co07s/co07sre.dat
t            = co07sall(:,1);
i            = find(t==t0):find(t==t1);

scasset      = co07sall(i,2);
sccogs       = co07sall(i,3);
sccomp       = co07sall(i,4);
scsalw       = co07sall(i,5);
scpena       = co07sall(i,6);
scempb       = co07sall(i,7);
scdepl       = co07sall(i,8);
scdepr       = co07sall(i,9);
scamor       = co07sall(i,10);
scintp       = co07sall(i,11);
scrent       = co07sall(i,12);
scneti       = co07sall(i,13);

scdeprr      = co07sre(i,9);
scamorr      = co07sre(i,10);
scintpr      = co07sre(i,11);
screntr      = co07sre(i,12);
scnetir      = co07sre(i,13);


%=======================================================================
%
% IRS Returns of Active Partnerships
%     Table 1. Total Assets, Trade or Bus Income and Deductions
%
%=======================================================================

load ../RawData/IRS/pa01/pa01all.dat
load ../RawData/IRS/pa01/pa01re.dat
t            = pa01all(:,1);
i            = find(t==t0):find(t==t1);

padepl       = pa01all(i,2);
padepr       = pa01all(i,3);
paintp       = pa01all(i,4);
parent       = pa01all(i,5);
paneti       = pa01all(i,6);
pareri       = pa01all(i,7);
pareot       = pa01all(i,8);

padeprr      = pa01re(i,3);
paintpr      = pa01re(i,4);
parentr      = pa01re(i,5);
panetir      = pa01re(i,6);
parerir      = pa01re(i,7);
pareotr      = pa01re(i,8);


%=======================================================================
%
% IRS Returns of Nonfarm Sole Proprietorships
%     Table 1. Business receipts, selected deductions, and net income
%
%=======================================================================

load ../RawData/IRS/sp02is/sp02isall.dat
load ../RawData/IRS/sp02is/sp02isre.dat

t            = sp02isall(:,1);
i            = find(t==t0):find(t==t1);

spdepl       = sp02isall(i,2);
spdepr       = sp02isall(i,3);
spmint       = sp02isall(i,4);
spoint       = sp02isall(i,5);
sprpme       = sp02isall(i,6);
sprpbp       = sp02isall(i,7);
sphoff       = sp02isall(i,8);
spneti       = sp02isall(i,9);

spdeprr      = sp02isre(i,3);
spmintr      = sp02isre(i,4);
spointr      = sp02isre(i,5);
sprpmer      = sp02isre(i,6);
sprpbpr      = sp02isre(i,7);
sphoffr      = sp02isre(i,8);
spnetir      = sp02isre(i,9);

ptrent       = scrent+parent+sprpme+sprpbp-(screntr+parentr+sprpmer+sprpbpr);
ptintp       = scintp+paintp+spmint+spoint-(scintpr+paintpr+spmintr+spointr);
ptdepr       = scamor+scdepr+scdepl+padepr+padepl+spdepr+spdepl- ...
               (scamorr+scdeprr+padeprr+spdeprr);

%=======================================================================
%
% IRS federal corporate income tax rates
%
%=======================================================================

load ../RawData/IRS/citrate/citrate.dat
t            = citrate(:,1);
i            = find(t==t0):find(t==t1);

domrate      = citrate(i,2).*((ctaxfed+ctaxsl)./ctaxfed);
domrate2     = citrate(i,2)*(1-.09/3).*((ctaxfed+ctaxsl)./ctaxfed);

%=======================================================================
%
% KPMG global corporate income tax rates
%
%=======================================================================

di           = xlsread('../RawData/IRS/citrate/di');
[n,m]        = size(di);
t            = 1982:1981+m;
i            = find(t==t0):find(t==t1);
di           = di([3,27,64],i);

load ../RawData/IRS/citrate/kpmg.dat
t            = kpmg(1,:);
if ((t0<2003) & (t1<2003));
  i          = find(t==t0):find(t==t1);
  forrate    = ones(3,1)*kpmg(7,i);
end; 
if ((t0<2003) & (t1>=2003));
  i          = find(t==t0):find(t==2002);
  forrate    = ones(3,1)*kpmg(7,i);
  i          = find(t==2003):find(t==t1);
  forrate    = [forrate,kpmg([6,8,4],i)];
end; 
if ((t0>=2003) & (t1>=2003));
  i          = find(t==t0):find(t==t1);
  forrate    = kpmg([6,8,4],i);
end;
forrate      = (sum((di.*forrate)./sum(di)))';

txrate       = (domrate.*(prof-scneti)+forrate.*rowprof)./ ...
               (prof+rowprof-scneti);
txrate2      = (domrate2.*(prof-scneti)+forrate.*rowprof)./ ...
               (prof+rowprof-scneti);


%=======================================================================
%
% BLS-CPS Average hours per person and population 16-64 (totals)
%
%=======================================================================

cpu          = xlsread('../RawData/BLS/CPU_Data_construction','BLS-CPS');
cpu          = cpu(3:72,:);
t            = 1948:2017;
i            = find(t==t0):find(t==t1);

blsahr       = cpu(i,22);
pop1664      = cpu(i,20)/10^6;

%=======================================================================
%
% Census population by age (collected by OECDStats, 1950-2019)
%   All ages and 18 5-year categories (with 85+ last)
%
%=======================================================================

popbyage     = xlsread('../RawData/Census/PopByAge/OECDpopbyage');
t            = popbyage(1:70,1);
i            = find(t==t0):find(t==t1);
i07          = find(t==2007);

page         = zeros(70,19);
page(:)      = popbyage(:,2);
page         = [sum(page(:,2:5)')',page(:,6)+page(:,7),page(:,8)+ ...
                page(:,9),page(:,10)+page(:,11),page(:,12)+page(:,13), ...
                sum(page(:,14:19)')'];
page07       = page(i07,:)/sum(page(i07,:))*100;
page         = page(i,:)./(sum(page(i,:)')'*ones(1,6))*100;
page         = page(:);


%=======================================================================
%
% Survey of Business Owners 2007
%   Number of businesses, owners, hours
%
%=======================================================================

load CBO/2007Survey/pums2.out
wgt         = pums2(:,1);
numown      = pums2(:,2);
acqyr1      = pums2(:,3);
acqyr2      = pums2(:,4);
acqyr3      = pums2(:,5);
acqyr4      = pums2(:,6);
scsavings   = pums2(:,7);
scassets    = pums2(:,8);
scequity    = pums2(:,9);
sccredit    = pums2(:,10);
scgovtloan  = pums2(:,11);
scgovtguar  = pums2(:,12);
scbankloan  = pums2(:,13);
scfamloan   = pums2(:,14);
scventure   = pums2(:,15);
scgrant     = pums2(:,16);
scother     = pums2(:,17);
scdontknow  = pums2(:,18);
scnotneed   = pums2(:,19);
scnotrep    = pums2(:,20);
scamount    = pums2(:,21);
prminc1     = pums2(:,22);
prminc2     = pums2(:,23);
prminc3     = pums2(:,24);
prminc4     = pums2(:,25);
pct1        = pums2(:,26);
pct2        = pums2(:,27);
pct3        = pums2(:,28);
pct4        = pums2(:,29);
hours1      = pums2(:,30);
hours2      = pums2(:,31);
hours3      = pums2(:,32);
hours4      = pums2(:,33);
age1        = pums2(:,34);
age2        = pums2(:,35);
age3        = pums2(:,36);
age4        = pums2(:,37);
estab       = pums2(:,38);
status      = pums2(:,39);
payroll     = pums2(:,40);
operate     = pums2(:,41);
sold        = pums2(:,42);
retired     = pums2(:,43);
dead        = pums2(:,44);

%
% Number of businesses, reporting/not reporting info on owners
%
i           = find(numown>0 & numown<8);
j           = find(numown==0  | numown==8 | isnan(numown));
nbus        = sum(wgt)/1000;
nbrepo      = sum(wgt(i))/1000;
nbnrepo     = sum(wgt(j))/1000;

%
% Owner counts (thousands) and hours, 
%   Totals & owners with business as main income
%
hvec        = [0,10,30,40,50,70];

%
%  -- 1-owner businesses, all owners
%
i           = find(numown==1);
j           = find((numown==0 | numown==8 | isnan(numown)) & pct2==0);
k           = find(numown==1 & hours1>0); 
nbrep1o     = sum(wgt(i))/1000;
nbnrep1o    = sum(wgt(j))/1000;
sum1        = 0;
for l=2:6;
  m         = find(numown==1 & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs1o      = sum1/sum(wgt(k));
town1o      = nbrep1o+nbnrep1o;

%
%  -- 1-owner businesses, owners with business as primary income
%
i           = find(numown==1 & prminc1==1);
k           = find(numown==1 & prminc1==1 & hours1>0); 
norep1omi   = sum(wgt(i))/1000;
nonrep1omi  = norep1omi/nbrep1o * nbnrep1o;
sum1        = 0;
for l=2:6;
  m         = find(numown==1 & prminc1==1 & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs1omi    = sum1/sum(wgt(k));
town1omi    = norep1omi+nonrep1omi;

%
%  -- 2-owner businesses, all owners
%
i           = find(numown==2);
j           = find((numown==0 | numown==8 | isnan(numown)) & pct2>0 & pct3==0);
k1          = find(numown==2 & hours1>0); 
k2          = find(numown==2 & hours2>0); 
nbrep2o     = sum(wgt(i))/1000;
nbnrep2o    = sum(wgt(j))/1000;
sum1        = 0;
sum2        = 0;
for l=2:6;
  m1        = find(numown==2 & hours1==l);
  m2        = find(numown==2 & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m1));
  sum2      = sum2+hvec(l)*sum(wgt(m2));
end;
ahrs2o1     = sum1/sum(wgt(k1));
ahrs2o2     = sum2/sum(wgt(k2));
ahrs2o      = (ahrs2o1+ahrs2o2)/2;
town2o      = 2*(nbrep2o+nbnrep2o);

%
%  -- 2-owner businesses, owners with business as primary income
%
i1          = find(numown==2 & prminc1==1);
i2          = find(numown==2 & prminc2==1);
k1          = find(numown==2 & prminc1==1 & hours1>0); 
k2          = find(numown==2 & prminc2==1 & hours2>0); 
norep2omi   = (sum(wgt(i1))+sum(wgt(i2)))/1000;
nonrep2omi  = norep2omi/nbrep2o * nbnrep2o;
sum1        = 0;
sum2        = 0;
for l=2:6;
  m1        = find(numown==2 & prminc1==1 & hours1==l);
  m2        = find(numown==2 & prminc2==1 & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m1));
  sum2      = sum2+hvec(l)*sum(wgt(m2));
end;
ahrs2o1mi   = sum1/sum(wgt(k1));
ahrs2o2mi   = sum2/sum(wgt(k2));
ahrs2omi    = (ahrs2o1mi+ahrs2o2mi)/2;
town2omi    = norep2omi+nonrep2omi;

%
%  -- 3-owner businesses, all owners
%
i           = find(numown==3);
j           = find((numown==0 | numown==8 | isnan(numown)) & pct3>0 & pct4==0);
k1          = find(numown==3 & hours1>0); 
k2          = find(numown==3 & hours2>0); 
k3          = find(numown==3 & hours3>0); 
nbrep3o     = sum(wgt(i))/1000;
nbnrep3o    = sum(wgt(j))/1000;
sum1        = 0;
sum2        = 0;
sum3        = 0;
for l=2:6;
  m1        = find(numown==3 & hours1==l);
  m2        = find(numown==3 & hours2==l);
  m3        = find(numown==3 & hours3==l);
  sum1      = sum1+hvec(l)*sum(wgt(m1));
  sum2      = sum2+hvec(l)*sum(wgt(m2));
  sum3      = sum3+hvec(l)*sum(wgt(m3));
end;
ahrs3o1     = sum1/sum(wgt(k1));
ahrs3o2     = sum2/sum(wgt(k2));
ahrs3o3     = sum3/sum(wgt(k3));
ahrs3o      = (ahrs3o1+ahrs3o2+ahrs3o3)/3;
town3o      = 3*(nbrep3o+nbnrep3o);

%
%  -- 3-owner businesses, owners with business as primary income
%
i1          = find(numown==3 & prminc1==1);
i2          = find(numown==3 & prminc2==1);
i3          = find(numown==3 & prminc3==1);
k1          = find(numown==3 & prminc1==1 & hours1>0); 
k2          = find(numown==3 & prminc2==1 & hours2>0); 
k3          = find(numown==3 & prminc3==1 & hours3>0); 
norep3omi   = (sum(wgt(i1))+sum(wgt(i2))+sum(wgt(i3)))/1000;
nonrep3omi  = norep3omi/nbrep3o * nbnrep3o;
sum1        = 0;
sum2        = 0;
sum3        = 0;
for l=2:6;
  m1        = find(numown==3 & prminc1==1 & hours1==l);
  m2        = find(numown==3 & prminc2==1 & hours2==l);
  m3        = find(numown==3 & prminc3==1 & hours3==l);
  sum1      = sum1+hvec(l)*sum(wgt(m1));
  sum2      = sum2+hvec(l)*sum(wgt(m2));
  sum3      = sum3+hvec(l)*sum(wgt(m3));
end;
ahrs3o1mi   = sum1/sum(wgt(k1));
ahrs3o2mi   = sum2/sum(wgt(k2));
ahrs3o3mi   = sum3/sum(wgt(k3));
ahrs3omi    = (ahrs3o1mi+ahrs3o2mi+ahrs3o3mi)/3;
town3omi    = norep3omi+nonrep3omi;


%
%  -- 4-owner businesses, all owners
%
i           = find(numown==4);
j           = find((numown==0 | numown==8 | isnan(numown)) & pct4>0);
k1          = find(numown==4 & hours1>0); 
k2          = find(numown==4 & hours2>0); 
k3          = find(numown==4 & hours3>0); 
k4          = find(numown==4 & hours4>0); 
nbrep4o     = sum(wgt(i))/1000;
nbnrep4o    = sum(wgt(j))/1000;
sum1        = 0;
sum2        = 0;
sum3        = 0;
sum4        = 0;
for l=2:6;
  m1        = find(numown==4 & hours1==l);
  m2        = find(numown==4 & hours2==l);
  m3        = find(numown==4 & hours3==l);
  m4        = find(numown==4 & hours4==l);
  sum1      = sum1+hvec(l)*sum(wgt(m1));
  sum2      = sum2+hvec(l)*sum(wgt(m2));
  sum3      = sum3+hvec(l)*sum(wgt(m3));
  sum4      = sum3+hvec(l)*sum(wgt(m4));
end;
ahrs4o1     = sum1/sum(wgt(k1));
ahrs4o2     = sum2/sum(wgt(k2));
ahrs4o3     = sum3/sum(wgt(k3));
ahrs4o4     = sum4/sum(wgt(k4));
ahrs4o      = (ahrs4o1+ahrs4o2+ahrs4o3+ahrs4o4)/4;
town4o      = 4*(nbrep4o+nbnrep4o);

%
%  -- 4-owner businesses, owners with business as primary income
%
i1          = find(numown==4 & prminc1==1);
i2          = find(numown==4 & prminc2==1);
i3          = find(numown==4 & prminc3==1);
i4          = find(numown==4 & prminc4==1);
k1          = find(numown==4 & prminc1==1 & hours1>0); 
k2          = find(numown==4 & prminc2==1 & hours2>0); 
k3          = find(numown==4 & prminc3==1 & hours3>0); 
k4          = find(numown==4 & prminc4==1 & hours4>0); 
norep4omi   = (sum(wgt(i1))+sum(wgt(i2))+sum(wgt(i3))+sum(wgt(i4)))/1000;
nonrep4omi  = norep4omi/nbrep4o * nbnrep4o;
sum1        = 0;
sum2        = 0;
sum3        = 0;
sum4        = 0;
for l=2:6;
  m1        = find(numown==4 & prminc1==1 & hours1==l);
  m2        = find(numown==4 & prminc2==1 & hours2==l);
  m3        = find(numown==4 & prminc3==1 & hours3==l);
  m4        = find(numown==4 & prminc4==1 & hours4==l);
  sum1      = sum1+hvec(l)*sum(wgt(m1));
  sum2      = sum2+hvec(l)*sum(wgt(m2));
  sum3      = sum3+hvec(l)*sum(wgt(m3));
  sum4      = sum4+hvec(l)*sum(wgt(m4));
end;
ahrs4o1mi   = sum1/sum(wgt(k1));
ahrs4o2mi   = sum2/sum(wgt(k2));
ahrs4o3mi   = sum3/sum(wgt(k3));
ahrs4o4mi   = sum4/sum(wgt(k4));
ahrs4omi    = (ahrs4o1mi+ahrs4o2mi+ahrs4o3mi+ahrs4o4mi)/4;
town4omi    = norep4omi+nonrep4omi;

town        = town1o+town2o+town3o+town4o;
townmi      = town1omi+town2omi+town3omi+town4omi;
ahrs        = (ahrs1o*town1o+ahrs2o*town2o+ahrs3o*town3o+ahrs4o*town4o)/town;
ahrsmi      = (ahrs1omi*town1omi+ahrs2omi*town2omi+ahrs3omi*town3omi+ ...
              ahrs4omi*town4omi)/townmi;

%=======================================================================
%
% Survey of Business Owners 2007
%   Business acquisition profiles
%
%=======================================================================

acqyrs1    = [0,1,2,3,8,18,28];
acqyrs2    = [0,1,2,7,17,27,47];
for m=7:-1:1;
  i        = find(acqyr1==m);
  j        = find(acqyr2==m);
  k        = find(acqyr3==m);
  l        = find(acqyr4==m);
  s(8-m)   = (sum(wgt(i))+sum(wgt(j))+sum(wgt(k))+sum(wgt(l)))/1000;
end;
acqpct     = [s(1:3),s(4)/5,s(5:6)/10,s(7)/20]/sum(s)*100;

for m=7:-1:1;
  i        = find(prminc1==1 & acqyr1==m);
  j        = find(prminc2==1 & acqyr2==m);
  k        = find(prminc3==1 & acqyr3==m);
  l        = find(prminc4==1 & acqyr4==m);
  s(8-m)   = (sum(wgt(i))+sum(wgt(j))+sum(wgt(k))+sum(wgt(l)))/1000;
end;
acqpctmi   = [s(1:3),s(4)/5,s(5:6)/10,s(7)/20]/sum(s)*100;

%=======================================================================
%
% Survey of Business Owners 2007
%   Owner age vs Business Age (since establishment or acquisition)
%
%=======================================================================

estyrs1    = [0,1,2,3,4,5,8,18,28];
estyrs2    = [0,1,2,3,4,7,17,27,47];
estyrs     = .5*estyrs1+.5*estyrs2;
s          = 0;
k1         = find(estab>=1 & estab<=9 & age1>=2 & age1<=6);
k2         = find(estab>=1 & estab<=9 & age2>=2 & age2<=6);
k3         = find(estab>=1 & estab<=9 & age3>=2 & age3<=6);
k4         = find(estab>=1 & estab<=9 & age4>=2 & age4<=6);
for l=2:6;
  i1       = find(estab>=1 & estab<=9 & age1==l);
  i2       = find(estab>=1 & estab<=9 & age2==l);
  i3       = find(estab>=1 & estab<=9 & age3==l);
  i4       = find(estab>=1 & estab<=9 & age4==l);
  for m=9:-1:1;
    j1     = find(estab==m & age1==l);
    j2     = find(estab==m & age2==l);
    j3     = find(estab==m & age3==l);
    j4     = find(estab==m & age4==l);
    s(10-m)= (sum(wgt(j1))+sum(wgt(j2))+sum(wgt(j3))+sum(wgt(j4)))/ ...
             (sum(wgt(i1))+sum(wgt(i2))+sum(wgt(i3))+sum(wgt(i4)))*100;
  end;
  bage(l)  = sum(estyrs.*s)/sum(s);
  j1       = find(estab>=1 & estab<=9 & age1==l);
  j2       = find(estab>=1 & estab<=9 & age2==l);
  j3       = find(estab>=1 & estab<=9 & age3==l);
  j4       = find(estab>=1 & estab<=9 & age4==l);
  nage(l)  = (sum(wgt(j1))+sum(wgt(j2))+sum(wgt(j3))+sum(wgt(j4)))/ ...
             (sum(wgt(k1))+sum(wgt(k2))+sum(wgt(k3))+sum(wgt(k4)))*100;
end;    
bage       = bage(:);
nage       = nage(:);


acqyrs1    = [0,1,2,3,8,18,28];
acqyrs2    = [0,1,2,7,17,27,47];
acqyrs     = .5*acqyrs1+.5*acqyrs2;
s          = 0;
k1         = find(acqyr1>=1 & acqyr1<=7 & age1>=2 & age1<=6);
k2         = find(acqyr2>=1 & acqyr2<=7 & age2>=2 & age2<=6);
k3         = find(acqyr3>=1 & acqyr3<=7 & age3>=2 & age3<=6);
k4         = find(acqyr4>=1 & acqyr4<=7 & age4>=2 & age4<=6);
for l=2:6;
  i1       = find(acqyr1>=1 & acqyr1<=7 & age1==l);
  i2       = find(acqyr2>=1 & acqyr2<=7 & age2==l);
  i3       = find(acqyr3>=1 & acqyr3<=7 & age3==l);
  i4       = find(acqyr4>=1 & acqyr4<=7 & age4==l);
  for m=7:-1:1;
    j1     = find(acqyr1==m & age1==l);
    j2     = find(acqyr2==m & age2==l);
    j3     = find(acqyr3==m & age3==l);
    j4     = find(acqyr4==m & age4==l);
    s(8-m) = (sum(wgt(j1))+sum(wgt(j2))+sum(wgt(j3))+sum(wgt(j4)))/ ...
             (sum(wgt(i1))+sum(wgt(i2))+sum(wgt(i3))+sum(wgt(i4)))*100;
  end;
  bage2(l) = sum(acqyrs.*s)/sum(s);
  j1       = find(acqyr1>=1 & acqyr1<=7 & age1==l);
  j2       = find(acqyr2>=1 & acqyr2<=7 & age2==l);
  j3       = find(acqyr3>=1 & acqyr3<=7 & age3==l);
  j4       = find(acqyr4>=1 & acqyr4<=7 & age4==l);
  nage2(l) = (sum(wgt(j1))+sum(wgt(j2))+sum(wgt(j3))+sum(wgt(j4)))/ ...
             (sum(wgt(k1))+sum(wgt(k2))+sum(wgt(k3))+sum(wgt(k4)))*100;
end;    
bage2      = bage2(:);
nage2      = nage2(:);

%=======================================================================
%
% Survey of Business Owners 2007
%   Sources of capital 
%     -- Percent using bank, venture, or govt loans and grants
%     -- Amounts in intervals [0,1mil+],[0,25K],[25K,100K],[100K+]
%
%=======================================================================

tem        = (min(pums2(:,[7:17,19])')==1)';
i          = find(tem==1);
j          = find((scbankloan==1 | scventure==1 | scgovtloan==1 |  ...
                  scgovtguar==1 | scgrant==1) & tem==1);                                   
extcap     = sum(wgt(j))/sum(wgt(i))*100;

%
%   -- percent with startup capital in [0,1mil or more]
%
i          = find(tem==1 & scamount>0 & scamount<9);
j          = find((scbankloan==1 | scventure==1 | scgovtloan==1 |  ...
             scgovtguar==1 | scgrant==1) & tem==1 & scamount>0 & scamount<9);
k          = find((scbankloan==1 | scventure==1 | scgovtloan==1 |  ...
             scgovtguar==1 | scgrant==1) & tem==1 & scamount>0 & scamount<=3);
l          = find((scbankloan==1 | scventure==1 | scgovtloan==1 |  ...
             scgovtguar==1 | scgrant==1) & tem==1 & scamount>3 & scamount<=5);
m          = find((scbankloan==1 | scventure==1 | scgovtloan==1 |  ...
             scgovtguar==1 | scgrant==1) & tem==1 & scamount>5 & scamount<9);

amtcap0    = sum(wgt(j))/sum(wgt(i))*100;
amtcap1    = sum(wgt(k))/sum(wgt(i))*100;
amtcap2    = sum(wgt(l))/sum(wgt(i))*100;
amtcap3    = sum(wgt(m))/sum(wgt(i))*100;



%=======================================================================
%
% Survey of Business Owners 2007
%   Concentration of business payroll
%     -- Conditioning on all businesses with < 10 owners
%     -- Sort payroll per owner
%     -- Compute % of payroll overseen by bottom 80% of owners
%
%=======================================================================

nowner     = 0*numown;
for i=1:length(numown);
  if (numown(i)==5);
    nowner(i,1) = 7;
  end;
  if (numown(i)==6);
    nowner(i,1) = 25;
  end;
  if (numown(i)==7);
    nowner(i,1) = 75;
  end;
end;
tab        = [wgt,nowner,payroll];

%
% -- All businesses
%
i          = find(nowner>=1 & nowner<=9);
tab1       = tab(i,:);
w          = tab1(:,1);
n          = tab1(:,2);
p          = tab1(:,3);
[pn,ind]   = sort(p./n);
wgti       = w(ind);
nowni      = n(ind);
payi       = p(ind);
pop        = [0;cumsum(wgti.*nowni)/sum(wgti.*nowni)];
pay        = [0;cumsum(wgti.*payi)/sum(wgti.*payi)];
i          = find(pop>=.8);
pay80      = pay(i(1))*100;


%
% -- All employer businesses
%
i          = find(nowner>=1 & nowner<=9 & status==1);
tab2       = tab(i,:);
w          = tab2(:,1);
n          = tab2(:,2);
p          = tab2(:,3);
[pn,ind]   = sort(p./n);
wgti       = w(ind);
nowni      = n(ind);
payi       = p(ind);
pop        = [0;cumsum(wgti.*nowni)/sum(wgti.*nowni)];
pay        = [0;cumsum(wgti.*payi)/sum(wgti.*payi)];
i          = find(pop>=.8);
pay80eb    = pay(i(1))*100;


%=======================================================================
%
% Survey of Business Owners 2007
%   Hours by age group, all owners
%
%=======================================================================

for j=1:6;

%
%  -- 1-owner businesses with owners reporting age and hours
%
i           = find(numown==1 & age1==j);
k           = find(numown==1 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==1 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs1o(j)   = sum1/sum(wgt(k));
nbrep1o(j)  = sum(wgt(i))/1000;

%
%  -- 2-owner businesses with owners reporting age and hours
%
i           = find(numown==2 & age1==j);
k           = find(numown==2 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==2 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs2o1(j)  = sum1/sum(wgt(k));
nbrep2o1(j) = sum(wgt(i))/1000;

i           = find(numown==2 & age2==j);
k           = find(numown==2 & age2==j & hours2>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==2 & age2==j & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs2o2(j)  = sum1/sum(wgt(k));
nbrep2o2(j) = sum(wgt(i))/1000;

%
%  -- 3-owner businesses with owners reporting age and hours
%
i           = find(numown==3 & age1==j);
k           = find(numown==3 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==3 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs3o1(j)  = sum1/sum(wgt(k));
nbrep3o1(j) = sum(wgt(i))/1000;

i           = find(numown==3 & age2==j);
k           = find(numown==3 & age2==j & hours2>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==3 & age2==j & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs3o2(j)  = sum1/sum(wgt(k));
nbrep3o2(j) = sum(wgt(i))/1000;

i           = find(numown==3 & age3==j);
k           = find(numown==3 & age3==j & hours3>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==3 & age3==j & hours3==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs3o3(j)  = sum1/sum(wgt(k));
nbrep3o3(j) = sum(wgt(i))/1000;

%
%  -- 4-owner businesses with owners reporting age and hours
%
i           = find(numown==4 & age1==j);
k           = find(numown==4 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o1(j)  = sum1/sum(wgt(k));
nbrep4o1(j) = sum(wgt(i))/1000;

i           = find(numown==4 & age2==j);
k           = find(numown==4 & age2==j & hours2>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & age2==j & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o2(j)  = sum1/sum(wgt(k));
nbrep4o2(j) = sum(wgt(i))/1000;

i           = find(numown==4 & age3==j);
k           = find(numown==4 & age3==j & hours3>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & age3==j & hours3==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o3(j)  = sum1/sum(wgt(k));
nbrep4o3(j) = sum(wgt(i))/1000;

i           = find(numown==4 & age4==j);
k           = find(numown==4 & age4==j & hours4>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & age4==j & hours4==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o4(j)  = sum1/sum(wgt(k));
nbrep4o4(j) = sum(wgt(i))/1000;
end;

nobyage     = nbrep1o+nbrep2o1+nbrep2o2+nbrep3o1+nbrep3o2+nbrep3o3+ ...
              nbrep4o1+nbrep4o2+nbrep4o3+nbrep4o4;
ahrsbyage   = (ahrs1o.*nbrep1o+ahrs2o1.*nbrep2o1+ahrs2o2.*nbrep2o2+ ...
               ahrs3o1.*nbrep3o1+ahrs3o2.*nbrep3o2+ahrs3o3.*nbrep3o3+ ...
               ahrs4o1.*nbrep4o1+ahrs4o2.*nbrep4o2+ahrs4o3.*nbrep4o3+ ...
               ahrs4o4.*nbrep4o4)./nobyage;
nobyage     = nobyage(:);
ahrsbyage   = ahrsbyage(:);


%=======================================================================
%
% Survey of Business Owners 2007
%   Hours by age group, owners with business as primary income
%
%=======================================================================

for j=1:6;

%
%  -- 1-owner businesses with owners reporting age and hours
%
i           = find(numown==1 & prminc1==1 & age1==j);
k           = find(numown==1 & prminc1==1 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==1 & prminc1==1 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs1o(j)   = sum1/sum(wgt(k));
nbrep1o(j)  = sum(wgt(i))/1000;

%
%  -- 2-owner businesses with owners reporting age and hours
%
i           = find(numown==2 & prminc1==1 & age1==j);
k           = find(numown==2 & prminc1==1 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==2 & prminc1==1 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs2o1(j)  = sum1/sum(wgt(k));
nbrep2o1(j) = sum(wgt(i))/1000;

i           = find(numown==2 & prminc2==1 & age2==j);
k           = find(numown==2 & prminc2==1 & age2==j & hours2>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==2 & prminc2==1 & age2==j & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs2o2(j)  = sum1/sum(wgt(k));
nbrep2o2(j) = sum(wgt(i))/1000;

%
%  -- 3-owner businesses with owners reporting age and hours
%
i           = find(numown==3 & prminc1==1 & age1==j);
k           = find(numown==3 & prminc1==1 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==3 & prminc1==1 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs3o1(j)  = sum1/sum(wgt(k));
nbrep3o1(j) = sum(wgt(i))/1000;

i           = find(numown==3 & prminc2==1 & age2==j);
k           = find(numown==3 & prminc2==1 & age2==j & hours2>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==3 & prminc2==1 & age2==j & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs3o2(j)  = sum1/sum(wgt(k));
nbrep3o2(j) = sum(wgt(i))/1000;

i           = find(numown==3 & prminc3==1 & age3==j);
k           = find(numown==3 & prminc3==1 & age3==j & hours3>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==3 & prminc3==1 & age3==j & hours3==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs3o3(j)  = sum1/sum(wgt(k));
nbrep3o3(j) = sum(wgt(i))/1000;

%
%  -- 4-owner businesses with owners reporting age and hours
%
i           = find(numown==4 & prminc1==1 & age1==j);
k           = find(numown==4 & prminc1==1 & age1==j & hours1>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & prminc1==1 & age1==j & hours1==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o1(j)  = sum1/sum(wgt(k));
nbrep4o1(j) = sum(wgt(i))/1000;

i           = find(numown==4 & prminc2==1 & age2==j);
k           = find(numown==4 & prminc2==1 & age2==j & hours2>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & prminc2==1 & age2==j & hours2==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o2(j)  = sum1/sum(wgt(k));
nbrep4o2(j) = sum(wgt(i))/1000;

i           = find(numown==4 & prminc3==1 & age3==j);
k           = find(numown==4 & prminc3==1 & age3==j & hours3>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & prminc3==1 & age3==j & hours3==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o3(j)  = sum1/sum(wgt(k));
nbrep4o3(j) = sum(wgt(i))/1000;

i           = find(numown==4 & prminc4==1 & age4==j);
k           = find(numown==4 & prminc4==1 & age4==j & hours4>0); 
sum1        = 0;
for l=2:6;
  m         = find(numown==4 & prminc4==1 & age4==j & hours4==l);
  sum1      = sum1+hvec(l)*sum(wgt(m));
end;
ahrs4o4(j)  = sum1/sum(wgt(k));
nbrep4o4(j) = sum(wgt(i))/1000;
end;

nobyagemi   = nbrep1o+nbrep2o1+nbrep2o2+nbrep3o1+nbrep3o2+nbrep3o3+ ...
              nbrep4o1+nbrep4o2+nbrep4o3+nbrep4o4;
ahrsbyagemi = (ahrs1o.*nbrep1o+ahrs2o1.*nbrep2o1+ahrs2o2.*nbrep2o2+ ...
               ahrs3o1.*nbrep3o1+ahrs3o2.*nbrep3o2+ahrs3o3.*nbrep3o3+ ...
               ahrs4o1.*nbrep4o1+ahrs4o2.*nbrep4o2+ahrs4o3.*nbrep4o3+ ...
               ahrs4o4.*nbrep4o4)./nobyagemi;
nobyagemi   = nobyagemi(:);
ahrsbyagemi = ahrsbyagemi(:);

%=======================================================================
%
% Survey of Business Owners 2007
%   Fraction of Exiting Owners that Sell 
%
%=======================================================================

i          = find(operate==2 & numown==1);
j          = find(operate==2 & numown==2);
k          = find(operate==2 & numown==3);
l          = find(operate==2 & numown==4);
m          = find(operate==2 & numown==5);
n          = find(operate==2 & numown==6);
o          = find(operate==2 & numown==7);
sum1       = sum(wgt(i))+sum(wgt(j))*2+sum(wgt(k))*3+sum(wgt(l))*4+ ...
             sum(wgt(m))*7+sum(wgt(n))*30+sum(wgt(o))*75;
i          = find(operate==2 & numown==1 & sold==1);
j          = find(operate==2 & numown==2 & sold==1);
k          = find(operate==2 & numown==3 & sold==1);
l          = find(operate==2 & numown==4 & sold==1);
m          = find(operate==2 & numown==5 & sold==1);
n          = find(operate==2 & numown==6 & sold==1);
o          = find(operate==2 & numown==7 & sold==1);
sum2       = sum(wgt(i))+sum(wgt(j))*2+sum(wgt(k))*3+sum(wgt(l))*4+ ...
             sum(wgt(m))*7+sum(wgt(n))*30+sum(wgt(o))*75;
i          = find(operate==2 & numown==1 & retired==1);
j          = find(operate==2 & numown==2 & retired==1);
k          = find(operate==2 & numown==3 & retired==1);
l          = find(operate==2 & numown==4 & retired==1);
m          = find(operate==2 & numown==5 & retired==1);
n          = find(operate==2 & numown==6 & retired==1);
o          = find(operate==2 & numown==7 & retired==1);
sum3       = sum(wgt(i))+sum(wgt(j))*2+sum(wgt(k))*3+sum(wgt(l))*4+ ...
             sum(wgt(m))*7+sum(wgt(n))*30+sum(wgt(o))*75;
i          = find(operate==2 & numown==1 & dead==1);
j          = find(operate==2 & numown==2 & dead==1);
k          = find(operate==2 & numown==3 & dead==1);
l          = find(operate==2 & numown==4 & dead==1);
m          = find(operate==2 & numown==5 & dead==1);
n          = find(operate==2 & numown==6 & dead==1);
o          = find(operate==2 & numown==7 & dead==1);
sum4       = sum(wgt(i))+sum(wgt(j))*2+sum(wgt(k))*3+sum(wgt(l))*4+ ...
             sum(wgt(m))*7+sum(wgt(n))*30+sum(wgt(o))*75;
fracsold1  = sum2/(sum1-sum3-sum4)*100;

i          = find(operate==2 & numown==1 & prminc1==1);
j1         = find(operate==2 & numown==2 & prminc1==1);
j2         = find(operate==2 & numown==2 & prminc2==1);
k1         = find(operate==2 & numown==3 & prminc1==1);
k2         = find(operate==2 & numown==3 & prminc2==1);
k3         = find(operate==2 & numown==3 & prminc3==1);
l1         = find(operate==2 & numown==4 & prminc1==1);
l2         = find(operate==2 & numown==4 & prminc2==1);
l3         = find(operate==2 & numown==4 & prminc3==1);
l4         = find(operate==2 & numown==4 & prminc4==1);
sum1       = sum(wgt(i))+sum(wgt(j1))+sum(wgt(j2))+sum(wgt(k1))+ ...
             sum(wgt(k2))+sum(wgt(k3))+sum(wgt(l1))+sum(wgt(l2))+ ...
             sum(wgt(l3))+sum(wgt(l4));
i          = find(operate==2 & numown==1 & prminc1==1 & sold==1);
j1         = find(operate==2 & numown==2 & prminc1==1 & sold==1);
j2         = find(operate==2 & numown==2 & prminc2==1 & sold==1);
k1         = find(operate==2 & numown==3 & prminc1==1 & sold==1);
k2         = find(operate==2 & numown==3 & prminc2==1 & sold==1);
k3         = find(operate==2 & numown==3 & prminc3==1 & sold==1);
l1         = find(operate==2 & numown==4 & prminc1==1 & sold==1);
l2         = find(operate==2 & numown==4 & prminc2==1 & sold==1);
l3         = find(operate==2 & numown==4 & prminc3==1 & sold==1);
l4         = find(operate==2 & numown==4 & prminc4==1 & sold==1);
sum2       = sum(wgt(i))+sum(wgt(j1))+sum(wgt(j2))+sum(wgt(k1))+ ...
             sum(wgt(k2))+sum(wgt(k3))+sum(wgt(l1))+sum(wgt(l2))+ ...
             sum(wgt(l3))+sum(wgt(l4));
i          = find(operate==2 & numown==1 & prminc1==1 & retired==1);
j1         = find(operate==2 & numown==2 & prminc1==1 & retired==1);
j2         = find(operate==2 & numown==2 & prminc2==1 & retired==1);
k1         = find(operate==2 & numown==3 & prminc1==1 & retired==1);
k2         = find(operate==2 & numown==3 & prminc2==1 & retired==1);
k3         = find(operate==2 & numown==3 & prminc3==1 & retired==1);
l1         = find(operate==2 & numown==4 & prminc1==1 & retired==1);
l2         = find(operate==2 & numown==4 & prminc2==1 & retired==1);
l3         = find(operate==2 & numown==4 & prminc3==1 & retired==1);
l4         = find(operate==2 & numown==4 & prminc4==1 & retired==1);
sum3       = sum(wgt(i))+sum(wgt(j1))+sum(wgt(j2))+sum(wgt(k1))+ ...
             sum(wgt(k2))+sum(wgt(k3))+sum(wgt(l1))+sum(wgt(l2))+ ...
             sum(wgt(l3))+sum(wgt(l4));
i          = find(operate==2 & numown==1 & prminc1==1 & dead==1);
j1         = find(operate==2 & numown==2 & prminc1==1 & dead==1);
j2         = find(operate==2 & numown==2 & prminc2==1 & dead==1);
k1         = find(operate==2 & numown==3 & prminc1==1 & dead==1);
k2         = find(operate==2 & numown==3 & prminc2==1 & dead==1);
k3         = find(operate==2 & numown==3 & prminc3==1 & dead==1);
l1         = find(operate==2 & numown==4 & prminc1==1 & dead==1);
l2         = find(operate==2 & numown==4 & prminc2==1 & dead==1);
l3         = find(operate==2 & numown==4 & prminc3==1 & dead==1);
l4         = find(operate==2 & numown==4 & prminc4==1 & dead==1);
sum4       = sum(wgt(i))+sum(wgt(j1))+sum(wgt(j2))+sum(wgt(k1))+ ...
             sum(wgt(k2))+sum(wgt(k3))+sum(wgt(l1))+sum(wgt(l2))+ ...
             sum(wgt(l3))+sum(wgt(l4));
fracsold2  = sum2/(sum1-sum3-sum4)*100;


%=======================================================================
%
% Census Business Patterns (Annual)
%
%=======================================================================

ccnest      = NaN*ones(nyr,10);
ccnemp      = NaN*ones(nyr,10);
ccpayr      = NaN*ones(nyr,10);
ptnest      = NaN*ones(nyr,10);
ptnemp      = NaN*ones(nyr,10);
ptpayr      = NaN*ones(nyr,10);
for t=t0:t1;
  if ((t>=2009) & (t<=2013));
    tab         = xlsread(['Census/BP/BP',int2str(t)]);
    i           = t-t0+1;
    tab         = tab/10^6;
    ccnest(i,:) = tab(11:20,9)';
    ccnemp(i,:) = tab(11:20,10)';
    ccpayr(i,:) = tab(11:20,11)';
    ptnest(i,:) = tab(21:30,9)'+tab(31:40,9)'+tab(41:50,9)';
    ptnemp(i,:) = tab(21:30,10)'+tab(31:40,10)'+tab(41:50,10)';
    ptpayr(i,:) = tab(21:30,11)'+tab(31:40,11)'+tab(41:50,11)';
  end;
end;


%=======================================================================
%
% Adjusted data:
%
%=======================================================================

xintang      = 0.12*gdp;
kintang      = 1.7*gdp;
dintang      = 0.05*kintang;

impkcds      = .04*kcd;
impkgs       = .04*kg;
trusccomp    = sccomp/(1-.18);
truscneti    = scneti/(1-.18);
misrepsc     = trusccomp+truscneti-sccomp-scneti;

adjgdi       = gdi+desc-stax+cddep+impkcds+impkgs-gpdinripp+xintang+misrepsc;
adjcomp      = comp-sccomp-dhhcomp-dnpcomp-dopbuscomp-dgebuscomp-dggcomp;
sclabc       = labc.*sccogs./cogs;
csw2inc      = labc+salw+pena+empb+coff-sccomp;
scw2inc      = sclabc+scsalw+scpena+scempb;
adjcccomp    = (dcbuscomp-sccomp).*(1-scw2inc./csw2inc); 
adjptcomp    = dppbuscomp+(dcbuscomp-sccomp).*scw2inc./csw2inc; 
adjcapi      = (prof+piadj)+(rent-dhhrent-dnprent-dopbusrent)+ ...
               (nint-dhhnint-dnpnint-dopbusnint)+(ibt-dhhibt-dnpibt- ...
               dopbusibt-stax)-gpdinripp+depipp+xintang-dintang-scneti;
sweat        = pinc-dopbuspinc-piadj+trusccomp+truscneti+desc;
adjdep       = dep+cddep-depres-depgov-depipp+dintang;
nonbus       = dhh+dnp+dopbus+dgebus+dgg+depres+depgov+impkcds+impkgs;

adjgdp       = gdp-stax+cddep+impkcds+impkgs-gpdinripp+xintang+misrepsc;
dstax        = pced.*stax./pce;
ndstax       = stax-dstax;
adjcons      = pce-pced-ndstax+cddep+impkcds+impkgs+misrepsc;

adjinv       = gpdi-gpdinripp+xintang +pced-dstax+govi+netxm;
xcands       = xpnrcseq+xpnrcsstr+xprescs+gpdiinv;
adjxcc       = (1-scasset./asset).*(xcands+xpnrcsipp./xpnripp.*xintang);
adjxsc       = (scasset./asset).*(xcands+xpnrcsipp./xpnripp.*xintang);
adjxsp       = xpnrspeq+xpnrspstr+xpressp+xpnrspipp./xpnripp.*xintang;
adjxpa       = xpnrpaeq+xpnrpastr+xpnrpaipp./xpnripp.*xintang;
adjxhh       = xpnrhheq+xpnrhhstr+xpreshh+xpnrhhipp./xpnripp.*xintang;
adjxnp       = xpnrnpeq+xpnrnpstr+xpnrteeq+xpnrtestr+xpresnp+xpnrnpipp./ ...
               xpnripp.*xintang+xpnrteipp./xpnripp.*xintang;
scale        = (gpdi-gpdinripp+xintang)./(gpdiinv+xpnreq+xpnrstr+xpres+ ...
               xintang);
adjxcc       = scale.*adjxcc;
adjxsc       = scale.*adjxsc;
adjxpa       = scale.*adjxpa;
adjxsp       = scale.*adjxsp;
adjxhh       = scale.*adjxhh;
adjxnp       = scale.*adjxnp;
adjxnb       = pced-dstax+govi+adjxhh+adjxnp+netxm;
adjxpt       = adjxsc+adjxsp+adjxpa;

kcands       = kpnrcs+kprescs;
ksp          = kpnrsp+kpressp;
kpa          = kpnrpa;
khh          = kpnrhh+kpreshh+kcd;
knp          = kpnrnp+kpnrte+kpresnp;
knb          = kcd+kg+kpnrhh+kpnrnp+kpnrte+kpreshh+kpresnp;
land         = cland+ncland+hhland;
ktot         = kpfix+kg+kcd+kinv+land;

adjkcands    = kcands-kpnrcsipp+kcinv;
adjkcc       = (1-scasset./asset).*(adjkcands+kpnrcsipp./kpnripp.*kintang);
adjksc       = (scasset./asset).*(adjkcands+kpnrcsipp./kpnripp.*kintang);
adjksp       = ksp-kpnrspipp+kncinv+kpnrspipp./kpnripp.*kintang;
adjkpa       = kpa-kpnrpaipp+kpnrpaipp./kpnripp.*kintang;
adjknb       = knb-kpnrhhipp-kpnrnpipp-kpnrteipp+(kpnrhhipp+kpnrnpipp+ ...
               kpnrteipp)./kpnripp.*kintang;
adjkpt       = adjksc+adjksp+adjkpa;
adjk         = adjkcc+adjkpt+adjknb;

%
% Remove capital income from sweat income (all interest towards FA)
%
ptcapownshr  = 1-(ptintp/.042+ptrent/.093)./adjkpt;
ptcapinc     = parerir+pareotr+ptcapownshr.*ptdepr*.042/.051;
adjsweat     = sweat-ptcapinc;

ptcapownshr2 = 1-(0*ptintp/.042+ptrent/.093)./adjkpt; %No interest to FA
ptcapinc2    = parerir+pareotr+ptcapownshr2.*ptdepr*.042/.051;
adjsweat2    = sweat-ptcapinc2;
ptcapinc3    = .042*adjkpt-ptintp-.042/.093*ptrent; %Alternative imputation
adjsweat3    = sweat-ptcapinc3;

ptcapinc1    = .5*(ptcapinc+ptcapinc2);
adjsweat1    = sweat-ptcapinc1;


t            = [t0:t1]';
y            = adjgdp;
%=======================================================================
%
% Print relevant tables
%
%=======================================================================

if (iprint==1);
clc
fprintf('NATIONAL ACCOUNTS AND FIXED ASSETS\n')
fprintf('==================================\n')
fprintf('  * NIPA without adjustments or normalizations\n')
fprintf('  * NIPA with adjustments but no normalizations\n')
fprintf('  * FA without adjustments or normalizations\n')
fprintf('  * FA with adjustments but no normalizations\n')
fprintf('  * NIPA with adjustments and normalizations \n')
fprintf('  * FA with adjustments and normalizations \n')
fprintf('  * Other moments needed for parameterization\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause
end;

clc
fprintf('IMPUTATIONS\n')
fprintf('===========\n')
fprintf('  X_intang/GDP                     %7.2f\n',mean(xintang./gdp))
fprintf('  K_intang/GDP                     %7.2f\n',mean(kintang./gdp))
fprintf('  Dep_intang/K_intang              %7.2f\n',mean(dintang./kintang))
fprintf('  CD services/K_cd                 %7.2f\n',mean(impkcds./kcd))
fprintf('  Govt inv services/K_g            %7.2f\n',mean(impkgs./kg))
fprintf('  Noncorp share of inventories     %7.2f\n',mean(kncinv./kinv))
fprintf('  Reported S corp comp/True        %7.2f\n',mean(sccomp./trusccomp))
fprintf('  Reported S corp net inc/True     %7.2f\n',mean(scneti./truscneti))
fprintf('  S-corp share of C&S capital      %7.2f\n',mean(scasset./asset))
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause

if (iprint==1);
clc
fprintf('                  TABLE 1\n')
fprintf('UNADJUSTED DOMESTIC ACCOUNTS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('DOMESTIC INCOME                    %7.0f\n',mean(gdi))
fprintf('  Compensation                     %7.0f\n',mean(comp))
fprintf('    Corporate                      %7.0f\n',mean(dcbuscomp))
fprintf('    Noncorporate                   %7.0f\n',mean(dncbuscomp))
fprintf('      Proprietors & Partners       %7.0f\n',mean(dppbuscomp))
fprintf('      Other private                %7.0f\n',mean(dopbuscomp))
fprintf('      Govt Enterprises             %7.0f\n',mean(dgebuscomp))
fprintf('    Households                     %7.0f\n',mean(dhhcomp))
fprintf('    Nonprofit                      %7.0f\n',mean(dnpcomp))
fprintf('    Government                     %7.0f\n',mean(dggcomp))
fprintf('  Corporate profit                 %7.0f\n',mean(prof))
fprintf('  Proprietors'' income              %7.0f\n',mean(pinc))
fprintf('  Rental income                    %7.0f\n',mean(rent))
fprintf('  Net interest                     %7.0f\n',mean(nint))
fprintf('  Government enterprises           %7.0f\n',mean(gprof))
fprintf('  Consumption of fixed capital     %7.0f\n',mean(dep))
fprintf('  Indirect business taxes          %7.0f\n',mean(ibt))

fprintf(' \n')
fprintf('DOMESTIC PRODUCT                   %7.0f\n',mean(gdp))
fprintf('  Personal consumption expend      %7.0f\n',mean(pce))
fprintf('    Durable                        %7.0f\n',mean(pced))
fprintf('    Nondurable                     %7.0f\n',mean(pcend))
fprintf('    Services                       %7.0f\n',mean(pces))
fprintf('  Grosss private domestic inv      %7.0f\n',mean(gpdi))
fprintf('    Nonres structures              %7.0f\n',mean(gpdinrstr))
fprintf('    Nonres equipment               %7.0f\n',mean(gpdinreq))
fprintf('    Nonres IPP                     %7.0f\n',mean(gpdinripp))
fprintf('    Residential                    %7.0f\n',mean(gpdires))
fprintf('    Change in inventories          %7.0f\n',mean(gpdiinv))
fprintf('  Government C&I                   %7.0f\n',mean(gov))
fprintf('    Consumption                    %7.0f\n',mean(govc))
fprintf('    Investment                     %7.0f\n',mean(govi))
fprintf('  Net exports                      %7.0f\n',mean(netxm))
fprintf('    Exports                        %7.0f\n',mean(netx))
fprintf('    Imports                        %7.0f\n',mean(netm))
fprintf('============================================\n')
fprintf('\n')
fprintf('Hit return to continue...\n')
fprintf('\n')
pause

clc

fprintf('                  TABLE 2\n')
fprintf('ADJUSTED DOMESTIC ACCOUNTS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('ADJUSTED GDI                       %7.0f\n',mean(adjgdi))
fprintf('  NIPA GDI                         %7.0f\n',mean(gdi))
fprintf('    +NIPA statistical desc.        %7.0f\n',mean(desc))
fprintf('    +FA durable depreciation       %7.0f\n',mean(cddep))
fprintf('    +Misreported S corp inc        %7.0f\n',mean(misrepsc))
fprintf('    +Imputed Kcd services          %7.0f\n',mean(impkcds))
fprintf('    +Imputed Kg services           %7.0f\n',mean(impkgs))
fprintf('    +Imputed private IPP inv       %7.0f\n',mean(xintang))
fprintf('    -NIPA private IPP inv          %7.0f\n',mean(gpdinripp))
fprintf('    -NIPA sales tax                %7.0f\n',mean(stax))
fprintf('\n')
fprintf('= Sweat income                     %7.0f\n',mean(adjsweat1))
fprintf('   NIPA proprietors'' income        %7.0f\n',mean(pinc))
fprintf('    -NIPA proprietors'' IVA,CCadj   %7.0f\n',mean(piadj))
fprintf('    -NIPA prop inc, other private  %7.0f\n',mean(dopbuspinc))
fprintf('    +IRS S corp comp (true)        %7.0f\n',mean(trusccomp))
fprintf('    +IRS S corp bus inc (true)     %7.0f\n',mean(truscneti))
fprintf('    +NIPA statistical desc.        %7.0f\n',mean(desc))
fprintf('    -Estimate of P-T cap inc       %7.0f\n',mean(ptcapinc1))
fprintf('\n')
fprintf('+ Nonsweat labor income            %7.0f\n',mean(adjcomp))
fprintf('   NIPA compensation               %7.0f\n',mean(comp))
fprintf('    -IRS S corp comp (reported)    %7.0f\n',mean(sccomp))
fprintf('    -NIPA household comp           %7.0f\n',mean(dhhcomp))
fprintf('    -NIPA nonprof comp             %7.0f\n',mean(dnpcomp))
fprintf('    -NIPA other private comp       %7.0f\n',mean(dopbuscomp))
fprintf('    -NIPA govt ent comp            %7.0f\n',mean(dgebuscomp))
fprintf('    -NIPA gen govt comp            %7.0f\n',mean(dggcomp))
fprintf('    ==> C-corp share               %7.0f\n',mean(adjcccomp))
fprintf('        Pass-thru share            %7.0f\n',mean(adjptcomp))
fprintf('\n')
fprintf('+ Capital income                   %7.0f\n',mean(adjcapi+ptcapinc1))
fprintf('   NIPA corporate profits          %7.0f\n',mean(prof))
fprintf('    +NIPA prop inc iva,ccadj       %7.0f\n',mean(piadj))
fprintf('    +NIPA rental income            %7.0f\n',mean(rent))
fprintf('    -NIPA household rental         %7.0f\n',mean(dhhrent))
fprintf('    -NIPA nonprof rental           %7.0f\n',mean(dnprent))
fprintf('    -NIPA other bus rental         %7.0f\n',mean(dopbusrent))
fprintf('    +NIPA net interest             %7.0f\n',mean(nint))
fprintf('    -NIPA household net int        %7.0f\n',mean(dhhnint))
fprintf('    -NIPA nonprof net int          %7.0f\n',mean(dnpnint))
fprintf('    -NIPA other bus net int        %7.0f\n',mean(dopbusnint))
fprintf('    +NIPA IBT                      %7.0f\n',mean(ibt))
fprintf('    -NIPA household IBT            %7.0f\n',mean(dhhibt))
fprintf('    -NIPA nonprofit IBT            %7.0f\n',mean(dnpibt))
fprintf('    -NIPA other bus IBT            %7.0f\n',mean(dopbusibt))
fprintf('    -NIPA sales tax                %7.0f\n',mean(stax))
fprintf('    -NIPA private IPP inv          %7.0f\n',mean(gpdinripp))
fprintf('    +Imputed private IPP inv       %7.0f\n',mean(xintang))
fprintf('    +NIPA IPP depreciation         %7.0f\n',mean(depipp))
fprintf('    -Imputed private IPP dep       %7.0f\n',mean(dintang))
fprintf('    -IRS S corp bus inc (reported) %7.0f\n',mean(scneti))
fprintf('    +Estimate of P-T cap inc       %7.0f\n',mean(ptcapinc1))
fprintf('\n')
fprintf('+ Capital depreciation             %7.0f\n',mean(adjdep))
fprintf('   NIPA cons of fixed capital      %7.0f\n',mean(dep))
fprintf('    +FA durable depreciation       %7.0f\n',mean(cddep))
fprintf('    -NIPA residential dep          %7.0f\n',mean(depres))
fprintf('    -NIPA government dep           %7.0f\n',mean(depgov))
fprintf('    -NIPA IPP depreciation         %7.0f\n',mean(depipp))
fprintf('    +Imputed private IPP dep       %7.0f\n',mean(dintang))
fprintf('\n')
fprintf('+ Nonbusiness income               %7.0f\n',mean(nonbus))
fprintf('   NIPA household income           %7.0f\n',mean(dhh))
fprintf('    +NIPA nonprofit income         %7.0f\n',mean(dnp))
fprintf('    +NIPA other bus income         %7.0f\n',mean(dopbus))
fprintf('    +NIPA govt ent income          %7.0f\n',mean(dgebus))
fprintf('    +NIPA gen govt income          %7.0f\n',mean(dgg))
fprintf('    +NIPA residential dep          %7.0f\n',mean(depres))
fprintf('    +NIPA government dep           %7.0f\n',mean(depgov))
fprintf('    +Imputed Kcd services          %7.0f\n',mean(impkcds))
fprintf('    +Imputed Kg services           %7.0f\n',mean(impkgs))
fprintf(' \n')
fprintf(' \n')
fprintf('ADJUSTED GDP                       %7.0f\n',mean(adjgdp))
fprintf('  NIPA GDP                         %7.0f\n',mean(gdp))
fprintf('    +FA durable depreciation       %7.0f\n',mean(cddep))
fprintf('    +Misreported S corp inc        %7.0f\n',mean(misrepsc))
fprintf('    +Imputed Kcd services          %7.0f\n',mean(impkcds))
fprintf('    +Imputed Kg services           %7.0f\n',mean(impkgs))
fprintf('    +Imputed private IPP inv       %7.0f\n',mean(xintang))
fprintf('    -NIPA private IPP inv          %7.0f\n',mean(gpdinripp))
fprintf('    -NIPA sales tax                %7.0f\n',mean(stax))
fprintf('= Private consumption              %7.0f\n',mean(adjcons))
fprintf('   NIPA PCE                        %7.0f\n',mean(pce))
fprintf('    -NIPA PCE durables             %7.0f\n',mean(pced))
fprintf('    +FA durable depreciation       %7.0f\n',mean(cddep))
fprintf('    +Misreported S corp inc        %7.0f\n',mean(misrepsc))
fprintf('    +Imputed Kcd services          %7.0f\n',mean(impkcds))
fprintf('    +Imputed Kg services           %7.0f\n',mean(impkgs))
fprintf('    -NIPA pro-rata sales tax       %7.0f\n',mean(ndstax))
fprintf('+ Public consumption               %7.0f\n',mean(govc))
fprintf('+ Investments                      %7.0f\n',mean(adjinv))
fprintf('   NIPA GDPI                       %7.0f\n',mean(gpdi))
fprintf('    -NIPA GDPI nonres IPP          %7.0f\n',mean(gpdinripp))
fprintf('    +Imputed private IPP inv       %7.0f\n',mean(xintang))
fprintf('    +NIPA PCE durables             %7.0f\n',mean(pced))
fprintf('    -NIPA pro-rata sales tax       %7.0f\n',mean(dstax))
fprintf('    +NIPA government investment    %7.0f\n',mean(govi))
fprintf('    +NIPA net exports              %7.0f\n',mean(netxm))
fprintf('    ==> C-corp share               %7.0f\n',mean(adjxcc))
fprintf('        Pass-thru share            %7.0f\n',mean(adjxpt))
fprintf('          S corporations           %7.0f\n',mean(adjxsc))
fprintf('          Sole proprietorships     %7.0f\n',mean(adjxsp))
fprintf('          Partnerships             %7.0f\n',mean(adjxpa))
fprintf('        Nonbusiness (incl NX)      %7.0f\n',mean(adjxnb))
fprintf('============================================\n')
fprintf('\n')
fprintf('Hit return to continue...\n')
fprintf('\n')
pause

clc
fprintf('                  TABLE 3\n')
fprintf('DETAILS OF P-T CAPITAL INCOME, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('ASSUMPITONS:\n')
fprintf('  %% Rate of return (r)             %7.1f\n',4.2)
fprintf('  %% Rate of depreciation (del)     %7.1f\n',5.1)
fprintf('IRS DATA: \n')
fprintf('  Real estate inc (rei)            %7.0f\n',mean(parerir+pareotr))
fprintf('  Depreciation (dep)               %7.0f\n',mean(ptdepr))
fprintf('  Interest paid (int)              %7.0f\n',mean(ptintp))
fprintf('  Rent paid (rnt)                  %7.0f\n',mean(ptrent))
fprintf('  Capital stock (k)                %7.0f\n',mean(adjkpt))
fprintf('  %% Ownership share (oshr1):           \n')
fprintf('    1-(int/r+rnt/(r+del))/k        %7.1f\n',mean(ptcapownshr*100))
fprintf('  %% Ownership share (oshr2):           \n')
fprintf('    1-(rnt/(r+del))/k              %7.1f\n',mean(ptcapownshr2*100))
fprintf('\n')
fprintf('PASS-THROUGH CAPITAL INCOME ESTIMATES: \n')
fprintf('  rei+r*oshr1*dep/del              %7.0f\n',mean(ptcapinc))
fprintf('  rei+r*oshr2*dep/del              %7.0f\n',mean(ptcapinc2))
fprintf('  r*k-int-r/(r+del)*rnt            %7.0f\n',mean(ptcapinc3))
fprintf('\n')
fprintf('ESTIMATES RELATIVE TO ADJ GDP:\n')
fprintf('  rei+r*oshr1*dep/del              %7.3f\n',mean(ptcapinc./y))
fprintf('  rei+r*oshr2*dep/del              %7.3f\n',mean(ptcapinc2./y))
fprintf('  r*k-int-r/(r+del)*rnt            %7.3f\n',mean(ptcapinc3./y))
fprintf('============================================\n')
fprintf('\n')
fprintf('Hit return to continue...\n')
fprintf('\n')

clc
fprintf('                  TABLE 4\n')
fprintf('UNADJUSTED CAPITAL STOCKS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('DOMESTIC CAPITAL                   %7.0f\n',mean(ktot))
fprintf('  FIXED ASSETS & CD                %7.0f\n',mean(kfix+kcd))
fprintf('    Corporate                      %7.0f\n',mean(kcands))
fprintf('      Structures                   %7.0f\n',mean(kpnrcsstr))
fprintf('      Equipment                    %7.0f\n',mean(kpnrcseq))
fprintf('      IPP                          %7.0f\n',mean(kpnrcsipp))
fprintf('      Residential                  %7.0f\n',mean(kprescs))
fprintf('    Sole proprietors               %7.0f\n',mean(ksp))
fprintf('      Structures                   %7.0f\n',mean(kpnrspstr))
fprintf('      Equipment                    %7.0f\n',mean(kpnrspeq))
fprintf('      IPP                          %7.0f\n',mean(kpnrspipp))
fprintf('      Residential                  %7.0f\n',mean(kpressp))
fprintf('    Partnerships                   %7.0f\n',mean(kpa))
fprintf('      Structures                   %7.0f\n',mean(kpnrpastr))
fprintf('      Equipment                    %7.0f\n',mean(kpnrpaeq))
fprintf('      IPP                          %7.0f\n',mean(kpnrpaipp))
fprintf('    Nonbusiness                    %7.0f\n',mean(knb))
fprintf('      Households                   %7.0f\n',mean(khh))
fprintf('      Nonprofits                   %7.0f\n',mean(knp))
fprintf('      Government                   %7.0f\n',mean(kg))
fprintf('\n')
fprintf('  INVENTORIES                      %7.0f\n',mean(kinv))
fprintf('    Corporate                      %7.0f\n',mean(kcinv))
fprintf('    Noncorporate                   %7.0f\n',mean(kncinv))
fprintf('\n')
fprintf('  LAND                             %7.0f\n',mean(land))
fprintf('    Corporate                      %7.0f\n',mean(cland))
fprintf('    Noncorporate                   %7.0f\n',mean(ncland))
fprintf('    Household/Nonprofits           %7.0f\n',mean(hhland))
fprintf('============================================\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause

clc
fprintf('                  TABLE 5\n')
fprintf('ADJUSTED CAPITAL STOCKS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('ADJUSTED CAPITAL                   %7.0f\n',mean(adjk+land))
fprintf('  Adjusted FA, CD, and inventories %7.0f\n',mean(adjk))
fprintf('    FA fixed assets & CD           %7.0f\n',mean(ktot-kinv-land))
fprintf('    +NIPA inventories              %7.0f\n',mean(kinv))
fprintf('    -FA IPP                        %7.0f\n',mean(kpnripp))
fprintf('    +Imputed private IPP capital   %7.0f\n',mean(kintang))
fprintf(' +Land                             %7.0f\n',mean(land))
fprintf('\n')
fprintf(' = C-CORP CAPITAL                  %7.0f\n',mean(adjkcc+cland))
fprintf('     C corporations                %7.0f\n',mean(adjkcc))
fprintf('     Corporate land                %7.0f\n',mean(cland))
fprintf(' + PASS-THROUGH CAPITAL            %7.0f\n',mean(adjkpt+ncland))
fprintf('     S corporations                %7.0f\n',mean(adjksc))
fprintf('     Sole proprietors              %7.0f\n',mean(adjksp))
fprintf('     Partnerships                  %7.0f\n',mean(adjkpa))
fprintf('     Noncorporate land             %7.0f\n',mean(ncland))
fprintf(' + NONBUSINESS CAPITAL             %7.0f\n',mean(adjknb+hhland))
fprintf('============================================\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause
end;


clc
if (iprint==0);
fprintf('                  TABLE 1\n')
else
fprintf('                  TABLE 6\n')
end;
fprintf('ADJUSTED DOMESTIC ACCOUNTS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('ADJUSTED GDI                       %7.3f\n',mean(adjgdi./y))
fprintf('  NIPA GDI                         %7.3f\n',mean(gdi./y))
fprintf('    +NIPA statistical desc.        %7.3f\n',mean(desc./y))
fprintf('    +FA durable depreciation       %7.3f\n',mean(cddep./y))
fprintf('    +Misreported S corp inc        %7.3f\n',mean(misrepsc./y))
fprintf('    +Imputed Kcd services          %7.3f\n',mean(impkcds./y))
fprintf('    +Imputed Kg services           %7.3f\n',mean(impkgs./y))
fprintf('    +Imputed private IPP inv       %7.3f\n',mean(xintang./y))
fprintf('    -NIPA private IPP inv          %7.3f\n',mean(gpdinripp./y))
fprintf('    -NIPA sales tax                %7.3f\n',mean(stax./y))
fprintf('\n')
fprintf('= Sweat income                     %7.3f\n',mean(adjsweat1./y))
fprintf('   NIPA proprietors'' income        %7.3f\n',mean(pinc./y))
fprintf('    -NIPA proprietors'' IVA,CCadj   %7.3f\n',mean(piadj./y))
fprintf('    -NIPA prop inc, other private  %7.0f\n',mean(dopbuspinc./y))
fprintf('    +IRS S corp comp (true)        %7.3f\n',mean(trusccomp./y))
fprintf('    +IRS S corp bus inc (true)     %7.3f\n',mean(truscneti./y))
fprintf('    +NIPA statistical desc.        %7.3f\n',mean(desc./y))
fprintf('    -Estimate of P-T cap inc       %7.3f\n',mean(ptcapinc1./y))
fprintf('\n')
fprintf('+ Nonsweat labor income            %7.3f\n',mean(adjcomp./y))
fprintf('   NIPA compensation               %7.3f\n',mean(comp./y))
fprintf('    -IRS S corp comp (reported)    %7.3f\n',mean(sccomp./y))
fprintf('    -NIPA household comp           %7.3f\n',mean(dhhcomp./y))
fprintf('    -NIPA nonprof comp             %7.3f\n',mean(dnpcomp./y))
fprintf('    -NIPA other bus comp           %7.3f\n',mean(dopbuscomp./y))
fprintf('    -NIPA govt ent comp            %7.3f\n',mean(dgebuscomp./y))
fprintf('    -NIPA gen govt comp            %7.3f\n',mean(dggcomp./y))
fprintf('    ==> C-corp share               %7.3f\n',mean(adjcccomp./y))
fprintf('        Pass-thru share            %7.3f\n',mean(adjptcomp./y))
fprintf('\n')
fprintf('+ Nonsweat capital income          %7.3f\n',mean((adjcapi+ptcapinc1)./y))
fprintf('   NIPA corporate profits          %7.3f\n',mean(prof./y))
fprintf('    +NIPA prop inc iva,ccadj       %7.3f\n',mean(piadj./y))
fprintf('    +NIPA rental income            %7.3f\n',mean(rent./y))
fprintf('    -NIPA household rental         %7.3f\n',mean(dhhrent./y))
fprintf('    -NIPA nonprof rental           %7.3f\n',mean(dnprent./y))
fprintf('    -NIPA other bus rental         %7.3f\n',mean(dopbusrent./y))
fprintf('    +NIPA net interest             %7.3f\n',mean(nint./y))
fprintf('    -NIPA household net int        %7.3f\n',mean(dhhnint./y))
fprintf('    -NIPA nonprof net int          %7.3f\n',mean(dnpnint./y))
fprintf('    -NIPA other bus net int        %7.3f\n',mean(dopbusnint./y))
fprintf('    +NIPA IBT                      %7.3f\n',mean(ibt./y))
fprintf('    -NIPA household IBT            %7.3f\n',mean(dhhibt./y))
fprintf('    -NIPA nonprofit IBT            %7.3f\n',mean(dnpibt./y))
fprintf('    -NIPA other bus IBT            %7.3f\n',mean(dopbusibt./y))
fprintf('    -NIPA sales tax                %7.3f\n',mean(stax./y))
fprintf('    -NIPA private IPP inv          %7.3f\n',mean(gpdinripp./y))
fprintf('    +Imputed private IPP inv       %7.3f\n',mean(xintang./y))
fprintf('    +NIPA IPP depreciation         %7.3f\n',mean(depipp./y))
fprintf('    -Imputed private IPP dep       %7.3f\n',mean(dintang./y))
fprintf('    -IRS S corp bus inc (reported) %7.3f\n',mean(scneti./y))
fprintf('    +Estimate of P-T cap inc       %7.3f\n',mean(ptcapinc1./y))
fprintf('\n')
fprintf('+ Capital depreciation             %7.3f\n',mean(adjdep./y))
fprintf('   NIPA cons of fixed capital      %7.3f\n',mean(dep./y))
fprintf('    +FA durable depreciation       %7.3f\n',mean(cddep./y))
fprintf('    -NIPA residential dep          %7.3f\n',mean(depres./y))
fprintf('    -NIPA government dep           %7.3f\n',mean(depgov./y))
fprintf('    -NIPA IPP depreciation         %7.3f\n',mean(depipp./y))
fprintf('    +Imputed private IPP dep       %7.3f\n',mean(dintang./y))
fprintf('\n')
fprintf('+ Nonbusiness income               %7.3f\n',mean(nonbus./y))
fprintf('   NIPA household income           %7.3f\n',mean(dhh./y))
fprintf('    +NIPA nonprofit income         %7.3f\n',mean(dnp./y))
fprintf('    +NIPA other bus income         %7.3f\n',mean(dopbus./y))
fprintf('    +NIPA govt ent income          %7.3f\n',mean(dgebus./y))
fprintf('    +NIPA gen govt income          %7.3f\n',mean(dgg./y))
fprintf('    +NIPA residential dep          %7.3f\n',mean(depres./y))
fprintf('    +NIPA government dep           %7.3f\n',mean(depgov./y))
fprintf('    +Imputed Kcd services          %7.3f\n',mean(impkcds./y))
fprintf('    +Imputed Kg services           %7.3f\n',mean(impkgs./y))
fprintf(' \n')
fprintf(' \n')
fprintf('ADJUSTED GDP                       %7.3f\n',mean(adjgdp./y))
fprintf('  NIPA GDP                         %7.3f\n',mean(gdp./y))
fprintf('    +FA durable depreciation       %7.3f\n',mean(cddep./y))
fprintf('    +Misreported S corp inc        %7.3f\n',mean(misrepsc./y))
fprintf('    +Imputed Kcd services          %7.3f\n',mean(impkcds./y))
fprintf('    +Imputed Kg services           %7.3f\n',mean(impkgs./y))
fprintf('    +Imputed private IPP inv       %7.3f\n',mean(xintang./y))
fprintf('    -NIPA private IPP inv          %7.3f\n',mean(gpdinripp./y))
fprintf('    -NIPA sales tax                %7.3f\n',mean(stax./y))
fprintf('= Private consumption              %7.3f\n',mean(adjcons./y))
fprintf('   NIPA PCE                        %7.3f\n',mean(pce./y))
fprintf('    -NIPA PCE durables             %7.3f\n',mean(pced./y))
fprintf('    +FA durable depreciation       %7.3f\n',mean(cddep./y))
fprintf('    +Misreported S corp inc        %7.3f\n',mean(misrepsc./y))
fprintf('    +Imputed Kcd services          %7.3f\n',mean(impkcds./y))
fprintf('    +Imputed Kg services           %7.3f\n',mean(impkgs./y))
fprintf('    -NIPA pro-rata sales tax       %7.3f\n',mean(ndstax./y))
fprintf('+ Public consumption               %7.3f\n',mean(govc./y))
fprintf('+ Investments                      %7.3f\n',mean(adjinv./y))
fprintf('   NIPA GDPI                       %7.3f\n',mean(gpdi./y))
fprintf('    -NIPA GDPI nonres IPP          %7.3f\n',mean(gpdinripp./y))
fprintf('    +Imputed private IPP inv       %7.3f\n',mean(xintang./y))
fprintf('    +NIPA PCE durables             %7.3f\n',mean(pced./y))
fprintf('    -NIPA pro-rata sales tax       %7.3f\n',mean(dstax./y))
fprintf('    +NIPA government investment    %7.3f\n',mean(govi./y))
fprintf('    +NIPA net exports              %7.3f\n',mean(netxm./y))
fprintf('    ==> C-corp share               %7.3f\n',mean(adjxcc./y))
fprintf('        Pass-thru share            %7.3f\n',mean(adjxpt./y))
fprintf('          S corporations           %7.3f\n',mean(adjxsc./y))
fprintf('          Sole proprietorships     %7.3f\n',mean(adjxsp./y))
fprintf('          Partnerships             %7.3f\n',mean(adjxpa./y))
fprintf('        Nonbusiness (incl NX)      %7.3f\n',mean(adjxnb./y))
fprintf(' \n')
fprintf(' Addendum: \n')
fprintf('    Govt net interest payment      %7.3f\n',mean((intpay-intrec)./y))
fprintf('    Tax rate on consumption        %7.3f\n',mean(stax./pce*100))
fprintf('    Tax rate on C-corp profits     %7.3f\n',mean(txrate))
fprintf('    Tax rate on C-corp profits*    %7.3f\n',mean(txrate2))
fprintf('    Retirement financing           %7.3f\n',mean(retfinance./y))
fprintf('      Social security              %7.3f\n',mean(socsec./y))
fprintf('      Medicare                     %7.3f\n',mean(medicare./y))
fprintf('      Railroad retirement          %7.3f\n',mean(rrret./y))
fprintf('      Pension benefit guaranty     %7.3f\n',mean(pbg./y))
fprintf('      Veteran pension+disability   %7.3f\n',mean(vetpen./y))
fprintf('    Sweat income estimates:\n')
fprintf('      Interest for fixed assets    %7.3f\n',mean(adjsweat./y))
fprintf('      No interest for FA           %7.3f\n',mean(adjsweat2./y))
fprintf('      Alternative imputation       %7.3f\n',mean(adjsweat3./y))
fprintf('============================================\n')
fprintf('* with domestic production deduction\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause

clc
if (iprint==0);
fprintf('                  TABLE 2\n')
else
fprintf('                  TABLE 7\n')
end;
fprintf('ADJUSTED CAPITAL STOCKS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('ADJUSTED CAPITAL                   %7.3f\n',mean((adjk+land)./y))
fprintf('  Adjusted FA, CD, and inventories %7.3f\n',mean(adjk./y))
fprintf('    FA fixed assets & CD           %7.3f\n',mean((ktot-kinv-land)./y))
fprintf('    +NIPA inventories              %7.3f\n',mean(kinv./y))
fprintf('    -FA IPP                        %7.3f\n',mean(kpnripp./y))
fprintf('    +Imputed private IPP capital   %7.3f\n',mean(kintang./y))
fprintf(' +Land                             %7.3f\n',mean(land./y))
fprintf('\n')
fprintf(' = C-CORP CAPITAL                  %7.3f\n',mean((adjkcc+cland)./y))
fprintf('     C corporations                %7.3f\n',mean(adjkcc./y))
fprintf('     Corporate land                %7.3f\n',mean(cland./y))
fprintf(' + PASS-THROUGH CAPITAL            %7.3f\n',mean((adjkpt+ncland)./y))
fprintf('     S corporations                %7.3f\n',mean(adjksc./y))
fprintf('     Sole proprietors              %7.3f\n',mean(adjksp./y))
fprintf('     Partnerships                  %7.3f\n',mean(adjkpa./y))
fprintf('     Non-corporate land            %7.3f\n',mean(ncland./y))
fprintf(' + NONBUSINESS CAPITAL             %7.3f\n',mean((adjknb+hhland)./y))
fprintf('============================================\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause

clc
if (iprint==0);
fprintf('                  TABLE 3\n')
else
fprintf('                  TABLE 8\n')
end;
fprintf('POPULATION, EMPLOYMENT, & HOURS, AVG %4.0f:%4.0f\n',t0,t1)
fprintf('============================================\n')
fprintf('BLS-CPS aggregates                      \n')
fprintf('  Noninst population, 16-64 (mil)  %7.1f\n',mean(pop1664))
fprintf('  Annual hours per person          %7.1f\n',mean(blsahr))
fprintf('\n')
fprintf('BEA persons engaged (mil)               \n')
fprintf('  Total persons engaged            %7.1f\n',mean(pereng))
fprintf('    FTE employees                  %7.1f\n',mean(ftempl))
fprintf('      Private                      %7.1f\n',mean(ftemplp))
fprintf('      Government                   %7.1f\n',mean(ftemplg))
fprintf('    FTE sole props & partners      %7.1f\n',mean(selfemp))
fprintf('\n')
fprintf('Survey of Business Owners, 2007         \n')
fprintf('  All owners (private business)         \n')
fprintf('    Number (mil)                   %7.1f\n',town/1000)
fprintf('    Annual hours per owner         %7.1f\n',52*ahrs)
fprintf('  Owners, business is primary income    \n')
fprintf('    Number (mil)                   %7.1f\n',townmi/1000)
fprintf('    Annual hours per owner         %7.1f\n',52*ahrsmi)
fprintf('\n')
fprintf('BEA non-bus comp shares (for HH+NP estimates)\n')
fprintf('  Household                        %7.1f\n',mean(dhhcomp./nbuscomp)*100)
fprintf('  Non-profit                       %7.1f\n',mean(dnpcomp./nbuscomp)*100)
fprintf('  General government               %7.1f\n',mean(dggcomp./nbuscomp)*100)
fprintf('  Government enterprises           %7.1f\n',mean(dgebuscomp./nbuscomp)*100)
fprintf('\n')
fprintf('==> if 5200 hours of discretionary time,\n')
fprintf('    then fraction of work time is:      \n')
fprintf('      -- for total population      %7.1f\n',mean(blsahr)/52)
fprintf('      -- for business owners:           \n')
fprintf('           -- all owners           %7.1f\n', ...
  ahrs*town/(mean(pop1664)*1000))
fprintf('           -- bus is primary inc   %7.1f\n', ...
  ahrsmi*townmi/(mean(pop1664)*1000))
fprintf('      -- for nonbusiness:               \n')
fprintf('           -- government           %7.1f\n', ...
  mean(blsahr.*ftemplg./pereng)/52)
fprintf('           -- households+nonprofit %7.1f\n', ...
  mean(blsahr.*ftemplg./pereng .*(dhhcomp+dnpcomp)./(dgebuscomp+dggcomp)/52))
fprintf('============================================\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause

clc
if (iprint==0);
fprintf('                  TABLE 4\n')
else
fprintf('                  TABLE 9\n')
end;
fprintf('PERCENT OF OWNERS BY YEAR OF ACQUISITION, 2007\n')
fprintf('==============================================\n')
fprintf('                        %% of Owners          \n')
fprintf('                 -----------------------------\n')
fprintf('  Year since     All SBO    Owners w/ business\n')
fprintf('  Acquisition    owners     as primary income \n')
fprintf('----------------------------------------------\n')
for i=1:7;
  tem = [acqyrs1(i),acqyrs2(i),acqpct(i),acqpctmi(i)];
  fprintf('  %3.0f-%3.0f    %10.2f   %10.2f \n',tem)
end;
fprintf('==============================================\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause

clc
if (iprint==0);
fprintf('                  TABLE 5\n')
else
fprintf('                  TABLE 10\n')
end;
fprintf('BANK, VENTURE, AND GOVT START-UP FINANCING, 2007\n')
fprintf('==============================================\n')
fprintf('  Fraction of SBO owners           %7.1f\n',amtcap0)
fprintf('    0 to 25K                       %7.1f\n',amtcap1)
fprintf('    25 to 100K                     %7.1f\n',amtcap2)
fprintf('    100K or more                   %7.1f\n',amtcap3)
fprintf('==============================================\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')
pause


clc
if (iprint==0);
fprintf('                  TABLE 6\n')
else
fprintf('                  TABLE 11\n')
end;
fprintf('CONCENTRATION OF BUSINESS PAYROLL, SBO, 2007*\n')
fprintf('============================================\n')
fprintf('  %% of payroll with bottom 80%% of owners  \n')
fprintf('    All businesses                 %7.1f\n',pay80)
fprintf('    All employer businesses        %7.1f\n',pay80eb)
fprintf('============================================\n')
fprintf('* Business with <10 owners sorted by payroll per owner\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')

clc
if (iprint==0);
fprintf('                  TABLE 7\n')
else
fprintf('                  TABLE 12\n')
end;
tab = [page(2:6)/sum(page(2:6))*100,nage(2:6)/sum(nage(2:6))*100, ...
       bage(2:6)];
fprintf('LIFE CYCLE STATISTICS FOR BUSINESSES, SBO, 2007\n')
fprintf('==================================================\n')
fprintf('    Age    %% US Adults  %% SBO Owners   Avg age \n')
fprintf('  Bracket  in bracket    in bracket  of Business\n')
fprintf('   25-34   %7.0f  %10.0f  %13.2f\n',tab(1,:))
fprintf('   35-44   %7.0f  %10.0f  %13.2f\n',tab(2,:))
fprintf('   45-54   %7.0f  %10.0f  %13.2f\n',tab(3,:))
fprintf('   55-64   %7.0f  %10.0f  %13.2f\n',tab(4,:))
fprintf('   65+     %7.0f  %10.0f  %13.2f\n',tab(5,:))
fprintf('==================================================\n')
fprintf('* Source of population distribution is Census\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')

clc
if (iprint==0);
fprintf('                  TABLE 8\n')
else
fprintf('                  TABLE 13\n')
end;
tab = [page(2:6)/sum(page(2:6))*100,nage2(2:6)/sum(nage2(2:6))*100, ...
       bage2(2:6)];
fprintf('LIFE CYCLE STATISTICS FOR BUSINESSES, SBO, 2007\n')
fprintf('==================================================\n')
fprintf('    Age    %% US Adults  %% SBO Owners   Avg time\n')
fprintf('  Bracket  in bracket    in bracket  Running Bus.\n')
fprintf('   25-34   %7.0f  %10.0f  %13.2f\n',tab(1,:))
fprintf('   35-44   %7.0f  %10.0f  %13.2f\n',tab(2,:))
fprintf('   45-54   %7.0f  %10.0f  %13.2f\n',tab(3,:))
fprintf('   55-64   %7.0f  %10.0f  %13.2f\n',tab(4,:))
fprintf('   65+     %7.0f  %10.0f  %13.2f\n',tab(5,:))
fprintf(' \n')
fprintf('Addendum:\n ')
fprintf('  Pct. of exiting owners that sell: \n')
fprintf('     All owners                      %8.2f\n',fracsold1)
fprintf('     Bus is primary inc              %8.2f\n',fracsold2)
fprintf('==================================================\n')
fprintf('* Source of population distribution is Census\n')
fprintf('  Retired and deceased are excluded from exiting group\n')
fprintf(' \n')
fprintf('Hit return to continue...\n')
fprintf(' \n')

clc
if (iprint==0);
fprintf('                  TABLE 9\n')
else
fprintf('                  TABLE 14\n')
end;
tab = [nobyage(2:6)/sum(nobyage(2:6))*100,ahrsbyage(2:6)*52, ...
       nobyagemi(2:6)/sum(nobyagemi(2:6))*100,ahrsbyagemi(2:6)*52];

fprintf('LIFE CYCLE ANNUAL HOURS, SBO, 2007\n')
fprintf('================================================\n')
fprintf('              All Owners      Bus is Primary Inc\n')
fprintf('          ------------------  ------------------\n')
fprintf('    Age   %% Owners    Avg     %% Owners    Avg\n')
fprintf('  Bracket reporting  Hours    reporting  Hours \n')
fprintf('   25-34  %7.0f  %7.0f    %7.0f  %7.0f\n',tab(1,:))
fprintf('   35-44  %7.0f  %7.0f    %7.0f  %7.0f\n',tab(2,:))
fprintf('   45-54  %7.0f  %7.0f    %7.0f  %7.0f\n',tab(3,:))
fprintf('   55-64  %7.0f  %7.0f    %7.0f  %7.0f\n',tab(4,:))
fprintf('   65+    %7.0f  %7.0f    %7.0f  %7.0f\n',tab(5,:))
fprintf('================================================\n')
