%WEDGES  Given data on y,k,l,c,i, compute a,g,taun,taux as follows:
% 
%            a[t]      = y[t]/(k[t-1]^theta * l[t]^(1-theta))
%            g[t]      = y[t]-i[t]-c[t] 
%            1-taun[t] = [(1+tauc)U_2(c[t],1-l[t])/U_1(c[t],1-l[t])]/
%                        [a[t] * (1-theta)*(k[t-1]/l[t])^theta ]
%            1+taux[t] = {(1+taux[t-1])*U_1[t-1]/(betg * U_1[t])
%                          - (1-tauk)* a[t]* theta (l[t]/k[t-1])^(1-theta)
%                          - tauk * delta } / (1-delta)
%
%            where a[t] = z[t]^(1-theta) 
%
%        See aea.tex for details.
%

%        Ellen McGrattan, 12-4-01
%        Revised, ERM, 4-23-02

t0  = 1929;
a1  = 0;
a2  = 1;
a3  = 0;
a4  = 1;
a5  = 0;
%
% Assumptions:
%     A1.  Durables 
%          a1 = 0   if included with investment
%             = 1   if included with private consumption
%     A2.  Net foreign investment
%          a2 = 0   if included with investment
%             = 1   if excluded from accounts
%     A3.  Government investment
%          a3 = 0   if military included with government consumption, 
%                   the remainder in investment
%             = 1   if all included with government consumption 
%             = 2   if all included with investment
%     A4.  Owner-occupied housing 
%          a4 = 0   if included in accounts
%             = 1   if excluded from accounts
%
% NOTE: Automatic adjustments made:
%
%      1.  Sales tax removed from IBT and from PCE
%      2.  Returns imputed to government capital at 4% rather than 0%

%
% Load in data, 1929-2000, unless noted:
%     * nominal National Accounts data 
%     * nominal capital stocks and investments
%     * deflators
%     * hours
%     * population over 16
%

t    = [1929:2000]';
t1   = find(t==2000);
tby  = find(t==1996);
load nipa/101
load nipa/102
load nipa/109
load nipa/110
load nipa/114
load nipa/305
load nipa/307
load nipa/309
load nipa/501
load nipa/511a
load nipa/605a
load nipa/605b
load nipa/605c
load nipa/608a
load nipa/608b
load nipa/608c
load nipa/609b
load nipa/609c
load nipa/701
load nipa/719
load nipa/821
load nipa/fa101
load nipa/fa102
load nipa/fa103
load nipa/fa104
load nipa/fa105
load nipa/fa106
load nipa/fa501
load nipa/fa502
load nipa/fa504
load nipa/fa505
load nipa/fa507
load nipa/fa508
load nipa/fa701
load nipa/fa702
load nipa/fa703
load nipa/fa704
load nipa/fa705
load nipa/fa706
load nipa/def82.dat
load labor/pop16.dat
load labor/mhrs_us.dat
X605a  = X605a*1000;
X605b  = X605b*1000;
X605c  = X605c*1000;
X608a  = X608a*1000;
X608b  = X608b*1000;
X608c  = X608c*1000;
X609b  = X609b*10^6 *1.1;    %1.1 factor to scale to Kendrick's estimates
X609c  = X609c*10^6 *1.1;
X605   = [X605a([1:4,80],1:19),X605b([1:4,81],1:39),X605c([1:4,81],:)]*10^3;
X608   = [X608a([1:4,80],1:19),X608b([1:4,81],1:39),X608c([1:4,81],:)]*10^3;
X609   = [NaN*ones(4,19),      X609b([1:4],1:39),   X609c([1:4],:);
          NaN*ones(1,19), 1700*X605b(81,1:39), 1700*X605c(81,:)];
khrs   = mhrs_us(141:212,[1,2,2,7,5])*10^6;
pop    = pop16(30:101,2)*1000;
fpop   = [pop16(31:101,2);pop16(101,2)^2/pop16(100,2)]*1000;

