%ACCOUNTS   Tables in "On Financing Retirement with an Aging Population"
%

%           Ellen McGrattan, 11-17-11
%           Revised, ERM, 9-4-12


%
% Sample t0:t1
%

t0           = 2000;
t1           = 2010;

%
% NIPA Table 1.1.5 Gross domestic product (Annual, Bil$)
%

load BEA/NIPA/table1_1_5A.dat
[n,m]        = size(table1_1_5A);
t            = table1_1_5A(1,:);
i            = find(t==t0):find(t==t1);
table1_1_5A  = table1_1_5A(2:n,i);

gdp          = table1_1_5A(1,:)';
pce          = table1_1_5A(2,:)';
pced         = table1_1_5A(4,:)';
pcend        = table1_1_5A(5,:)';
pces         = table1_1_5A(6,:)';
gpdi         = table1_1_5A(7,:)';
gpdif        = table1_1_5A(8,:)';
gpdii        = table1_1_5A(13,:)';
netxm        = table1_1_5A(14,:)';
netx         = table1_1_5A(15,:)';
netm         = table1_1_5A(18,:)';
gov          = table1_1_5A(21,:)';


%
% NIPA Table 1.1.6 Real GDP, Chained 2005 Dollars
%

load BEA/NIPA/table1_1_6A.dat
[n,m]        = size(table1_1_6A);
t            = table1_1_6A(1,:);
i            = find(t==t0):find(t==t1);
table1_1_6A  = table1_1_6A(2:n,i);

rgdp         = table1_1_6A(1,:)';
rpce         = table1_1_6A(2,:)';
rgpdi        = table1_1_6A(7,:)';
rnetxm       = table1_1_6A(14,:)';
rnetx        = table1_1_6A(15,:)';
rnetm        = table1_1_6A(18,:)';

%
% NIPA Table 1.7.5. Relation of GDP, GNP, NNP, NI, PI (Annual, Bil$)
%

load BEA/NIPA/table1_7_5A.dat
[n,m]        = size(table1_7_5A);
t            = table1_7_5A(1,:);
i            = find(t==t0):find(t==t1);
table1_7_5A  = table1_7_5A(2:n,i);

nfir         = table1_7_5A(2,:)';
nfip         = table1_7_5A(3,:)';          
gnp          = table1_7_5A(4,:)';
row          = nfir-nfip;


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

load BEA/NIPA/table1_10A.dat
[n,m]        = size(table1_10A);
t            = table1_10A(1,:);
i            = find(t==t0):find(t==t1);
table1_10A   = table1_10A(2:n,i);

gdi          = table1_10A(1,:)';
comp         = table1_10A(2,:)';
wage         = table1_10A(3,:)';
supp         = table1_10A(8,:)';
txpi         = table1_10A(9,:)';
subs         = table1_10A(10,:)';
bustra       = table1_10A(14,:)';
ibt          = txpi-subs+bustra;
nint         = table1_10A(13,:)';
pinc         = table1_10A(15,:)';
rent         = table1_10A(16,:)';
prof         = table1_10A(17,:)';
ctax         = table1_10A(18,:)';
ddiv         = table1_10A(20,:)';
dre          = table1_10A(21,:)';
gprof        = table1_10A(22,:)';
dep          = table1_10A(23,:)';
desc         = table1_10A(26,:)';
labinc       = comp+.7*pinc;

%
% NIPA Table 1.12 National Income by Type of Income (Annual, Bil$)
%

load BEA/NIPA/table1_12A.dat
[n,m]        = size(table1_12A);
t            = table1_12A(1,:);
i            = find(t==t0):find(t==t1);
table1_12A   = table1_12A(2:n,i);

ncomp        = table1_12A(2,:)';
nprof        = table1_12A(13,:)';
nctax        = table1_12A(14,:)';
ndiv         = table1_12A(16,:)';
nre          = table1_12A(17,:)';
nnint        = table1_12A(18,:)';

%
% NIPA Table 1.14  Gross VA of domestic corporate business (Annual, Bil$)
%

load BEA/NIPA/table1_14A.dat
[n,m]        = size(table1_14A);
t            = table1_14A(1,:);
i            = find(t==t0):find(t==t1);
table1_14A   = table1_14A(2:n,i);

cgdi         = table1_14A(1,:)';
cdep         = table1_14A(2,:)';
ccomp        = table1_14A(4,:)';
cibt         = table1_14A(7,:)'+table1_14A(10,:)';
cnint        = table1_14A(9,:)';
cprof        = table1_14A(11,:)';
ctax         = table1_14A(12,:)';
ddiv2        = table1_14A(14,:)';
rete         = table1_14A(15,:)';
pbtd         = table1_14A(32,:)';
patd         = table1_14A(33,:)';
ccadj        = table1_14A(35,:)';

%
% NIPA Table 3.1 Govt current receipts and expenditures (Annual, Bil$)
%

load BEA/NIPA/table3_1A.dat
[n,m]        = size(table3_1A);
t            = table3_1A(1,:);
i            = find(t==t0):find(t==t1);
t3_1A        = table3_1A(1,i)';
table3_1A    = table3_1A(2:n,i);

hhdt         = table3_1A(3,:)';
hhss         = table3_1A(7,:)';

govr         = table3_1A(30,:)';
govtax       = table3_1A(2,:)';
govtaxpc     = table3_1A(3,:)';
govtaxibt    = table3_1A(4,:)';
govtaxci     = table3_1A(5,:)';
govtaxrow    = table3_1A(6,:)';
govrss       = table3_1A(7,:)';
govrin       = table3_1A(8,:)';
govrge       = table3_1A(14,:)';
govrtrcur    = table3_1A(11,:)';
govrtrcap    = table3_1A(32,:)';
govrtr       = govrtrcur+govrtrcap;

gove         = table3_1A(33,:)';
govc         = table3_1A(16,:)';
govi         = table3_1A(35,:)';
govetrcur    = table3_1A(17,:)';
govetrsub    = table3_1A(25,:)';
govetrcap    = table3_1A(36,:)';
govenp       = table3_1A(37,:)';
govetr       = govetrcur+govetrcap+govetrsub;
govein       = table3_1A(22,:)';

govdep       = table3_1A(38,:)';
deficit      = -table3_1A(39,:)';


