load NIPA/income.dat
load NIPA/natinc.dat
load NIPA/corpinc.dat
load NIPA/product.dat
load NIPA/deflator.dat
load NIPA/inventory.dat
load NIPA/hours.dat
load NIPA/ibtdat.dat
load FIXASSET/capqind.dat
load FOF/atabs/corpinv.dat
load FOF/atabs/depcd.dat
load FOF/btabs/hhbal.dat
load FOF/btabs/nfbbal.dat
load FOF/btabs/cbal.dat
load BLS/men16.dat
load BLS/women16.dat
load prescott/prescott.csv
load FIXASSET/capital.dat
capital = [capital(:,5:78),NaN*ones(34,1)];
capqind = [capqind(:,5:78),NaN*ones(34,1)];

t       = income(1,:)';
gdi     = income(2,:)';
dep     = income(3,:)';
ibt     = income(4,:)'+income(5,:)';
comp    = income(6,:)';
prof    = income(7,:)';
gprof   = income(8,:)';
natprof = natinc(6,:)';
pinc    = income(9,:)';
rent    = income(10,:)';
nint    = income(11,:)';
desc    = income(12,:)';
pbexp   = income(13,:)';
nccadj  = income(14,:)';
cgdi    = corpinc(2,:)';
cdep    = corpinc(3,:)';
cibt    = corpinc(4,:)'+corpinc(5,:)';
ccomp   = corpinc(6,:)';
ctax    = corpinc(8,:)';
fedres  = corpinc(11,:)';
ndiv    = corpinc(12,:)';
rete    = corpinc(13,:)';
pbtd    = corpinc(16,:)';
patd    = corpinc(17,:)';
gains   = corpinc(18,:)';
ccadj   = corpinc(9,:)';
cnint   = corpinc(10,:)';
gdp     = product(2,:)';
pce     = product(3,:)';
pced    = product(4,:)';
pcend   = product(5,:)';
pces    = product(6,:)';
gpdi    = product(7,:)';
gpdif   = product(8,:)';
gpdii   = product(9,:)';
netxm   = product(10,:)';
netx    = product(11,:)';
netm    = product(12,:)';
gov     = product(13,:)';
govc    = product(14,:)';
govi    = product(15,:)';
kreal   = capqind.* (capital(:,72)*ones(1,75))/100;
kdefl   = capital./ kreal;
kfix    = capital(2,:)';
kpfix   = capital(7,:)';
kgfix   = capital(12,:)';
kcfix   = capital(19,:)'+capital(29,:)';
kncfix  = capital(22,:)'+capital(30,:)'+kgfix;
kcd     = capital(17,:)';
hhland  = [NaN*ones(23,1);hhbal(4:4:208,4)- ...
           hhbal(4:4:208,44)-hhbal(4:4:208,48)];
nfbland = [NaN*ones(23,1);nfbbal(4:4:208,4)- ...
           nfbbal(4:4:208,33)-nfbbal(4:4:208,34)];
cland   = [NaN*ones(23,1);cbal(4:4:208,4)- ...
           cbal(4:4:208,34)-cbal(4:4:208,35)];