Efa101 = fa101(:,5:76)/1000;
Bfa101 = fa101(:,4:75)/1000;
Dfa103 = fa103(:,5:76)/1000;
Ifa105 = fa105(:,29:100)/1000;
Efa102 = fa102(:,5:76).*(fa101(:,72)*ones(1,72))/100/1000;
Bfa102 = fa102(:,4:75).*(fa101(:,72)*ones(1,72))/100/1000;
Dfa104 = fa104(:,5:76).*(fa103(:,72)*ones(1,72))/100/1000;
Ifa106 = fa106(:,29:100).*(fa105(:,96)*ones(1,72))/100/1000;
Efa501 = fa501(:,5:76)/1000;
Bfa501 = fa501(:,4:75)/1000;
Dfa504 = fa504(:,5:76)/1000;
Ifa507 = fa507(:,29:100)/1000;
Efa502 = fa502(:,5:76).*(fa501(:,72)*ones(1,72))/100/1000;
Bfa502 = fa502(:,4:75).*(fa501(:,72)*ones(1,72))/100/1000;
Dfa505 = fa505(:,5:76).*(fa504(:,72)*ones(1,72))/100/1000;
Ifa508 = fa508(:,29:100).*(fa507(:,96)*ones(1,72))/100/1000;
Efa701 = fa701(:,5:76)/1000;
Bfa701 = fa701(:,4:75)/1000;
Dfa703 = fa703(:,5:76)/1000;
Ifa705 = fa705(:,16:87)/1000;
Efa702 = fa702(:,5:76).*(fa701(:,72)*ones(1,72))/100/1000;
Bfa702 = fa702(:,4:75).*(fa701(:,72)*ones(1,72))/100/1000;
Dfa704 = fa704(:,5:76).*(fa703(:,72)*ones(1,72))/100/1000;
Ifa706 = fa706(:,16:87).*(fa705(:,83)*ones(1,72))/100/1000;

%Product
ngnp   = X109(5,:)';   %=nc+ng+ngdi+nnfi = nc+ngc+ngdi+ngi+nnfi
nc     = X101(3,:)';
ncd    = X101(4,:)';
ncnd   = X101(5,:)';
ncs    = X101(6,:)';
nch    = X821(113,:)'-X821(116,:)'-X821(150,:)';
ngdi   = X101(7,:)';
nif    = X101(8,:)';
nifh   = X821(150,:)';
nii    = X101(13,:)';
ng     = X101(21,:)';
ngc    = X309(2,:)';
ngi    = X501(21,:)';
ngim   = X307(12,:)';
nnfi   = ngnp-nc-ngdi-ng;

%Income
ngni   = X109(5,:)';
ndep   = X109(6,:)';
nibt   = X109(14,:)'+X109(15,:)'+X109(16,:)'-X109(17,:)';
nstx   = X305(4,:)'+X305(18,:)'+X305(34,:)';
nwl    = X114(3,:)'+.8*X114(10,:)';   
nrk    = X114(2,:)'-nwl; 
ndeph  = X821(126,:)';
nibth  = X821(122,:)'-X821(125,:)';
nwlh   = .8*X821(129,:)';
nrkh   = .2*X821(129,:)'+X821(130,:)'+X821(119,:)';
Tauc   = nstx./(nc-nstx);

%Capital/Investment
nlk    = Bfa101(2,:)';
nk     = Efa101(2,:)';
nlkif  = Bfa101(4,:)';
nkif   = Efa101(4,:)';
nlkifh = Bfa501(15,:)';
nkifh  = Efa501(15,:)';
nlkgi  = Bfa101(9,:)';
nkgi   = Efa101(9,:)';
nlkgim = Bfa701(19,:)'+Bfa701(30,:)';
nkgim  = Efa701(19,:)'+Efa701(30,:)';
nlkcd  = Bfa101(14,:)';
nkcd   = Efa101(14,:)';
niif   = Ifa105(4,:)';
niifh  = Ifa507(15,:)';
nigi   = Ifa105(9,:)';
nicd   = Ifa105(14,:)';
ndgi   = Dfa103(9,:)';
ndcd   = Dfa103(14,:)';
ndifh  = Dfa504(15,:)';
ndgim  = Dfa703(19,:)'+Dfa703(30,:)';
nkfor(1)  = nnfi(1);  % estimated foreign capital
nlkfor(1) = 0;
for j=1:t1-1 
 nkfor(j+1,1) = (1-.05)*nkfor(j)+nnfi(j+1);
 nlkfor(j+1,1) = nkfor(j);