%
% NIPA Table 3.2 Federal govt current receipts and expenditures (Annual, Bil$)
%

load BEA/NIPA/table3_2A.dat
[n,m]        = size(table3_2A);
t            = table3_2A(1,:);
i            = find(t==t0):find(t==t1);
table3_2A    = table3_2A(2:n,i);

fedgove      = table3_2A(40,:)';
fedres       = table3_2A(8,:)';

%
% NIPA Table 3.3 S&Lgovt current receipts and expenditures (Annual, Bil$)
%

load BEA/NIPA/table3_3A.dat
[n,m]        = size(table3_3A);
t            = table3_3A(1,:);
i            = find(t==t0):find(t==t1);
table3_3A    = table3_3A(2:n,i);

slgove       = table3_3A(33,:)';

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

load BEA/NIPA/table3_5A.dat
[n,m]        = size(table3_5A);
t            = table3_5A(1,:);
i            = find(t==t0):find(t==t1);
table3_5A    = table3_5A(2:n,i);

stax         = sum(table3_5A([3,11,12,14,28,29,30,31],:))';
fedexcise    = table3_5A(3,:)';
fedcustom    = table3_5A(11,:)';
slsales      = table3_5A(14,:)';
slmvl        = table3_5A(28,:)';
slsev        = table3_5A(29,:)';
slspec       = table3_5A(30,:)';
txpioth      = table3_5A(12,:)'+table3_5A(31,:)';
ptax         = table3_5A(27,:)';

%
% NIPA Table 3.9.5 Govt consumption exp. and gross investment (Annual, Bil$)
%

load BEA/NIPA/table3_9_5A.dat
[n,m]        = size(table3_9_5A);
t            = table3_9_5A(1,:);
i            = find(t==t0):find(t==t1);
table3_9_5A  = table3_9_5A(2:n,i);

govcdef      = table3_9_5A(12,:)';
govidef      = table3_9_5A(13,:)';
govdef       = govcdef+govidef;

%
% NIPA Table 3.12 Govt social benefits (Annual, Bil$)
%
load BEA/NIPA/table3_12A.dat
[n,m]        = size(table3_12A);
t            = table3_12A(1,:);
i            = find(t==t0):find(t==t1);
table3_12A   = table3_12A(2:n,i);

govetrsb     = table3_12A(1,:)';
govetrsbss   = table3_12A(5,:)';
govetrsbmc   = table3_12A(6,:)';
govetrsbor   = sum(table3_12A([12,13,14,17,22,23,36],:))';
govetrsbui   = table3_12A(7,:)';
govetrsbwc   = table3_12A(15,:)'+table3_12A(28,:)';
govetrsbpa   = sum(table3_12A([21,25,31],:))'-table3_12A(36,:)';
govetrsbed   = table3_12A(40,:)'+table3_12A(41,:)';
govetrsbo    = govetrcur-govetrsbss-govetrsbmc-govetrsbor-govetrsbui- ...
               govetrsbwc-govetrsbpa-govetrsbed;
govetroth    = govetrcur-govetrsbss-govetrsbmc;

%
% NIPA Table 3.13 Govt subsidies (Annual, Bil$)
%
load BEA/NIPA/table3_13A.dat
[n,m]        = size(table3_13A);
t            = table3_13A(1,:);
i            = find(t==t0):find(t==t1);
table3_13A   = table3_13A(2:n,i);

govetrsubag  = table3_13A(3,:)';
govetrsubhh  = table3_13A(4,:)';
govetrsubo   = sum(table3_13A(5:8,:))';

%
% NIPA Table 3.15.5 Govt consumption and gross inv. by function (Annual, Bil$)
%
load BEA/NIPA/table3_15_5A.dat
[n,m]        = size(table3_15_5A);
t            = table3_15_5A(1,:);
i            = find(t==t0):find(t==t1);
table3_15_5A = table3_15_5A(2:n,i);

govtotal     = table3_15_5A(1,:)';
genpubser    = table3_15_5A(2,:)';
natdefense   = table3_15_5A(6,:)';
puborder     = table3_15_5A(7,:)';
econaffairs  = table3_15_5A(12,:)';
housing      = table3_15_5A(26,:)';
health       = table3_15_5A(27,:)';
recreation   = table3_15_5A(28,:)';
education    = table3_15_5A(29,:)';
welfare      = table3_15_5A(35,:)';

%
% NIPA Table 3.17 Govt consumption and gross inv. by function (Annual, Bil$)
%
load BEA/NIPA/table3_17A.dat
[n,m]        = size(table3_17A);
t            = table3_17A(1,:);
i            = find(t==t0):find(t==t1);
table3_17A   = table3_17A(2:n,i);

govtotalc    = table3_17A(1,:)';
genpubserc   = table3_17A(2,:)';
natdefensec  = table3_17A(3,:)';
puborderc    = table3_17A(4,:)';
econaffairsc = table3_17A(5,:)';
housingc     = table3_17A(6,:)';
healthc      = table3_17A(7,:)';
recreationc  = table3_17A(8,:)';
educationc   = table3_17A(9,:)';
welfarec     = table3_17A(10,:)';
govtotali    = table3_17A(98,:)';
genpubseri   = table3_17A(99,:)';
natdefensei  = table3_17A(100,:)';
puborderi    = table3_17A(101,:)';
econaffairsi = table3_17A(102,:)';
housingi     = table3_17A(103,:)';
healthi      = table3_17A(104,:)';
recreationi  = table3_17A(105,:)';
educationi   = table3_17A(106,:)';
welfarei     = table3_17A(107,:)';


%
% NIPA Table 4.1 Foreign transactions in NIPA (Annual, Bil$)
%

load BEA/NIPA/table4_1A.dat
[n,m]        = size(table4_1A);
t            = table4_1A(1,:);
i            = find(t==t0):find(t==t1);
table4_1A    = table4_1A(2:n,i);

rowinvrec    = table4_1A(9,:)';
rowinvpay    = table4_1A(21,:)';

%
% NIPA Table 5.6.5A Change in private inventories by industry (Annual, Bil$)
%

load BEA/NIPA/table5_6_5AA.dat
[n,m]        = size(table5_6_5AA);
load BEA/NIPA/table5_6_5BA.dat
table5_6_5A  = [table5_6_5AA(1:3,1:m-1),table5_6_5BA(1:3,:)];
t            = table5_6_5A(1,:);
i            = find(t==t0):find(t==t1);
table5_6_5A  = table5_6_5A(2:3,i);