land    = hhland+nfbland+cland;
kinv    = inventory(2,:)'-.9*inventory(3,:)';
kncinv  = .9*inventory(3,:)';
z       = zeros(4,57);
z(:)    = kinv(1:228);
kcinv   = [NaN*ones(18,1);z(4,:)'];
z       = zeros(4,57);
z(:)    = kncinv(1:228);
kncinv  = [NaN*ones(18,1);z(4,:)'];
kc      = kcfix+kcinv+cland;
knc     = kncfix+kncinv+kcd+hhland+nfbland;
k       = kfix+kcd+kcinv+kncinv+kcd+land;
dgdp    = deflator(2,:)';
dpce    = deflator(3,:)';
dpced   = deflator(4,:)';
dpcend  = deflator(5,:)';
dpces   = deflator(6,:)';
dgpdi   = deflator(7,:)';
dgpdif  = deflator(8,:)';
dgpdii  = deflator(9,:)';
dnetxm  = deflator(10,:)';
dnetx   = deflator(11,:)';
dnetm   = deflator(12,:)';
dgov    = deflator(13,:)';
dgovc   = deflator(14,:)';
dgovi   = deflator(15,:)';
ibtonly = ibtdat(2,:)';
stax    = sum(ibtdat([4,5,7,9],:))';
ptax    = ibtdat(8,:)';

cgpdif  = sum(corpinv([6,9,12],:) )'+.1*corpinv(8,:)';
z       = zeros(4,52);
z(:)    = cgpdif(1:length(cgpdif)-1);
cgpdif  = [NaN*ones(23,1);mean(z)'];
cgpdii  = corpinv(16,:)'+.1*corpinv(18,:)';
z       = zeros(4,52);
z(:)    = cgpdii(1:length(cgpdii)-1);
cgpdii  = [NaN*ones(23,1);mean(z)'];

cddep   = depcd(2,:)/1000;
z       = zeros(4,52);
z(:)    = cddep(1:length(cddep)-1);
cddep   = [NaN*ones(23,1);mean(z)'];

z       = zeros(4,45);
z(:)    = prescott(:,9);
hcps    = [NaN*ones(30,1);mean(z)']/5000;
z(:)    = prescott(:,4);
pcps    = [NaN*ones(30,1);mean(z)'];


t0      = find(t==1990);
t1      = find(t==1999);
disp('1990-1999 AVERAGES')
disp('------------------')
disp(sprintf(' Corporate income         %5.3f', ...
  mean(cgdi(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Capital consumption    %5.3f', ...
  mean(cdep(t0:t1)./gdp(t0:t1))))
disp(sprintf('   IBT+transfers-subs.    %5.3f', ...
  mean(cibt(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Compensation           %5.3f', ...
  mean(ccomp(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Profits (adjusted)     %5.3f', ...
  mean(prof(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Profits tax liab.    %5.3f', ...
  mean(ctax(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Federal Reserve      %5.3f', ...
  mean(fedres(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Cap.Cons.Adj.        %5.3f', ...
  mean(ccadj(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Net interest           %5.3f', ...
  mean(cnint(t0:t1)./gdp(t0:t1))))
disp(' ')
disp(sprintf(' Noncorporate income      %5.3f', ...
  mean((gdi(t0:t1)-cgdi(t0:t1))./gdp(t0:t1))))
disp(sprintf('   Capital consumption    %5.3f', ...
  mean((dep(t0:t1)-cdep(t0:t1))./gdp(t0:t1))))
disp(sprintf('   IBT+transfers-subs.    %5.3f', ...
  mean((ibt(t0:t1)-cibt(t0:t1))./gdp(t0:t1))))
disp(sprintf('   Compensation           %5.3f', ...
  mean((comp(t0:t1)+.7*pinc(t0:t1)-ccomp(t0:t1))./gdp(t0:t1))))
disp(sprintf('     Wages and salaries   %5.3f', ...
  mean((comp(t0:t1)-ccomp(t0:t1))./gdp(t0:t1))))
disp(sprintf('     70%% prop. income     %5.3f', ...
  mean(.7*pinc(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Profits                %5.3f', ...
  mean((rent(t0:t1)+.3*pinc(t0:t1)+gprof(t0:t1))./gdp(t0:t1))))
disp(sprintf('     30%% prop. income     %5.3f', ...
  mean(.3*pinc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Rental income        %5.3f', ...
  mean(rent(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Govt. Enterprises    %5.3f', ...
  mean(gprof(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Net interest           %5.3f', ...
  mean((nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1))))
disp(' ')
disp(sprintf(' Statistical desc.        %5.3f', ...
  mean(desc(t0:t1)./gdp(t0:t1))))
disp(' ')
disp(sprintf(' GDI                      %5.3f', ...
  mean(gdi(t0:t1)./gdp(t0:t1))))
disp(' ')
disp(' ')
disp(' Domestic Product')
disp(sprintf('   Private consumption    %5.3f', ...
  mean((pce(t0:t1)-pced(t0:t1))./gdp(t0:t1))))
disp(sprintf('   Govt. consumption      %5.3f', ...
  mean(govc(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Investment             %5.3f', ...
  mean((pced(t0:t1)+govi(t0:t1)+gpdi(t0:t1)+netxm(t0:t1))./gdp(t0:t1))))
disp(sprintf('     Corporate tangible   %5.3f', ...
  mean((cgpdif(t0:t1)+cgpdii(t0:t1))./gdp(t0:t1))))
disp(sprintf('       Fixed              %5.3f', ...
  mean(cgpdif(t0:t1)./gdp(t0:t1))))
disp(sprintf('       Inventory          %5.3f', ...
  mean(cgpdii(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Noncorporate         %5.3f', ...
  mean((pced(t0:t1)+govi(t0:t1)+gpdi(t0:t1)+netxm(t0:t1)- ...
        cgpdif(t0:t1)-cgpdii(t0:t1))./gdp(t0:t1))))
disp(sprintf('       Noncorp fixed      %5.3f', ...
  mean((gpdif(t0:t1)-cgpdif(t0:t1))./gdp(t0:t1))))
disp(sprintf('       Noncorp inventory  %5.3f', ...
  mean((gpdii(t0:t1)-cgpdii(t0:t1))./gdp(t0:t1))))
disp(sprintf('       Consumer durables  %5.3f', ...
  mean(pced(t0:t1)./gdp(t0:t1))))
disp(sprintf('       Government         %5.3f', ...
  mean(govi(t0:t1)./gdp(t0:t1))))
disp(sprintf('       Net Exports       %6.3f', ...
  mean(netxm(t0:t1)./gdp(t0:t1))))
disp(' ')
disp(' Domestic Capital')
disp(sprintf('   Corporate tangible     %5.3f', ...
  mean(kc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Land, nonfarm        %5.3f', ...
  mean(cland(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Noncorporate           %5.3f', ...
  mean(knc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Government           %5.3f', ...
  mean(kgfix(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Consumer durables    %5.3f', ...
  mean(kcd(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Land, HH+Nonprofit   %5.3f', ...
  mean(hhland(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Land, nonfarm bus.   %5.3f', ...
  mean(nfbland(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Other                %5.3f', ...
  mean((knc(t0:t1)-kgfix(t0:t1)-kcd(t0:t1)- ...
        land(t0:t1))./gdp(t0:t1))))
disp(' ')
disp(' Addendum:')
disp(sprintf(' Foreign After-tax Prof.  %5.3f', ...
  mean((natprof(t0:t1)-prof(t0:t1))./gdp(t0:t1))))
disp(sprintf(' Consumer durable dep.    %5.3f', ...
  mean(cddep(t0:t1)./gdp(t0:t1))))
disp(sprintf(' PCE                      %5.3f', ...
  mean(pce(t0:t1)./gdp(t0:t1))))
disp(sprintf(' PCE durable              %5.3f', ...
  mean(pced(t0:t1)./gdp(t0:t1))))
disp(sprintf(' IBT only                 %5.3f', ...
  mean(ibtonly(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Sales+other            %5.3f', ...
  mean(stax(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Property               %5.3f', ...
  mean(ptax(t0:t1)./gdp(t0:t1))))
disp(sprintf(' Intermediate Fin. Ser.   %5.3f', ...
  mean(pbexp(t0:t1)./gdp(t0:t1))))
disp(sprintf(' CCadj, Dom. noncorp.     %5.3f', ...
  mean(nccadj(t0:t1)./gdp(t0:t1))))
disp(sprintf(' Net dividends            %5.3f', ...
  mean(ndiv(t0:t1)./gdp(t0:t1))))
disp(sprintf(' Retained Earnings        %5.3f', ...
  mean(rete(t0:t1)./gdp(t0:t1))))
disp(sprintf(' Dom. BT Profits (wo/adj) %5.3f', ...
  mean(pbtd(t0:t1)./gdp(t0:t1))))
disp(sprintf(' Dom. AT Profits (wo/adj) %5.3f', ...
  mean(patd(t0:t1)./gdp(t0:t1))))
disp(' ')
disp(' ')

i     = .041;
gamma = .02;
eta   = .01;
k1u   = .65;
del1u = 0;
x1u   = ( (1+gamma)*(1+eta)-1+del1u )*k1u;
mu1   = mean(cibt(t0:t1)./ibt(t0:t1));
mu2   = 1-mu1;
cons  = mean(pce(t0:t1)./gdp(t0:t1)) ...
        +i*mean(kcd(t0:t1)./gdp(t0:t1)) ...
        +mean(cddep(t0:t1)./gdp(t0:t1));
tem   = mean(stax(t0:t1)./gdp(t0:t1))/cons;
tauc  = tem/(1-tem);
cd    = mean(pced(t0:t1)./gdp(t0:t1))/(1+tauc);
tauccd= tauc*cd;
taucc = mean(stax(t0:t1)./gdp(t0:t1))-tauccd;


disp('1990-1999 ADJUSTED AVERAGES')
disp('---------------------------')
disp(sprintf(' Corporate income         %5.3f', ...
  mean(cgdi(t0:t1)./gdp(t0:t1)) ...
    -mu1* mean(stax(t0:t1)./gdp(t0:t1))  ...
    +x1u-mean(pbexp(t0:t1)./gdp(t0:t1)) ))
disp(sprintf('   Capital consumption    %5.3f', ...
  mean(cdep(t0:t1)./gdp(t0:t1))))
disp(sprintf('   IBT+transfers-subs.    %5.3f', ...
  mean(cibt(t0:t1)./gdp(t0:t1)) ...
    -mu1* mean(stax(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Compensation           %5.3f', ...
  mean(ccomp(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Profits (adjusted)     %5.3f', ...
  mean((prof(t0:t1)+cnint(t0:t1)-pbexp(t0:t1))./ ...
       gdp(t0:t1))+x1u))
disp(sprintf('     Profits tax liab.    %5.3f', ...
  mean(ctax(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Federal Reserve      %5.3f', ...
  mean(fedres(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Cap.Cons.Adj.        %5.3f', ...
  mean(ccadj(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Net interest         %5.3f', ...
  mean((cnint(t0:t1)-pbexp(t0:t1))./gdp(t0:t1))))
disp(' ')
disp(sprintf(' Noncorporate income      %5.3f', ...
  mean((gdi(t0:t1)-cgdi(t0:t1))./gdp(t0:t1)) ...
     -mu2* mean(stax(t0:t1)./gdp(t0:t1))  ...
     +mean(cddep(t0:t1)./gdp(t0:t1)) ...
     +i*mean((kgfix(t0:t1)+kcd(t0:t1))./gdp(t0:t1)) ...
     -0*mean((nint(t0:t1)-cnint(t0:t1))./ ...
        gdp(t0:t1)) ...
     +mean(desc(t0:t1)./gdp(t0:t1)) ))
disp(sprintf('   Capital consumption    %5.3f', ...
  mean((dep(t0:t1)-cdep(t0:t1)+cddep(t0:t1))./ ...
        gdp(t0:t1))))
disp(sprintf('   IBT+transfers-subs.    %5.3f', ...
  mean((ibt(t0:t1)-cibt(t0:t1))./gdp(t0:t1)) ...
     -mu2* mean(stax(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Compensation           %5.3f', ...
  mean((comp(t0:t1)+.7*pinc(t0:t1)-ccomp(t0:t1) ...
     +desc(t0:t1))./gdp(t0:t1))))
disp(sprintf('     Wages and salaries   %5.3f', ...
  mean((comp(t0:t1)-ccomp(t0:t1))./gdp(t0:t1))))
disp(sprintf('     70%% prop. income     %5.3f', ...
  mean(.7*pinc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Statistical desc.    %5.3f', ...
  mean(desc(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Profits                %5.3f', ...
  mean((rent(t0:t1)+.3*pinc(t0:t1)+ ...
        gprof(t0:t1))./gdp(t0:t1)) ...
     +i*mean((kgfix(t0:t1)+kcd(t0:t1))./ ...
        gdp(t0:t1)) ...
 +mean((nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1))))
disp(sprintf('     30%% prop. income     %5.3f', ...
  mean(.3*pinc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Rental income        %5.3f', ...
  mean(rent(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Govt. enterprises    %5.3f', ...
  mean(gprof(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Imputed services     %5.3f', ...
  i*mean((kgfix(t0:t1)+kcd(t0:t1))./gdp(t0:t1))))
disp(sprintf('     Net interest         %5.3f', ...
  mean((nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1))))
disp(' ')
disp(sprintf(' Total income             %5.3f', ...
  mean(cgdi(t0:t1)./gdp(t0:t1)) ...
    -mu1* mean(stax(t0:t1)./gdp(t0:t1))  ...
    +x1u-mean(pbexp(t0:t1)./gdp(t0:t1)) ...
 +mean((gdi(t0:t1)-cgdi(t0:t1))./gdp(t0:t1)) ...
     -mu2* mean(stax(t0:t1)./gdp(t0:t1))  ...
     +mean(cddep(t0:t1)./gdp(t0:t1)) ...
     +i*mean((kgfix(t0:t1)+kcd(t0:t1))./ ...
        gdp(t0:t1)) ...
     -0*mean((nint(t0:t1)-cnint(t0:t1))./ ...
        gdp(t0:t1)) ...
     +mean(desc(t0:t1)./gdp(t0:t1)) ))
disp(' ')
disp(' ')
disp(' Domestic Product')
disp(sprintf('   Private consumption    %5.3f', ...
  mean((pce(t0:t1)-pced(t0:t1))./gdp(t0:t1)) ...
    +i*mean(kcd(t0:t1)./gdp(t0:t1)) ...
    +mean(cddep(t0:t1)./gdp(t0:t1)) ...
    -taucc  ...
    -mean(pbexp(t0:t1)./gdp(t0:t1)) ...
    -0*mean((nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1))))
disp(sprintf('   Govt. consumption      %5.3f', ...
  mean(govc(t0:t1)./gdp(t0:t1)) ...
    +i*mean(kgfix(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Investment             %5.3f', ...
  mean((pced(t0:t1)+govi(t0:t1)+gpdi(t0:t1) ...
    +netxm(t0:t1))./gdp(t0:t1)) ...
    -tauccd+x1u))
disp(sprintf('     Corporate tangible   %5.3f', ...
  mean((cgpdif(t0:t1)+cgpdii(t0:t1))./gdp(t0:t1))))
disp(sprintf('     Corporate intangible %5.3f',x1u))
disp(sprintf('     Noncorporate         %5.3f', ...
  mean((pced(t0:t1)+govi(t0:t1)+gpdi(t0:t1)+netxm(t0:t1)- ...
        cgpdif(t0:t1)-cgpdii(t0:t1))./gdp(t0:t1))-tauccd))
disp(' ')
disp(sprintf(' Total product            %5.3f', ...
  mean((pce(t0:t1)-pced(t0:t1))./gdp(t0:t1)) ...
    +i*mean(kcd(t0:t1)./gdp(t0:t1)) ...
    +mean(cddep(t0:t1)./gdp(t0:t1)) ...
    -taucc  ...
    -mean(pbexp(t0:t1)./gdp(t0:t1)) ...
    -0*mean((nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1)) ...
 +mean(govc(t0:t1)./gdp(t0:t1)) ...
    +i*mean(kgfix(t0:t1)./gdp(t0:t1)) ...
 +mean((pced(t0:t1)+govi(t0:t1)+gpdi(t0:t1) ...
    +netxm(t0:t1))./gdp(t0:t1)) ...
    -tauccd+x1u));
   
disp(' ')
disp(' Domestic Capital')
disp(sprintf('   Corporate tangible     %5.3f', ...
  mean(kc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Land, nonfarm        %5.3f', ...
  mean(cland(t0:t1)./gdp(t0:t1))))
disp(sprintf('   Noncorporate           %5.3f', ...
  mean(knc(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Government           %5.3f', ...
  mean(kgfix(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Consumer durables    %5.3f', ...
  mean(kcd(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Land, HH+Nonprofit   %5.3f', ...
  mean(hhland(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Land, nonfarm bus.   %5.3f', ...
  mean(nfbland(t0:t1)./gdp(t0:t1))))
disp(sprintf('     Other                %5.3f', ...
  mean((knc(t0:t1)-kgfix(t0:t1)-kcd(t0:t1)- ...
        land(t0:t1))./gdp(t0:t1))))
disp(' ')
disp(' Per-capita hours')
disp(sprintf('   Total                  %5.3f', ...
  mean(hcps(t0:t1))))
disp(sprintf('   Corporate              %5.3f', ...
  mean(hcps(t0:t1)./(1+(comp(t0:t1)+.7*pinc(t0:t1)- ...
       ccomp(t0:t1)+desc(t0:t1))./ccomp(t0:t1))) ))
disp(sprintf('   Noncorporate           %5.3f', ...
  mean(hcps(t0:t1)-hcps(t0:t1)./(1+(comp(t0:t1)+ ...
       .7*pinc(t0:t1)-ccomp(t0:t1)+desc(t0:t1))./ ...
       ccomp(t0:t1))) ))