end;
nlk    = nlk+nlkfor;
nk     = nk+nkfor;

%Population/Hours
pop    = pop/10^9;  % put in billions
fpop   = fpop/10^9;
hrsa   = X609(2,:)'/10^9;
hrsd   = X609(3,:)'/10^9;
hrsp   = X609(4,:)'/10^9;
hrsm   = X609(5,:)'/10^9;
empa   = X605(2,:)'/10^9;
empd   = X605(3,:)'/10^9;
empp   = X605(4,:)'/10^9;
empm   = X605(5,:)'/10^9;
enga   = X608(2,:)'/10^9;
engd   = X608(3,:)'/10^9;
engp   = X608(4,:)'/10^9;
engm   = X608(5,:)'/10^9;
mhrsa  = hrsa+hrsa .*(enga-empa)./empa;
mhrsd  = hrsd+hrsd .*(engd-empd)./empd;
mhrsp  = hrsp+hrsp .*(engp-empp)./empp;
mhrsm  = hrsm+hrsm .*(engm-empm)./empm;
khrsa  = khrs(:,2)/10^9;
khrsd  = khrs(:,3)/10^9;
khrsp  = khrs(:,4)/10^9;
khrsm  = khrs(:,5)/10^9;
mhrsa(1:19) = khrsa(1:19);
mhrsd(1:19) = khrsd(1:19);
mhrsp(1:19) = khrsp(1:19);
mhrsm(1:19) = khrsm(1:19);

%Automatic adjustments
ik     = .04;
coef   = 1;
ngnp   = ngnp-coef*nstx+ik*nkgi;
ngni   = ngni-coef*nstx+ik*nkgi;
ncd    = ncd-coef*nstx.*ncd./nc;
ncnd   = ncnd-coef*nstx.*ncnd./nc;
ncs    = ncs-coef*nstx.*ncs./nc;
nc     = nc-coef*nstx;
ng     = ng+ik*nkgi;
ngc    = ngc+ik*nkgi;
nibt   = nibt-coef*nstx;
nrk    = nrk+ik*nkgi;

%Requested adjustments
if a1==0;     % if durables included in investment
  ngnp  = ngnp+ik*nkcd+ coef*ndcd;
  ngni  = ngni+ik*nkcd+ coef*ndcd;
  nc    = nc+ik*nkcd+ coef*ndcd-ncd;
  ncs   = ncs+ik*nkcd+ coef*ndcd;
  nif   = nif+ncd;
  ngdi  = ngdi+ncd;
  ncd   = ncd-ncd;
  ndep  = ndep+ coef*ndcd;
  nlkif = nlkif+nlkcd;
  nkif  = nkif +nkcd;
  niif  = niif+nicd;
  nicd  = nicd-nicd;
  nlkcd = nlkcd-nlkcd;
  nkcd  = nkcd -nkcd;
end;
if a1==1;     % if durables included in consumption
  nlk   = nlk-nlkcd;
  nlkcd = nlkcd-nlkcd;
  nk    = nk-nkcd;
  nkcd  = nkcd-nkcd;
end;
if a2==1;    % if net foreign investment is excluded
  ngnp  = ngnp-nnfi;
  ngni  = ngni-nnfi;
  nrk   = nrk-nnfi;
  nnfi  = nnfi-nnfi; 
  nlk   = nlk-nlkfor;
  nk    = nk-nkfor;