gpdiifarm    = table5_6_5A(2,:)';
gpdiinfarm   = gpdii-gpdiifarm;
gpdiicfarm   = .13*gpdiifarm;
               %Note: the fraction of farm inventories assigned to the
               %      corporate sector is 13% based on the ratio of corporate
               %      farmland and buildings to total values in the U.S.
               %      Statistical Abstract 2012

%
% NIPA Table 5.7.5A Private inventories and final sales (Quarterly, Bil$)
%

load BEA/NIPA/table5_7_5BQ.dat
load BEA/NIPA/table5_7_5AQ.dat
table5_7_5BQ = [table5_7_5BQ(1,:)+table5_7_5BQ(2,:)/10;
                table5_7_5BQ(3:4,:)];
table5_7_5AQ = [table5_7_5AQ(1,:)+table5_7_5AQ(2,:)/10;
                table5_7_5AQ(3:4,:)];
i            = find(table5_7_5AQ(1,:)==table5_7_5BQ(1,1));
table5_7_5Q  = [table5_7_5AQ(1:3,1:i-1),table5_7_5BQ(1:3,:)];
table5_7_5A  = (qtr2ann(table5_7_5Q',1))';
[n,m]        = size(table5_7_5A);
t            = table5_7_5A(1,:);
i            = find(t==t0):find(t==t1);
table5_7_5A  = table5_7_5A(2:n,i);

kcinv        = table5_7_5A(1,:)'-.87*table5_7_5A(2,:)';
kncinv       = .87*table5_7_5A(2,:)';

%
% NIPA Table 7.16 Relation of NIPA and IRS Corporate Statistics (Annual, Bil$)
%
load BEA/NIPA/table7_16A.dat
[n,m]        = size(table7_16A);
t            = table7_16A(1,:);
i            = find(t==t0):find(t==t1);
table7_16A   = table7_16A(2:n,i);

divirs       = table7_16A(30,:)';
divnipa      = table7_16A(38,:)';    


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

load BEA/FA/table1_1A.dat
[n,m]        = size(table1_1A);
t            = table1_1A(1,:);
i            = find(t==t0):find(t==t1);
table1_1A    = table1_1A(2:n,i);

kfix         = table1_1A(2,:)';
kfixres      = table1_1A(7,:)';
kgfix        = table1_1A(8,:)';
kcd          = table1_1A(13,:)';

%
% FA Table 1.5 Historical-cost investments (Annual, Bil$, year-end) 
%

load BEA/FA/table1_5A.dat
[n,m]        = size(table1_5A);
t            = table1_5A(1,:);
i            = find(t==t0):find(t==t1);
table1_5A    = table1_5A(2:n,i);

xfix         = table1_5A(2,:)';
xgfix        = table1_5A(8,:)';
xcd          = table1_5A(13,:)';
xpfix        = xfix-xgfix;

%
% FA Tables 6.1 Current-cost net stocks by owner (Annual, Bil$, year-end) 
%

load BEA/FA/table6_1A.dat
[n,m]        = size(table6_1A);
t            = table6_1A(1,:);
i            = find(t==t0):find(t==t1);
table6_1A    = table6_1A(2:n,i);

kcfix        = table6_1A(2,:)';
kppfix       = table6_1A(6,:)'+table6_1A(7,:)';
khhfix       = table6_1A(9,:)';
knpfix       = table6_1A(8,:)';
ktefix       = table6_1A(10,:)';
kncfix       = kppfix+khhfix+knpfix+ktefix;

%
% FA Table 6.7 Historical-cost investments by owner (Annual, Bil$, year-end) 
%

load BEA/FA/table6_7A.dat
[n,m]        = size(table6_7A);
t            = table6_7A(1,:);
i            = find(t==t0):find(t==t1);
table6_7A    = table6_7A(2:n,i);

xcfix        = table6_7A(2,:)';
xppfix       = table6_7A(6,:)'+table6_7A(7,:)';
xhhfix       = table6_7A(9,:)';
xnpfix       = table6_7A(8,:)';
xtefix       = table6_7A(10,:)';
xncfix       = xppfix+xhhfix+xnpfix+xtefix;

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

load BEA/FA/table8_4A.dat
[n,m]        = size(table8_4A);
t            = table8_4A(1,:);
i            = find(t==t0):find(t==t1);
table8_4A    = table8_4A(2:n,i);

dcddep       = table8_4A(1,:)';

%
% ICI Defined contribution plan holdings (Annual, Bil$)
%
load ICI/dc_mf_type.dat
[n,m]        = size(dc_mf_type);
t            = dc_mf_type(:,1);
i            = find(t==t0):find(t==t1);
dcmftype     = dc_mf_type(i,2:m);

dcmfdomeq    = dcmftype(:,1);
dcmfforeq    = dcmftype(:,3);
dcmfhybrid   = dcmftype(:,5);
dcforeqrat   = dcmfforeq./(dcmfdomeq+dcmfforeq+.5*dcmfhybrid);


%
% ICI IRA plan holdings (Annual, Bil$)
%
load ICI/ira_mf_type.dat
[n,m]        = size(ira_mf_type);
t            = ira_mf_type(:,1);
i            = find(t==t0):find(t==t1);
iramftype    = ira_mf_type(i,2:m);

iramfdomeq   = iramftype(:,1);
iramfforeq   = iramftype(:,3);
iramfhybrid  = iramftype(:,5);
iramftotal   = iramftype(:,11);
iraforeqrat  = iramfforeq./(iramfdomeq+iramfforeq+.5*iramfhybrid);

%
% FOF/ICI IRA equity share (%)
%
load ICI/ira_eq_share.dat
[n,m]        = size(ira_eq_share);
t            = ira_eq_share(:,1);
i            = find(t==t0):find(t==t1);
ira_eq_share = ira_eq_share(i,2);


%
% FOF Table F6 Gross domestic product (Quarterly, Mil$, Averaged)
%

load FOF/atabs/atab6d.dat
atab6dA      = (qtr2ann(atab6d))';
[n,m]        = size(atab6dA);
t            = atab6dA(1,:);
i            = find(t==t0):find(t==t1);
atab6dA      = atab6dA(2:n,i)/1000;


cgpdif       = sum(atab6dA([10,12,15,17],:))';
npgpdif      = atab6dA(9,:)';
hhgpdif      = atab6dA(14,:)';
ncbgpdif     = sum(atab6dA([11,16],:))';
ncgpdif      = npgpdif+hhgpdif+ncbgpdif;
fofgpdif     = cgpdif+ncgpdif;

cgpdii       = atab6dA(19,:)';
ncbgpdii     = atab6dA(20,:)';
fofgpdii     = cgpdii+ncbgpdii;
fofgpdi      = fofgpdif+fofgpdii;
cgpdi        = cgpdif+cgpdii;

fofgdp       = atab6dA(1,:)';
fofpce       = atab6dA(2,:)';
fofpced      = atab6dA(3,:)';
fofpcend     = atab6dA(4,:)';
fofpces      = atab6dA(5,:)';
fofnx        = atab6dA(21,:)';
fofgovc      = atab6dA(25,:)';
fofgovi      = atab6dA(28,:)';

%
% FOF Table F10 Measures of personal saving (Quarterly, Mil$)
%
load FOF/atabs/atab10d.dat
atab10dA     = (qtr2ann(atab10d))';
[n,m]        = size(atab10dA);
t            = atab10dA(1,:);
i            = find(t==t0):find(t==t1);
atab10dA     = atab10dA(2:n,i)/1000;

cddep        = atab10dA(27,:)';

%
% FOF Table B100 Balance sheet of households (Quarterly, Mil$)
%

load FOF/gtabs/gtab3d.dat
gtab3dA      = (qtr2ann(gtab3d,1))';
[n,m]        = size(gtab3dA);
t            = gtab3dA(1,:);
i            = find(t==t0):find(t==t1);
gtab3dA      = gtab3dA(2:n,i)/1000;

govdebt      = gtab3dA(7,:)'+gtab3dA(8,:)';

%
% FOF Table B100 Balance sheet of households (Quarterly, Mil$)
%

load FOF/btabs/btab100d.dat
btab100dA    = (qtr2ann(btab100d,1))';
[n,m]        = size(btab100dA);
t            = btab100dA(1,:);
i            = find(t==t0):find(t==t1);
btab100dA    = btab100dA(2:n,i)/1000;

vhha         = btab100dA(1,:)';
vhhtang      = btab100dA(2,:)';
vhhfa        = btab100dA(8,:)';
vhhli        = btab100dA(31,:)';
hhland       = btab100dA(3,:)'-btab100dA(43,:)'-btab100dA(46,:)';
nwhh         = btab100dA(42,:)';

%
% FOF Table B102 Balance sheet of nonfarm nonfin. corporate (Quarterly, Mil$)
%

load FOF/btabs/btab102d.dat
btab102dA    = (qtr2ann(btab102d,1))';
[n,m]        = size(btab102dA);
t            = btab102dA(1,:);
i            = find(t==t0):find(t==t1);
btab102dA    = btab102dA(2:n,i)/1000;

vnfctang     = btab102dA(2,:)';
nfcland      = btab102dA(3,:)'-btab102dA(33,:)'-btab102dA(34,:)';

%
% FOF Table B103 Balance sheet of nonfarm noncorporate (Quarterly, Mil$)
%

load FOF/btabs/btab103d.dat
btab103dA    = (qtr2ann(btab103d,1))';
[n,m]        = size(btab103dA);
t            = btab103dA(1,:);
i            = find(t==t0):find(t==t1);
btab103dA    = btab103dA(2:n,i)/1000;

vnfnctang    = btab103dA(2,:)';
nfncland     = btab103dA(3,:)'-btab103dA(33,:)'-btab103dA(34,:)';
vbustang     = vnfctang+vnfnctang;

land         = hhland+nfcland+nfncland;
kinv         = kcinv+kncinv;
k            = kfix+kinv+kcd+land;
kc           = kcfix+kcinv+nfcland;

%
% FOF Table L104 State and Local Governments (Quarterly, Mil$)
%
load FOF/ltabs/ltab104d.dat
ltab104dA    = (qtr2ann(ltab104d,1))';
[n,m]        = size(ltab104dA);
t            = ltab104dA(1,:);
i            = find(t==t0):find(t==t1);
ltab104dA    = ltab104dA(2:n,i)/1000;

slmuni       = ltab104dA(20,:)';
sla          = ltab104dA(1,:)';
slli         = ltab104dA(18,:)';

%
% FOF Table L105 Federal Government (Quarterly, Mil$)
%
load FOF/ltabs/ltab105d.dat
ltab105dA    = (qtr2ann(ltab105d,1))';
[n,m]        = size(ltab105dA);
t            = ltab105dA(1,:);
i            = find(t==t0):find(t==t1);
ltab105dA    = ltab105dA(2:n,i)/1000;

fedtrs       = ltab105dA(20,:)'+ltab105dA(21,:)'+ltab105dA(22,:)';
fedagcy      = ltab105dA(23,:)';
feda         = ltab105dA(1,:)';
fedli        = ltab105dA(15,:)';
govsec       = slmuni+fedtrs+fedagcy;

%
% FOF Table L213 Corporate Equities (Quarterly, Mil$)
%
load FOF/ltabs/ltab213d.dat
ltab213dA    = (qtr2ann(ltab213d,1))';
[n,m]        = size(ltab213dA);
t            = ltab213dA(1,:);
i            = find(t==t0):find(t==t1);
ltab213dA    = ltab213dA(2:n,i)/1000;

dommktval    = ltab213dA(23,:)';

%
% FOF Table L225 Life Insurance and Pension Fund Reserves (Quarterly, Mil$)
%
load FOF/ltabs/ltab225d.dat
ltab225dA    = (qtr2ann(ltab225d,1))';
[n,m]        = size(ltab225dA);
t            = ltab225dA(1,:);
i            = find(t==t0):find(t==t1);
ltab225dA    = ltab225dA(2:n,i)/1000;

liannuity    = ltab225dA(3,:)';

%
% FOF Table L225i Individual Retirement Accounts (Annual, Mil$)
%
load FOF/stabs/ltab225i.dat
[n,m]        = size(ltab225i);
t            = ltab225i(:,1);
i            = find(t==t0):find(t==t1);
ltab225iA    = ltab225i(i,2:m)/1000;

ira          = ltab225iA(:,1);
irali        = ltab225iA(:,4);



%
% FOF Table B100e Balance sheet of households with equity detail (Annual, Mil$)
%
load FOF/stabs/btab100e.dat
[n,m]        = size(btab100e);
t            = btab100e(:,1);
i            = find(t==t0):find(t==t1);
btab100eA    = btab100e(i,2:m)/1000;

hhequity     = btab100eA(:,6);
ppequity     = btab100eA(:,10);
gpequity     = btab100eA(:,13)+btab100eA(:,14);
iraequity    = ira_eq_share.*(ira-irali)/100;
ppequityd    = ppequity.*(1-dcforeqrat);
iraequityd   = iraequity.*(1-iraforeqrat);
equntaxed    = (ppequity+gpequity+liannuity+iraequity)./hhequity;
equntaxedd   = (ppequityd+gpequity+liannuity+iraequityd)./dommktval;

%
% SCB BEA dividend income vs. IRS dividend income (Annual, Bil$)
%
load BEA/SCB/nipadiv.dat
[n,m]        = size(nipadiv);
t            = nipadiv(:,1);
i            = find(t==t0):find(t==t1);
nipadiv      = nipadiv(i,2:m);

divtaxed     = nipadiv(:,13)-nipadiv(:,10);
divuntaxed   = sum(nipadiv(:,5:8)')';
divtotal     = divtaxed+divuntaxed;
divtaxed     = divtaxed./divtotal;
divuntaxed   = divuntaxed./divtotal;

%
% IRS Business net income less deficit (Annual, Thous$)
%
load IRS/bustax.dat
[n,m]        = size(bustax);
t            = bustax(:,1);
i            = find(t==t0):find(t==t1);
bustax       = bustax(i,2:m)/10^9;

ni_corp      = bustax(:,3);
ni_ccorp     = bustax(:,7);
ni_scorp     = bustax(:,3)-bustax(:,7);  % includes S, RIC, REIT
da_corp      = bustax(:,1)-bustax(:,2);
da_ccorp     = bustax(:,5)-bustax(:,6);
da_scorp     = da_corp-da_ccorp;
ni_sshare    = ni_scorp./ni_corp;
ni_cshare    = ni_ccorp./ni_corp;
kc_sshare    = da_scorp./da_corp;
kc_cshare    = da_ccorp./da_corp;

%
% IRS Ordinary and qualified dividends (Annual starting 2003, $bil)
%
load IRS/qualdiv.dat
[n,m]        = size(qualdiv);
t            = qualdiv(:,1);
i            = find(t==t0):find(t==t1);
qualdiv      = qualdiv(i,2:m);

qualdivrat   = qualdiv(:,2)./qualdiv(:,1);

%
% IRS Wage-like tax rates of Barro-Redlick using SOI data
%
load IRS/barro.dat
[n,m]        = size(barro);
t            = barro(:,1);
i            = find(t==t0):find(t==t1);
barro        = barro(i,2:m)/100;

taulfed      = barro(:,2);
taulsl       = barro(:,3)-barro(:,2);
taulfica     = barro(:,4);
taul         = taulfed+taulsl+taulfica;

%
% CPS Hours 
%
load CPS/avghrsa.dat
[n,m]        = size(avghrsa);
t            = avghrsa(:,1);
i            = find(t==t0):find(t==t1);
avghrsa      = avghrsa(i,8);

%
% CPS Employment
%
load CPS/employed.dat
[n,m]        = size(employed);
t            = employed(:,1);
i            = find(t==t0):find(t==t1);
employed     = employed(i,2);

%
% BEA Full-time equivalent employees
%
load BEA/NIPA/fteemp.dat
[n,m]        = size(fteemp);
t            = fteemp(:,1);
i            = find(t==t0):find(t==t1);
fteemp       = fteemp(i,2);

%
% Census Population by age
%
load CPS/popbyage.dat
[n,m]        = size(popbyage);
t            = popbyage(:,1);
i            = find(t==t0):find(t==t1);
grpop        = [NaN;popbyage(2:n,2)./popbyage(1:n-1,2)];
grpop        = grpop(i);
grpop16_64   = [NaN;sum(popbyage(2:n,5:8)')'./sum(popbyage(1:n-1,5:8)')'];
grpop16_64   = grpop16_64(i);
poptot       = popbyage(i,2);
popunder5    = popbyage(i,3);
pop5_15      = popbyage(i,4);
pop16_19     = popbyage(i,5);
pop20_24     = popbyage(i,6);
pop25_44     = popbyage(i,7);
pop45_64     = popbyage(i,8);
pop65plus    = popbyage(i,9);
empratio     = employed./(employed+pop65plus);
fteratio     = fteemp./(fteemp+pop65plus);
pop25_64     = pop25_44+pop45_64;
pop20_64     = pop20_24+pop25_64;
pop16_64     = pop16_19+pop20_64;
popratio     = pop20_64./(pop20_64+pop65plus);

%
% Measures relevant for model
%

irate        = .041;
tem          = stax./(pce+irate*kcd+cddep);
tauc         = tem./(1-tem);
staxd        = tauc.*pced;
impcd        = irate*kcd;
impgc        = irate*kgfix;
impc         = impcd+impgc;
contax       = tauc.*(pce-pced);
contax       = stax-staxd;
con          = pce-pced-contax+cddep+impc+govc-govcdef;
agdi         = gdi+cddep-stax+impc;
agdp         = gdp+cddep-stax+impc;
agnp         = gnp+cddep-stax+impc;
inve         = gpdi+pced-staxd+govi-govidef+netxm+row;
capinc       = prof+.3*pinc+rent+gprof+nint+ibt-stax+irate*(kgfix+kcd)+ ...
               desc+row;
capincdep    = capinc+dep+cddep;


disp('=================================================================')
fprintf('  US DATA MEASURES TO COMPARE TO MODEL, AVGERAGED %g:%g\n',t0,t1)
disp('=================================================================')
disp(' ') 
disp('  TOTAL ADJUSTED INCOME:')
disp(' ') 
disp(sprintf('    Labor Income                                          %7.4f',mean(labinc./agnp)))
disp(sprintf('      Compensation of employees                           %7.4f',mean(comp./agnp)))
disp(sprintf('        Wage and salary accruals                          %7.4f',mean(wage./agnp)))
disp(sprintf('        Supplements to W&S                                %7.4f',mean(supp./agnp)))
disp(sprintf('      70%% Proprietors'' income with IVA,CCAdj              %7.4f',mean(.7*pinc./agnp)))
disp(sprintf('    Capital Income                                        %7.4f',mean(capincdep./agnp)))
disp(sprintf('      Corporate profits                                   %7.4f',mean(prof./agnp)))
disp(sprintf('      30%% Proprietors'' income with IVA,CCAdj              %7.4f',mean(.3*pinc./agnp)))
disp(sprintf('      Rental income of persons with CCadj                 %7.4f',mean(rent./agnp)))
disp(sprintf('      Surplus of govt. enterprises                        %7.4f',mean(gprof./agnp)))
disp(sprintf('      Net interest and misc payments, domestic ind.       %7.4f',mean(nint./agnp)))
disp(sprintf('      Indirect business tax                               %7.4f',mean(ibt./agnp)))
disp(sprintf('        Taxes on production and imports                   %7.4f',mean(txpi./agnp)))
disp(sprintf('        Business current transfer payments                %7.4f',mean(bustra./agnp)))
disp(sprintf('        Less: Subsidies                                   %7.4f',mean(subs./agnp)))
disp(sprintf('      Less: sales tax                                     %7.4f',mean(stax./agnp)))
disp(sprintf('        Federal excise taxes                              %7.4f',mean(fedexcise./agnp)))
disp(sprintf('        Federal customs duties                            %7.4f',mean(fedcustom./agnp)))
disp(sprintf('        S&L sales taxes                                   %7.4f',mean(slsales./agnp)))
disp(sprintf('        Motor vehicle licenses                            %7.4f',mean(slmvl./agnp)))
disp(sprintf('        Severance taxes                                   %7.4f',mean(slsev./agnp)))
disp(sprintf('        Special assessments                               %7.4f',mean(slspec./agnp)))
disp(sprintf('        Other taxes on production and imports             %7.4f',mean(txpioth./agnp)))
disp(sprintf('      Income from ROW                                     %7.4f',mean(row./agnp)))
disp(sprintf('        Receipts                                          %7.4f',mean(nfir./agnp)))
disp(sprintf('        Payments                                          %7.4f',mean(nfip./agnp)))
disp(sprintf('      Tangible depreciation                               %7.4f',mean((dep+cddep)./agnp)))
disp(sprintf('        Consumption of fixed capital (NIPA)               %7.4f',mean(dep./agnp)))
disp(sprintf('        Consumer durable depreciation (FOF)               %7.4f',mean(cddep./agnp)))
disp(sprintf('      Imputed capital services                            %7.4f',mean(impc./agnp)))
disp(sprintf('        Consumer durables services                        %7.4f',mean(impcd./agnp)))
disp(sprintf('        Government capital services                       %7.4f',mean(impgc./agnp)))
disp(sprintf('      Statistical descrepancy                             %7.4f',mean(desc./agnp)))
disp(        '                                                          =======')
disp(sprintf('    Adjusted GNP                                          %7.4f',1.))
disp(sprintf('      NIPA GDP                                            %7.4f',mean(gdp./agnp)))
disp(sprintf('      ROW income                                          %7.4f',mean(row./agnp)))
disp(sprintf('      Adjustments to GDP                                  %7.4f',mean((agdp-gdp)./agnp)))
disp(' ')
disp('  TOTAL ADJUSTED PRODUCT:')
disp(' ')
disp(sprintf('    Consumption                                           %7.4f',1-mean((govdef+inve)./agnp)))
disp(sprintf('      Personal consumption expenditures                   %7.4f',mean(pce./agnp)))
disp(sprintf('      Less: Consumer durable goods                        %7.4f',mean(pced./agnp)))
disp(sprintf('      Less: Imputed sales tax, nondurables and services   %7.4f',mean(contax./agnp)))
disp(sprintf('      Plus: Imputed capital services, durables            %7.4f',mean(impcd./agnp)))
disp(sprintf('      Govt consumption expenditures, nondefense           %7.4f',mean((govc-govcdef)./agnp)))
disp(sprintf('      Plus: Imputed capital services, govt capital        %7.4f',mean(impgc./agnp)))
disp(sprintf('      Consumer durable depreciation (FOF)                 %7.4f',mean(cddep./agnp)))
disp(sprintf('    Defense spending                                      %7.4f',mean(govdef./agnp)))
disp(sprintf('      Govt consumption expenditures, natl defense         %7.4f',mean(govcdef./agnp)))
disp(sprintf('      Govt gross investment, natl defense                 %7.4f',mean(govidef./agnp)))
disp(sprintf('    Tangible investment                                   %7.4f',mean(inve./agnp)))
disp(sprintf('      Gross private domestic investment                   %7.4f',mean(gpdi./agnp)))
disp(sprintf('        Schedule C corporation                            %7.4f',mean(kc_cshare.*cgpdi./agnp)))
disp(sprintf('        Other private business                            %7.4f',mean((gpdi-kc_cshare.*cgpdi)./agnp)))
disp(sprintf('      Consumer durable goods                              %7.4f',mean(pced./agnp)))
disp(sprintf('      Less: Imputed sales tax, durables                   %7.4f',mean(staxd./agnp)))
disp(sprintf('      Govt investment expenditures, nondefense            %7.4f',mean((govi-govidef)./agnp)))
disp(sprintf('      Net exports of goods and services                   %7.4f',mean(netxm./agnp)))
disp(sprintf('      Income from ROW                                     %7.4f',mean(row./agnp)))
disp(        '                                                          =======')
disp(sprintf('    GNP                                                   %7.4f',1.))
disp(' ')
disp(' ')
disp('  CAPITAL STOCKS (End of period):')
disp(' ')
disp(sprintf('    Tangible capital                                      %7.4f',mean(k./agnp)))
disp(sprintf('      Fixed assets, private                               %7.4f',mean((kfix-kgfix)./agnp)))
disp(sprintf('        Schedule C corporations                           %7.4f',mean(kc_cshare.*kcfix./agnp)))
disp(sprintf('        Other private business                            %7.4f',mean((kfix-kgfix-kc_cshare.*kcfix)./agnp)))
disp(sprintf('          Residential                                     %7.4f',mean(kfixres./agnp)))
disp(sprintf('      Fixed assets, public                                %7.4f',mean(kgfix./agnp)))
disp(sprintf('      Consumer durables                                   %7.4f',mean(kcd./agnp)))
disp(sprintf('      Inventories                                         %7.4f',mean(kinv./agnp)))
disp(sprintf('        Schedule C corporations                           %7.4f',mean(kc_cshare.*kcinv./agnp)))
disp(sprintf('        Other private business                            %7.4f',mean((kinv-kc_cshare.*kcinv)./agnp)))
disp(sprintf('      Land                                                %7.4f',mean(land./agnp)))
disp(sprintf('        Schedule C corporations                           %7.4f',mean(kc_cshare.*nfcland./agnp)))
disp(sprintf('        Other private business                            %7.4f',mean((land-kc_cshare.*nfcland)./agnp)))
disp(sprintf('          Nonfinancial corporate                          %7.4f',mean(kc_sshare.*nfcland./agnp)))
disp(sprintf('          Nonfinancial noncorporate                       %7.4f',mean(nfncland./agnp)))
disp(sprintf('          Households and nonprofits                       %7.4f',mean(hhland./agnp)))
disp(' ')
disp(sprintf('    Intangible capital                                    %7.4f',1.718))
disp(sprintf('      Plant-specific (McGrattan-Prescott 2010)            %7.4f',1.198))
disp(sprintf('      Technology capital (McGrattan-Prescott 2010)        %7.4f',.519))
disp(' ')
disp(sprintf('    Total                                                 %7.4f',mean(k./agnp)+1.718))
disp(' ')
disp(' ')
disp('  GOVERNMENT EXPENDITURES, RECEIPTS, AND DEFICIT:')
disp(' ')
disp(sprintf('    Govt consumption expenditures and gross investment    %7.4f',mean(govtotal./agnp)))
disp(sprintf('      National defense                                    %7.4f',mean(natdefense./agnp)))
disp(sprintf('      Nondefense spending                                 %7.4f',mean((govtotal-natdefense)./agnp)))
disp(sprintf('        General public service                            %7.4f',mean(genpubser./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(genpubserc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(genpubseri./agnp)))
disp(sprintf('        Public order and safety                           %7.4f',mean(puborder./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(puborderc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(puborderi./agnp)))
disp(sprintf('        Economic affairs                                  %7.4f',mean(econaffairs./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(econaffairsc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(econaffairsi./agnp)))
disp(sprintf('        Housing and community services                    %7.4f',mean(housing./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(housingc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(housingi./agnp)))
disp(sprintf('        Health                                            %7.4f',mean(health./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(healthc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(healthi./agnp)))
disp(sprintf('        Recreation and culture                            %7.4f',mean(recreation./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(recreationc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(recreationi./agnp)))
disp(sprintf('        Education                                         %7.4f',mean(education./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(educationc./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(educationi./agnp)))
disp(sprintf('        Income security                                   %7.4f',mean(welfare./agnp)))
disp(sprintf('          consumption                                     %7.4f',mean(welfarec./agnp)))
disp(sprintf('          investment                                      %7.4f',mean(welfarei./agnp)))
disp(sprintf('    Govt transfers                                        %7.4f',mean(govetr./agnp)))
disp(sprintf('      Current transfers                                   %7.4f',mean(govetrcur./agnp)))
disp(sprintf('        Social security                                   %7.4f',mean(govetrsbss./agnp)))
disp(sprintf('        Medicare                                          %7.4f',mean(govetrsbmc./agnp)))
disp(sprintf('        Other                                             %7.4f',mean(govetroth./agnp)))
disp(sprintf('      Subsidies                                           %7.4f',mean(govetrsub./agnp)))
disp(sprintf('      Capital transfers                                   %7.4f',mean(govetrcap./agnp)))
disp(sprintf('    Interest payments                                     %7.4f',mean(govein./agnp)))
disp(sprintf('    Net purchases of nonproductive assets                 %7.4f',mean(govenp./agnp)))
disp(sprintf('    Less: Consumption of fixed capital                    %7.4f',mean(govdep./agnp)))
disp(        '                                                          =======')
disp(sprintf('    Total expenditures                                    %7.4f',mean(gove./agnp)))
disp(' ')               
disp(sprintf('    Taxes                                                 %7.4f',mean(govtax./agnp)))
disp(sprintf('      Personal current taxes                              %7.4f',mean(govtaxpc./agnp)))
disp(sprintf('      Taxes on production and imports                     %7.4f',mean(govtaxibt./agnp)))
disp(sprintf('      Taxes on corporate income                           %7.4f',mean(govtaxci./agnp)))
disp(sprintf('      Taxes from ROW                                      %7.4f',mean(govtaxrow./agnp)))
disp(sprintf('    Contributions for govt social insurance               %7.4f',mean(govrss./agnp)))
disp(sprintf('    Other                                                 %7.4f',mean((govrtr+govrin+govrge)./agnp)))
disp(sprintf('      Income receipts from assets                         %7.4f',mean(govrin./agnp)))
disp(sprintf('      Current transfer receipts                           %7.4f',mean(govrtrcur./agnp)))
disp(sprintf('      Capital transfer receipts                           %7.4f',mean(govrtrcap./agnp)))
disp(sprintf('      Surplus of govt. enterprises                        %7.4f',mean(govrge./agnp)))
disp(        '                                                          =======')
disp(sprintf('    Total receipts                                        %7.4f',mean(govr./agnp)))
disp(' ')
disp(sprintf('    Deficit                                               %7.4f',mean(deficit./agnp)))
disp(' ')
disp(' ')
disp(' BALANCE SHEET ITEMS:')
disp(' ')
disp(sprintf('    Household Net Worth:                                  %7.4f',mean(nwhh./agnp)))
disp(sprintf('      Assets                                              %7.4f',mean(vhha./agnp)))
disp(sprintf('        Tangible                                          %7.4f',mean(vhhtang./agnp)))
disp(sprintf('        Financial                                         %7.4f',mean(vhhfa./agnp)))
disp(sprintf('      Liabilities                                         %7.4f',mean(vhhli./agnp)))
disp(' ')
disp(sprintf('    Government debt:                                      %7.4f',mean(govsec./agnp)))
disp(sprintf('      State and local municipal securities                %7.4f',mean(slmuni./agnp)))
disp(sprintf('      Federal treasury securities                         %7.4f',mean(fedtrs./agnp)))
disp(sprintf('      Federal budget agency securities                    %7.4f',mean(fedagcy./agnp)))
disp(' ') 
disp(' ')
disp(' POPULATION, HOURS, AND EMPLOYMENT:')
disp(' ')
disp(sprintf('    Population, all ages (mil)                            %7.2f',mean(poptot/1000)))
disp(sprintf('    Population, 16-64 (mil)                               %7.2f',mean(pop16_64/1000)))
disp(sprintf('    Average growth in population, 16-64 (%%)               %7.4f',mean((grpop16_64-1)*100)))
disp(sprintf('    Average growth in population (%%)                      %7.4f',mean((grpop-1)*100)))
disp(sprintf('    Average growth in population, 16-64 (%%)               %7.4f',mean((grpop16_64-1)*100)))
disp(sprintf('    Average hours worked (rel. to 5200)                   %7.4f',mean(avghrsa)))
disp(sprintf('    # employed/(# employed+population 65+)                %7.4f',mean(empratio)))
disp(sprintf('    # FTEs/(# FTEs+population 65+)                        %7.4f',mean(fteratio)))
disp(sprintf('    Popuation 20-64/Population >=20                       %7.4f',mean(popratio)))
disp(' ')
disp(' ')
disp(' TAX RATES ON WAGE-LIKE INCOMES (Barro-Redlick)')
disp(' ') 
disp(sprintf('    Federal                                               %7.2f',mean(taulfed)*100))
disp(sprintf('    State and Local                                       %7.2f',mean(taulsl)*100))
disp(sprintf('    FICA                                                  %7.2f',mean(taulfica)*100))
disp(' ')
disp(sprintf('    Average Marginal (fed+sl+fica)                        %7.2f',mean(taul)*100))
disp(' ')
disp(' TAX RATE ON CORPORATE DISTRIBUTIONS ')
disp(' ') 
mqualdiv   = avgnan(qualdivrat);
mequuntax  = avgnan(equntaxedd);
mdivuntax  = avgnan(divuntaxed);
taud       = (1-mequuntax)*(mqualdiv*.2353+(1-mqualdiv)*.3198);
disp(        '    TAXSIM tax rates (only 2013): ')
disp(sprintf('      Ordinary dividends                                  %7.2f',31.98));
disp(sprintf('      Qualified dividends                                 %7.2f',23.53));

disp(sprintf('    Fraction of dividends qualified (after 2002)          %7.2f',mqualdiv*100));

disp(sprintf('    Domestic equity in Pensions/LI/NP (Source: FOF,ICI)   %7.2f',mequuntax*100));
disp(sprintf('    NIPA dividends to Pensions/LI/NP (Source: NIPA,IRS)   %7.2f',mdivuntax*100));
disp(' ')
disp(sprintf('    Estimate of tax rate on dividends                     %7.2f',taud*100));
disp('=================================================================')
disp(   ' NOTE: NIPA and FOF values all relative to adjusted GNP')


return
disp(' ')
disp(' ')
disp('====================================================')
disp('Check constistency (values relative to adjusted GNP)') 
disp('====================================================')
disp('   GDP')
disp(sprintf('     * NIPA                                 %7.4f',mean(gdp./agnp)))
disp(sprintf('     * FOF                                  %7.4f',mean(fofgdp./agnp)))
disp(' ')
disp('   PCE, Total')
disp(sprintf('     * NIPA                                 %7.4f',mean(pce./agnp)))
disp(sprintf('     * FOF                                  %7.4f',mean(fofpce./agnp)))
disp('     o Durables')
disp(sprintf('       * NIPA                               %7.4f',mean(pced./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofpced./agnp)))
disp(sprintf('       * FA                                 %7.4f',mean(xcd./agnp)))
disp('     o Nondurables')
disp(sprintf('       * NIPA                               %7.4f',mean(pcend./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofpcend./agnp)))
disp('     o Services')
disp(sprintf('       * NIPA                               %7.4f',mean(pces./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofpces./agnp)))
disp(' ')
disp('   GPDI, Total')
disp(sprintf('     * NIPA                                 %7.4f',mean(gpdi./agnp)))
disp(sprintf('     * FOF                                  %7.4f',mean(fofgpdi./agnp)))
disp('     o Fixed investment')
disp(sprintf('       * NIPA                               %7.4f',mean(gpdif./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofgpdif./agnp)))
disp(sprintf('       * FA                                 %7.4f',mean(xpfix./agnp)))
disp('       - Corporate ')
disp(sprintf('         * FOF                              %7.4f',mean(cgpdif./agnp)))
disp(sprintf('         * FA                               %7.4f',mean(xcfix./agnp)))
disp('       - Noncorporate ')
disp(sprintf('         * FOF                              %7.4f',mean(ncgpdif./agnp)))
disp(sprintf('         * FA                               %7.4f',mean(xncfix./agnp)))
disp('     o Change in private inventories')
disp(sprintf('       * NIPA                               %7.4f',mean(gpdii./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofgpdii./agnp)))
disp('       - Corporate ')
disp(sprintf('         * FOF                              %7.4f',mean(cgpdii./agnp)))
disp('       - Noncorporate ')
disp(sprintf('         * FOF                              %7.4f',mean(ncbgpdii./agnp)))
disp(' ')
disp('   Govt. spending, Total')
disp(sprintf('     * NIPA                                 %7.4f',mean((govi+govc)./agnp)))
disp(sprintf('     * FOF                                  %7.4f',mean((fofgovi+fofgovc)./agnp)))
disp('     o Consumption')
disp(sprintf('       * NIPA                               %7.4f',mean(govc./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofgovc./agnp)))
disp('     o Investment')
disp(sprintf('       * NIPA                               %7.4f',mean(govi./agnp)))
disp(sprintf('       * FOF                                %7.4f',mean(fofgovi./agnp)))
disp(sprintf('       * FA                                 %7.4f',mean(xgfix./agnp)))


t            = [t0:t1]';
lt           = length(t);
