%
% PLOTBG    plot balanced growth results from bg*.f90. For details,
%           see the appendix (../appendix/appendix.pdf).
%

%           Ellen McGrattan, 12-19-10
%           Revised, ERM, 9-20-13

load bg.dat

iplot     = 0;
alb       = bg(1);
aub       = bg(2);
alpha     = bg(3);
beta      = bg(4);
delt1     = bg(5);
deli1     = bg(6);
delt2     = bg(7);
deli2     = bg(8);
gam       = bg(9);
gn        = bg(10);
gz        = bg(11);
lam       = bg(12);
tauc      = bg(13);
taud      = bg(14);
taul      = bg(15);
taupc     = bg(16);
taupu     = bg(17);
tfp       = bg(18);
thett1    = bg(19);
theti1    = bg(20);
thett2    = bg(21);
theti2    = bg(22);
B         = bg(23);
G         = bg(24);
alph1     = 1-alpha;
KT1       = bg(25);
KT2       = bg(26);
KI1       = bg(27);
KI2       = bg(28);
L1        = bg(29);
L2        = bg(30);
C         = bg(31);
XT1       = bg(32);
XT2       = bg(33);
XI1       = bg(34);
XI2       = bg(35);
Y         = bg(36);
irate     = bg(37);
Tran      = bg(38);
J         = bg(39);
n         = bg(40);
thetl1    = 1-thett1-theti1;
thetl2    = 1-thett2-theti2;
rt1       = (irate-1)/(1-taupc)+delt1;
rt2       = (irate-1)/(1-taupu)+delt2;
ri1       = irate-1+deli1;
ri2       = irate-1+deli2;
w         = thetl1*alpha*Y/L1;
agrid     = linspace(alb,aub,n)';
eps       = bg(41:40+J);
ps        = bg(41+J:40+2*J);
xi        = bg(41+2*J:40+3*J);
jpi       = zeros(J,1);
cj        = zeros(J,1);
lj        = zeros(J,1);
aj        = zeros(J,1);
apj       = zeros(J,1);
opt       = zeros(J,n);
jpi(:)    = bg(41+3*J:40+4*J);
cj(:)     = bg(41+4*J:40+5*J);
lj(:)     = bg(41+5*J:40+6*J);
aj(:)     = bg(41+6*J:40+7*J);
apj(:)    = bg(41+7*J:40+8*J);
opt(:)    = bg(41+8*J:40+J*(n+8));
beqp      = bg(length(bg)-1);
jbeq      = bg(length(bg));
mu(1,1)   = 1/(1+sum(cumprod(ps/(1+gn))));
for j=2:J;
  mu(j,1) = mu(j-1)*ps(j-1)/(1.+gn);
end;
Jr        = find(eps<1e-10);
Jr        = Jr(1);
Tranj     = sum(xi.*mu);
Tranwrk   = sum((Tran+xi(1:Jr-1)).*mu(1:Jr-1));
Tranret   = sum((Tran+xi(Jr:J)).*mu(Jr:J));

Y1        = tfp*KT1^thett1*KI1^theti1*L1^thetl1;
Y2        = tfp*KT2^thett2*KI2^theti2*L2^thetl2;
Ychk      = 2*(Y1)^alpha*(Y2)^alph1;
p1        = alpha*Y/Y1;
p2        = alph1*Y/Y2;
cprof     = p1*Y1-w*L1-delt1*KT1-XI1;
uprof     = p2*Y2-w*L2-delt2*KT2-XI2;
D         = p1*Y1-w*L1-XT1-XI1-taupc*cprof;
KT1eoy    = (1+gz)*(1+gn)*KT1;
KT2eoy    = (1+gz)*(1+gn)*KT2;
KI1eoy    = (1+gz)*(1+gn)*KI1;
KI2eoy    = (1+gz)*(1+gn)*KI2;
Beoy      = (1+gz)*(1+gn)*B;
vT1       = (1-taud)*KT1eoy;
vI1       = (1-taud)*(1-taupc)*KI1eoy;
vT2       = KT2eoy;
vI2       = (1-taupu)*KI2eoy;
v         = vT1+vT2+vI1+vI2;