end 
if a3==0;    % if military investment is put in govt consumption
  ngnp  = ngnp- coef*ndgim;
  ngni  = ngni- coef*ndgim;
  ng    = ng- coef*ndgim;
  ngc   = ngc+ngim- coef*ndgim;
  ngi   = ngi-ngim;
  ndep  = ndep- coef*ndgim;
  nlk   = nlk-nlkgim;
  nk    = nk-nkgim;
  nlkgi = nlkgi-nlkgim;
  nkgi  = nkgi-nkgim;
  nlkgim= nlkgim-nlkgim;
  nkgim = nkgim-nkgim;
end
if a3==1;    % if all govt investment is put in govt consumption
  ngnp  = ngnp-ndgi;
  ngni  = ngni-ndgi;
  ngc   = ngc+ngi-ndgi;
  ngi   = ngi-ngi;
  ndep  = ndep-ndgi;
  nlk   = nlk-nlkgi;
  nk    = nk-nkgi;
  nlkgi = nlkgi-nlkgi;
  nkgi  = nkgi-nkgi;
  nlkgim= nlkgim-nlkgim;
  nlkgim= nkgim-nkgim;
end;
if a4==1;    % if owner-occupied housing is excluded
  ngnp  = ngnp-nch-nifh;
  ngni  = ngni-ndeph-nibth-nrkh-nwlh;
  nc    = nc-nch;
  ncs   = ncs-nch;
  ngdi  = ngdi-nifh;
  nif   = nif-nifh;
  ndep  = ndep-ndeph;
  nibt  = nibt-nibth;
  nrk   = nrk-nrkh;
  nwl   = nwl-nwlh;
  nlk   = nlk-nlkifh;
  nk    = nk-nkifh;
  nlkif = nlkif-nlkifh;
  nkif  = nkif-nkifh;
  nlkifh= nlkifh-nlkifh;
  nkifh = nkifh-nkifh;
  niif  = niif-niifh;
  niifh = niifh-niifh;
end;
pgnp    = X109(5,:)'./X110(5,:)';
pgnp    = pgnp/pgnp(tby);
lpgnp   = [pgnp(1);pgnp(1:t1-1)];

%Real inventories:
iinv     = [X511a(2,:)';X102(13,(t1-2):t1)'];
kinv(t1) = 1505;
for j=t1-1:-1:1;
  kinv(j)    = kinv(j+1)-iinv(j+1);
  lkinv(j+1) = kinv(j);
end;
lkinv(1) = lkinv(2);
kinv     = kinv(:);
lkinv    = lkinv(:);


rgnp     = ngnp./pgnp;
rc       = (ncd+ncnd+ncs)./pgnp;
ri       = (ngdi+ngi+nnfi)./pgnp;
rg       = ngc./pgnp;
rlk      = nlk./lpgnp+lkinv;
rk       = nk./pgnp+kinv;
l        = (mhrsa-mhrsm)./pop/(50*100);


%
%Real, per-capita, detrended data:
%
gn       = .015;   % population growth in early 20th century
gz       = .016;   % technology growth in early 20th century
grate    = pop.*(1+gz).^(t-t0);
fgrate   = fpop.*(1+gz).^(t+1-t0);
k        = rk./fgrate;
lk       = rlk./grate;
i        = ri./grate;
c        = rc./grate;
g        = rg./grate;
gnp      = rgnp./grate;

clag     = [NaN;c(1:t1-1)];
llag     = [NaN;l(1:t1-1)];

%
%Parameters
%
util  = 1;          %util = 1:  u(c,1-l) = (c*(1-l)^psi)^phi/phi;
                    %     = 2:  u(c,1-l) = log(c) + psi*((1-l)^phi-1)/phi
