% FIGURES      plot figures shown in 
%
%                 Productivity and the Post-1990 U.S. Economy
%                 Ellen McGrattan and Edward C. Prescott
%                 Federal Reserve Bank of St. Louis Review, 2005
%


%
% Load in data
%
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 prescott/prescott.csv
load FIXASSET/capital.dat
capital  = [capital(:,5:78),NaN*ones(34,1)];
capqind  = [capqind(:,5:78),NaN*ones(34,1)];
grcap    = capital(:,74)./capital(:,73);
capital(:,75) = capital(:,74) .*grcap;

%
% Assign names to raw data series
%

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,:)';
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,:)';
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)'];

%
% Sample:
%
t0       = find(t==1990);
t1       = find(t==2003);

%
% Interest and technological growth rate:
%
i        = .041;
gamma    = .02;

%
% Adjustments to national accounts
%
mu1      = cibt(t0:t1)./ibt(t0:t1);
mu2      = 1-mu1;
cons     = pce(t0:t1)./gdp(t0:t1) ...
           +i*kcd(t0:t1)./gdp(t0:t1) ...
           +cddep(t0:t1)./gdp(t0:t1);
tem      = (stax(t0:t1)./gdp(t0:t1))./cons;
tauc     = tem./(1-tem);
cd       = (pced(t0:t1)./gdp(t0:t1))./(1+tauc);
tauccd   = tauc.*cd;
taucc    = stax(t0:t1)./gdp(t0:t1)-tauccd;

p1y1mx1u = cgdi(t0:t1)./gdp(t0:t1) ...
           -mu1.*(stax(t0:t1)./gdp(t0:t1))  ...
           -(pbexp(t0:t1)./gdp(t0:t1));
del1mk1m = cdep(t0:t1)./gdp(t0:t1);
tau1kk1m = cibt(t0:t1)./gdp(t0:t1) ...
           -mu1.*(stax(t0:t1)./gdp(t0:t1));
wn1      = ccomp(t0:t1)./gdp(t0:t1);
adjprof  = prof(t0:t1)./gdp(t0:t1);
fedprof  = fedres(t0:t1)./gdp(t0:t1);
profbtx  = pbtd(t0:t1)./gdp(t0:t1);
proftax  = ctax(t0:t1)./gdp(t0:t1);
ccadjp   = ccadj(t0:t1)./gdp(t0:t1);
tau1hold = (proftax-fedprof)./(adjprof-ccadjp-fedprof);
tau1     = (proftax-fedprof)./(profbtx-fedprof);
p2y2     = (gdi(t0:t1)-cgdi(t0:t1))./gdp(t0:t1) ...
            -mu2.* (stax(t0:t1)./gdp(t0:t1))  ...
            +cddep(t0:t1)./gdp(t0:t1) ...
            +i*(kgfix(t0:t1)+kcd(t0:t1))./gdp(t0:t1) ...
            -0*(nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1) ...
            +desc(t0:t1)./gdp(t0:t1);
del2mk2m = (dep(t0:t1)-cdep(t0:t1)+cddep(t0:t1))./gdp(t0:t1);
tau2kk2m = (ibt(t0:t1)-cibt(t0:t1))./gdp(t0:t1) ...
           -mu2.*(stax(t0:t1)./gdp(t0:t1));
wn2      = (comp(t0:t1)+.7*pinc(t0:t1)-ccomp(t0:t1)+ ...
            desc(t0:t1))./gdp(t0:t1);
ncprof   = (rent(t0:t1)+.3*pinc(t0:t1)+gprof(t0:t1))./gdp(t0:t1) ...
           +i*(kgfix(t0:t1)+kcd(t0:t1))./gdp(t0:t1) ...
           +(nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1);
yincmx1u = p1y1mx1u+p2y2;
c        = (pce(t0:t1)-pced(t0:t1))./gdp(t0:t1) ...
           +i*kcd(t0:t1)./gdp(t0:t1) ...
           +cddep(t0:t1)./gdp(t0:t1) ...
           -taucc  ...
           -pbexp(t0:t1)./gdp(t0:t1) ...
           -0*(nint(t0:t1)-cnint(t0:t1))./gdp(t0:t1);
g        = govc(t0:t1)./gdp(t0:t1) ...
           +i*kgfix(t0:t1)./gdp(t0:t1);
x1m      = (cgpdif(t0:t1)+cgpdii(t0:t1))./gdp(t0:t1);
x2m      = (pced(t0:t1)+govi(t0:t1)+gpdi(t0:t1)+netxm(t0:t1)- ...
            cgpdif(t0:t1)-cgpdii(t0:t1))./gdp(t0:t1)-tauccd;