Rev       = taud*D+taupc*cprof+taupu*uprof+taul*w*(L1+L2)+tauc*C;
GDP       = Y-XI1-XI2;
Int       = (irate-1)*B;
Beq       = G+Tran+Tranj-Rev-Beoy+irate*B;
GovR      = Beq+Rev+max(-Int,0);
GovE      = G+Tran+Tranj+max(Int,0);
Def       = GovE-GovR;
NWhh      = v+Beoy;
NWgov     = KT1eoy+KT2eoy+KI1eoy+KI2eoy-v-Beoy;

Cj        = cj.*mu;
Lj        = lj.*mu;
Aj        = apj.*mu;
U         = sum((beta.^[0:J-1]').*cumprod([1;ps(1:J-1)]).* ...
            (log(cj)+gam*log(1-lj)));
Disc      = sum((beta.^[0:J-1]').*cumprod([1;ps(1:J-1)]));

age = 20+[1:J];


disp('Results of OLG Model with Intangible Capital')
disp('--------------------------------------------')
disp('  Parameters')
disp('  ----------')
disp(sprintf('    Asset grid lower bound               %8.3f ',alb))
disp(sprintf('    Asset grid upper bound               %8.3f ',aub))
disp(sprintf('    Corporate share of output            %8.3f ',alpha))
disp(sprintf('    Discount factor                      %8.3f ',beta))
disp(sprintf('    Tangible depreciation, sector 1      %8.3f ',delt1))
disp(sprintf('    Intangible depreciation, sector 1    %8.3f ',deli1))
disp(sprintf('    Tangible depreciation, sector 2      %8.3f ',delt2))
disp(sprintf('    Intangible depreciation, sector 2    %8.3f ',deli2))
disp(sprintf('    Weight on leisure in utility         %8.3f ',gam))
disp(sprintf('    Bequest parameter                    %8.3f ',beqp))
disp(sprintf('    Bequest at J/c                       %8.3f ',beqp/mean(cj)))
disp(sprintf('    Growth of population (%%)             %8.3f ',gn*100))
disp(sprintf('    Growth of technology (%%)             %8.3f ',gz*100))
disp(sprintf('    Tax rate on consumption              %8.3f ',tauc))
disp(sprintf('    Tax rate on distributions            %8.3f ',taud))
disp(sprintf('    Tax rate on labor                    %8.3f ',taul))
disp(sprintf('    Tax rate on corporate profits        %8.3f ',taupc))
disp(sprintf('    Tax rate on noncorp. profits         %8.3f ',taupu))
disp(sprintf('    Tangible capital share, sector 1     %8.3f ',thett1))
disp(sprintf('    Intangible capital share, sector 1   %8.3f ',theti1))
disp(sprintf('    Tangible capital share, sector 2     %8.3f ',thett2))
disp(sprintf('    Intangible capital share, sector 2   %8.3f ',theti2))
disp(sprintf('    TFP parameter                        %8.3f ',tfp))
disp(sprintf('    Fraction working                     %8.3f ',sum(mu.*eps)))
disp(sprintf('    Fraction retired                     %8.3f ',1-sum(mu.*eps)))
disp(sprintf('    Age of first job                     %8g ',21))
disp(sprintf('    Age when retired                     %8g ',Jr+20))
disp(sprintf('    Maximum age                          %8g ',J+20))
disp(sprintf('    Bequest age                          %8g ',jbeq+20))
disp(sprintf('    Degree of annuitization              %8g ',lam*100))
disp(' ')
disp('  Steady State Values')
disp('  -------------------')
disp(sprintf('    Interest rate (%%)                    %8.3f ',100*(irate-1)))
disp(sprintf('    Wage rate                            %8.3f ',w))
disp(sprintf('    HH transfer                          %8.3f ',Tran))
disp(' ')
disp('  National Accounts (rel. to GDP)  ')
disp('  -------------------------------  ')
disp('  * Gross National Product/Income ')
disp(' ')  
disp(sprintf('    Consumption (C)                      %8.3f ',C/GDP))
disp(sprintf('    Tangible investment (XT)             %8.3f ',(XT1+XT2)/GDP))
disp(sprintf('      Corporate (XT1)                    %8.3f ',XT1/GDP))
disp(sprintf('      Noncorporate (XT2)                 %8.3f ',XT2/GDP))
disp(sprintf('    Government spending (G)              %8.3f ',G/GDP))
disp(     '                                            =====')
disp(sprintf('    GDP (Y-XI)                           %8.3f ',1.))
disp(' ')
disp(sprintf('    Labor income (wL)                    %8.3f ',w*(L1+L2)/GDP))
disp(sprintf('      Corporate (wL1)                    %8.3f ',w*L1/GDP))
disp(sprintf('      Noncorporate (wL2)                 %8.3f ',w*L2/GDP))
disp(sprintf('    Capital income (Y-wL-XI)             %8.3f ',(cprof+delt1*KT1+uprof+delt2*KT2)/GDP))
disp(sprintf('      Corporate (p1 Y1-wL1-XI1)          %8.3f ',(cprof+delt1*KT1)/GDP))
disp(sprintf('        Taxes                            %8.3f ',taupc*cprof/GDP))
disp(sprintf('        Net dividends                    %8.3f ',D/GDP))
disp(sprintf('        Undistributed                    %8.3f ',(XT1-delt1*KT1)/GDP))
disp(sprintf('        Depreciation                     %8.3f ',delt1*KT1/GDP))

disp(sprintf('      Noncorporate (p2 Y2-wL2-XI2)       %8.3f ',(uprof+delt2*KT2)/GDP))
disp(sprintf('        Taxes                            %8.3f ',taupu*uprof/GDP))
disp(sprintf('        Distributions before tax         %8.3f ',(1-taupu)*uprof/GDP))
disp(sprintf('        Depreciation                     %8.3f ',delt2*KT2/GDP))
disp(     '                                            =====')
disp(sprintf('    GDI = GDP                            %8.3f ',1.))
disp(sprintf('     Corporate value added               %8.3f ',(w*L1+cprof+delt1*KT1)/GDP))
disp(sprintf('     Noncorporate value added            %8.3f ',(w*L2+uprof+delt2*KT2)/GDP))
disp(' ')
disp(sprintf('    Addenda: intangible (XI)             %8.3f ',(XI1+XI2)/GDP))
disp(' ')
disp('  * Government Expenditures/Receipts ')
disp(' ')  
disp(sprintf('    Consumption of G&S                   %8.3f ',G/GDP)) 
disp(sprintf('    Transfers                            %8.3f ',(Tran+Tranj)/GDP)) 
disp(sprintf('      To retirees                        %8.3f ',Tranret/GDP))
disp(sprintf('      To workers                         %8.3f ',Tranwrk/GDP))
disp(sprintf('    Interest paid                        %8.3f ',max(Int/GDP,0))) 
disp(     '                                            =====')
disp(sprintf('    Expenditures                         %8.3f ',GovE/GDP))
disp(' ')
disp(sprintf('    Tax revenues                         %8.3f ',Rev/GDP))
disp(sprintf('      Consumption                        %8.3f ',tauc*C/GDP))
disp(sprintf('      Labor                              %8.3f ',taul*w*(L1+L2)/GDP))
disp(sprintf('      Dividends                          %8.3f ',taud*D/GDP))
disp(sprintf('      C-Profits                          %8.3f ',taupc*cprof/GDP))
disp(sprintf('      NC-Profits                         %8.3f ',taupu*uprof/GDP))
disp(sprintf('    Interest received                    %8.3f ',max(-Int/GDP,0))) 
disp(sprintf('    Accidental bequests                  %8.3f ',Beq/GDP))
disp(     '                                            =====')
disp(sprintf('    Receipts                             %8.3f ',GovR/GDP))
disp(     '                                            =====')
disp(sprintf('    Deficit=Expend.-Receipts             %8.3f ',Def/GDP))

disp(' ')
disp('  Fixed Assets (rel. to GDP)')
disp('  --------------------------')
disp(sprintf('    Tangible capital (KT'')               %8.3f ',(KT1eoy+KT2eoy)/GDP))
disp(sprintf('      Corporate (KT1'')                   %8.3f ',KT1eoy/GDP))
disp(sprintf('      Noncorporate (KT2'')                %8.3f ',KT2eoy/GDP))
disp(sprintf('    Addenda: intangible (KI'')            %8.3f ',(KI1eoy+KI2eoy)/GDP))
disp(sprintf('      Corporate (KI1'')                   %8.3f ',KI1eoy/GDP))
disp(sprintf('      Noncorporate (KI2'')                %8.3f ',KI2eoy/GDP))
disp(     '                                            =====')
disp(sprintf('    Total                                %8.3f ',(KT1eoy+KT2eoy+KI1eoy+KI2eoy)/GDP))
disp(' ')
disp('  Flow of Funds (rel. to GDP)')
disp('  ---------------------------')
disp('  * Household Balance Sheet')
disp(' ')
disp(sprintf('    Tangible capital (vT)                %8.3f ',(vT1+vT2)/GDP))
disp(sprintf('      Corporate (vT1)                    %8.3f ',vT1/GDP))
disp(sprintf('      Noncorporate (vT2)                 %8.3f ',vT2/GDP))
disp(sprintf('    Government Debt (B'')                 %8.3f ',Beoy/GDP))
disp(sprintf('    Addenda: intangible (vI)             %8.3f ',(vI1+vI2)/GDP))
disp(sprintf('      Corporate (vI1)                    %8.3f ',vI1/GDP))
disp(sprintf('      Noncorporate (vI2)                 %8.3f ',vI2/GDP))
disp(     '                                            =====')
disp(sprintf('    Net Worth                            %8.3f ',NWhh/GDP))
disp(' ')
disp('  * Government Balance Sheet')
disp(' ')
disp(sprintf('    Tangible capital (KT''-vT)            %8.3f ',(KT1eoy+KT2eoy-vT1-vT2)/GDP))
disp(sprintf('    Government Debt (-B'')                %8.3f ',-Beoy/GDP))
disp(sprintf('    Addenda: intangible (KI''-vI)         %8.3f ',(KI1eoy+KI2eoy-vI1-vI2)/GDP))
disp(     '                                            =====')
disp(sprintf('    Net Worth                            %8.3f ',NWgov/GDP))
disp(' ')
disp('  Welfare')
disp('  -------')
disp('  * Individual consumption')
disp(sprintf('      Average across ages                %8.3f ',mean(cj(1:J))))
disp(sprintf('      Sum mu_j*c_j                       %8.3f ',sum(mu.*cj)))
disp('  * Individual labor supply ')
disp(sprintf('      Average across working ages        %8.3f ',mean(lj(1:Jr-1))))
disp(sprintf('      Sum mu_j*l_j                       %8.3f ',sum(mu.*lj)))
disp('  * Values for computing welfare')
disp(sprintf('    Individual lifetime utility (U)      %8.3f ',U))
disp(sprintf('    Discount = sum beta^j*prod(ps(j)))   %8.3f ',Disc))
disp(' ')
disp('    % Welfare gain=(exp((U1-U2)/Discount)-1)*100 ')
disp(' ')
disp(' ')
disp(sprintf('  Detrended GDP:                         %8.3f ',GDP))


if iplot==1;

figure(1), plot(age,Aj)
           xlabel('Age')
           title('Cross-sectional asset holdings over the life cycle')
figure(2), plot(age,cj)
           xlabel('Age')
           title('Individual consumption over the life cycle')
figure(3), plot(age,lj)
           xlabel('Age')
           title('Individual labor supply over the life cycle')

end;