phi   = .000001;
tauks = .3;
tauns = .035;
tauxs = .0;
ir    = .045;
beta  = ( (1+gz)^(1-phi) )/(1+ir);
betg  = beta*(1+gz)^(phi-1);
delta = i(t0-1928)/lk(t0-1928)+1-(1+gz)*(1+gn);
theta = ((1/betg-1)/(1-tauks)+delta)*lk(t0-1928)/gnp(t0-1928); 
tauc  = mean(Tauc(t0-1928:t0-1928+5));
psi   = (1/l(t0-1928)-1)/(c(t0-1928)/gnp(t0-1928)*(1+tauc)/ ...
        (1-theta)/(1-tauns) );

%
% Using real, per-capita, detrended data, compute wedges:
%

if util==1;
  u     = (c.*(1-l).^psi).^phi/phi;
  ulag  = (clag.*(1-llag).^psi).^phi/phi;
  u1    = phi*u./c;
  u1lag = phi*ulag./clag;
  u2    = psi*phi*u./(1-l);
else
  u     = log(c)+psi*((1-l).^phi-1)/phi;
  ulag  = log(clag)+psi*((1-llag).^phi-1)/phi;
  u1    = ones(length(c),1)./c;
  u1lag = ones(length(c),1)./clag;
  u2    = psi*(1-l).^(phi-1);
end

a     = gnp./(lk.^theta .* l.^(1-theta));
taun  = 1-((1+tauc)*u2./u1)./( gnp./l *(1-theta));
taux  = tauxs*ones(1,t1);
tauxl = tauxs*ones(1,t1);
for j=(t0-1927):t1;
  taux(j) = ( (1+taux(j-1))*u1lag(j)/(betg*u1(j)) ....
                  -(1-tauks) *gnp(j)/lk(j)*theta  ...
                  -tauks*delta )/(1-delta) -1 ;
  tauxl(j)=taux(j-1);
end;
taux  = taux(:);
tauxl = tauxl(:);

gnp0   = gnp(t0-1928);
figure(1)
plot(t,gnp/gnp0,'k-',t,c/gnp0,'b-',t,i/gnp0,'r-',t,g/gnp0,'g-')
axis([t0,t0+10,0,1.2])

figure(2)
plot(t,k/gnp0)
axis([t0,t0+10,1,5])

figure(3)
plot(t,l,'k-')
axis([t0,t0+10,0,.4])

wedges = [t,g/gnp0,taux,taun,a.^(1/(1-theta))/gnp0];
aggs   = [t,gnp/gnp0,c/gnp0,i/gnp0,g/gnp0,lk/gnp0,k/gnp0,l];
disp('    Year    GNP        C        I        G      K(-1)      K        L')
disp(sprintf('%8g %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n', ...
     aggs(t0-1928:t0-1928+10,:)'))

disp('    Year    C_G      Tau_x    Tau_n    Z  ')          
disp(sprintf('%8g %8.4f %8.4f %8.4f %8.4f\n',wedges(t0-1928:t0-1928+10,:)'))

figure(4)
plot(t,gnp/gnp0,t,a/a(t0-1928))
title('Real Per-Capita Detrended GNP and A-wedge')
legend('GNP','A')
axis([t0,t0+10,.5,1.1])

figure(5)
plot(t,gnp/gnp0,t,(1-taun)/(1-taun(t0-1928)))
title('Real Per-Capita Detrended GNP and \tau_n-wedge ')
legend('GNP','1-\tau_n')
axis([t0,t0+10,.5,1.1])

figure(6)
plot(t,gnp/gnp0,t,(ones(t1,1)./(1+taux))*(1+taux(t0-1928)))
title('Real Per-Capita Detrended GNP and \tau_x-wedge ')
legend('GNP','1/(1+\tau_x)')
axis([t0,t0+10,.5,1.5])

figure(7)
plot(t,gnp/gnp0,t,a/a(t0-1928),t,(1-taun)/(1-taun(t0-1928)))
title('Real Per-Capita Detrended GNP, A, and \tau_n-wedge ')
legend('GNP','A','1-\tau_n')
axis([t0,t0+10,.5,1.1])

tem = [t,gnp/gnp0,c/gnp0,i/gnp0,g/gnp0,k/gnp0,l];
save nipa29_00.dat  tem -ascii