yprodmx1u= c+x1m+x2m+g;
   
k1mbop   = kc(t0-1:t1-1) ./gdp(t0:t1);
k1meop   = kc(t0:t1)     ./gdp(t0:t1);
k2mbop   = knc(t0-1:t1-1)./gdp(t0:t1);
k2meop   = knc(t0:t1)    ./gdp(t0:t1);
k1m      = .5*k1mbop+.5*k1meop;
k2m      = .5*k2mbop+.5*k2meop;
n        = hcps(t0:t1);
n2       = n./(1+wn1./wn2);
n1       = n-n2;
w        = wn2./n2;

%
% Time series:
%
gdp0     = gdp(t0);
pcps     = pcps(t0:t1)/pcps(t0);
dgdp     = dgdp(t0:t1)/dgdp(t0);
gdp      = gdp(t0:t1)/gdp(t0);
grate    = (1+gamma).^[0:t1-t0]';
ct       = c.*gdp       ./(dgdp.*grate.*pcps);
x1mt     = x1m.*gdp     ./(dgdp.*grate.*pcps);
k1mt     = k1m.*gdp     ./(dgdp.*grate.*pcps);
k2mt     = k2m.*gdp     ./(dgdp.*grate.*pcps);
x2mt     = x2m.*gdp     ./(dgdp.*grate.*pcps);
gt       = g.*gdp       ./(dgdp.*grate.*pcps);
wt       = w.*gdp       ./(dgdp.*grate.*pcps);
p2y2t    = p2y2.*gdp    ./(dgdp.*grate.*pcps);
p1y1mx1ut= p1y1mx1u.*gdp./(dgdp.*grate.*pcps);
gdpt     = gdp          ./(dgdp.*grate.*pcps);
tau1kkt  = tau1kk1m.*gdp./(dgdp.*grate.*pcps);
del1mkt  = del1mk1m.*gdp./(dgdp.*grate.*pcps);
wn1t     = wn1.*gdp     ./(dgdp.*grate.*pcps);
wn2t     = wn2.*gdp     ./(dgdp.*grate.*pcps);

%
% Capital shares:
%
alpha    = .496;
thetm    = .2825;
thetu    = .0473;
theta    = thetm+thetu;

%
% Two measures of intangible:
%
x1uta    = 1/(1-theta)*(wn1t)-p1y1mx1ut;
x1utb    = (1-alpha)/(1-theta)*(p2y2t.*wn1t./wn2t)-p1y1mx1ut;

%
% Implied economic outputs
%
yta      = ct+x1mt+x1uta+x2mt+gt;
ytb      = ct+x1mt+x1utb+x2mt+gt;
p1y1ta   = p1y1mx1ut+x1uta;
p1y1tb   = p1y1mx1ut+x1utb;

%
% Figures
%

figure(1)
plot(t(t0:t1),x1uta)
title('Intangible Investment (Measure 1)','Fontsize',14)
axis([1990,2003,0,.04])
set(gca,'Fontsize',14)

figure(2)
plot(t(t0:t1),[gdpt./n*n(1),yta./n*n(1)])
title('Productivity Relative to 2% trend (Measure 1)','Fontsize',14)
axis([1990,2003,.95,1.1])
set(gca,'Fontsize',14)
legend('NIPA','Adjusted',2)

figure(3)
plot(t(t0:t1),[(p1y1mx1ut-wn1t-del1mkt-tau1kkt)./gdpt, ...
               (p1y1ta-wn1t-del1mkt)./yta])
title('Corporate Profit Share (Measure 1)','Fontsize',14)
axis([1990,2003,.06,.13])
set(gca,'Fontsize',14)
legend('NIPA','Adjusted',2)

figure(4)
plot(t(t0:t1),[x1mt./gdpt,(x1mt+x1uta)./yta])
title('Corporate Investment Share (Measure 1)','Fontsize',14)
axis([1990,2003,.08,.14])
set(gca,'Fontsize',14)
legend('NIPA','Adjusted',2)

figure(5)
plot(t(t0:t1),x1utb)
title('Intangible Investment (Measure 2)','Fontsize',14)
axis([1990,2003,0,.06])
set(gca,'Fontsize',14)

figure(6)
plot(t(t0:t1),[gdpt./n*n(1),ytb./n*n(1)])
title('Productivity Relative to 2% trend (Measure 2)','Fontsize',14)
axis([1990,2003,.95,1.1])
set(gca,'Fontsize',14)
legend('NIPA','Adjusted',2)
