%MOD3SIM      Plot simulations generated by model3.f90.  See details about
%             model economy in Fortran code.

%             Ellen McGrattan, 3-17-06
%             Revised, ERM, 8-20-09   

load model3.dat
load model3p.dat
s       = model3(:,1);
kt      = model3(:,2);
ki      = model3(:,3);
ktt     = model3(:,4);
kit     = model3(:,5);
c       = model3(:,6);
h1      = model3(:,7);
h2      = model3(:,8);
xt      = model3(:,9);
xi      = model3(:,10);
q       = model3(:,11);
y       = model3(:,12);
rt      = model3(:,13);
ri      = model3(:,14);
w       = model3(:,15);
g       = model3(:,16);
tauc    = model3(:,17);
tauh    = model3(:,18);
taub    = model3(:,19);
taux    = model3(:,20);
tauk    = model3(:,21);
taup    = model3(:,22);
tauu    = model3(:,23);
taud    = model3(:,24);
z1      = model3(:,25);
z2      = model3(:,26);
hn      = model3(:,27);
h       = h1+h2+hn;
hb      = h1+h2;
xn      = model3(:,28);
x       = xt+xn;
yn      = model3(:,29);
v       = model3(:,30);
d       = model3(:,31);
k1t     = model3(:,32);
k2t     = model3(:,33);
alph1   = model3p(1);
alph2   = model3p(2);
beta    = model3p(3);
chi     = model3p(4);
deli    = model3p(5);
delt    = model3p(6);
gn      = model3p(7);
gz      = model3p(8);
phi     = model3p(9);
psi     = model3p(10);
thet1   = model3p(11);
thet2   = model3p(12);
txc     = tauc.*c;
txh     = tauh.*w.*h;
se      = (1-chi)*taub.*q.*xi;
txx     = taux.*xt;
txk     = tauk.*kt;
txp     = taup.*((rt-delt-tauk).*kt+ri.*ki-chi*q.*xi);
txu     = tauu.*(xt-delt*kt);
txd     = taud.*(rt.*kt-xt-txk-txx-txp-txu);
tx      = txc+txh-se+txx+txk+txp+txu+txd;
tr      = c+xt+q.*xi+xn-rt.*kt-ri.*ki-w.*(h1+h2)+txc+txh-se+txx+ ...
          txk+txp+txu+txd-yn;
gdp     = y+yn;
lt      = length(gdp);
toty    = gdp+q.*xi;
t       = [2004:2003+lt]';
gdp0    = gdp(1);
grate   = (1+gz)*(1+gn);
dv      = [NaN;v(2:lt)*grate-v(1:lt-1)];


%
% Load data from observables.dat
%

load observables.dat
tdat    = observables(:,1);
gdpdat  = observables(:,2);
ydat    = observables(:,2)-observables(:,3);
cdat    = observables(:,4);
xdat    = observables(:,5);
xtdat   = observables(:,5)-observables(:,6);
hdat    = observables(:,7);
hbdat   = observables(:,7)-observables(:,8);
taucdat = observables(:,9);
tauhdat = observables(:,10);
taupdat = observables(:,11);
taukdat = observables(:,13);
hdat0   = hdat(1);
gdpdat0 = gdpdat(1);
hbdat0  = hbdat(1);
ydat0   = ydat(1);
T       = length(tdat);
kdat(1) = 4.15;
del     = xdat(1)/kdat(1)-gn-gz;
for i=2:T
  kdat(i) = ((1-del)*kdat(i-1)+xdat(i-1))/((1+gn)*(1+gz));
end;
kdat    = kdat(:);
k(1)    = 4.15;
for i=2:length(t);
  k(i)    = ((1-del)*k(i-1)+x(i-1))/((1+gn)*(1+gz));
end;
k       = k(:);
tfpdat  = gdpdat./(kdat.^(thet1+alph1).*hdat.^(1-thet1-alph1));
tfp     = gdp./(k.^(thet1+alph1).*h.^(1-thet1-alph1));
aggpdat = gdpdat./hdat*hdat0/gdpdat0*100;
buspdat = ydat./hbdat*hbdat0/ydat0*100;
aggp    = gdp./h*h(1)/gdp(1)*100;
busp    = y./hb*hb(1)/y(1)*100;


%
% Plot results
%

figure(1)
plot(tdat,gdpdat,'Linewidth',2,'Color','b','Marker','o')
line(t,gdp/gdp0,'Linewidth',2,'Color','r','Marker','^')
line(t,gdp/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,.88,1.02])
title('GDP','Fontsize',16);
legend('Actual','Predicted','Location','Southwest')

figure(2)
plot(tdat,hdat,'Linewidth',2,'Color','b','Marker','o')
line(t,h,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,.25,.29])
title('Per Capita Hours','Fontsize',16);
legend('Actual','Predicted','Location','Southwest')

figure(3)
plot(tdat,cdat,'Linewidth',2,'Color','b','Marker','o')
line(t,c/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,.72,.8])
title('Total Consumption','Fontsize',16);
legend('Actual','Predicted','Location','Southwest')

figure(4)
plot(tdat,xdat,'Linewidth',2,'Color','b','Marker','o')
line(t,x/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,.14,.24])
title('Total Investment','Fontsize',16);
legend('Actual','Predicted','Location','Southwest')

figure(5)
plot(tdat,aggpdat,'Linewidth',2,'Color','b','Marker','o')
line(t,aggp,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,95,100])
title('Aggregate Labor Productivity','Fontsize',16);
ylabel('Index')
legend('Actual','Predicted','Location','Southwest')

figure(6)
plot(tdat,buspdat,'Linewidth',2,'Color','b','Marker','o')
line(t,busp,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,95,102])
title('Business Labor Productivity','Fontsize',16);
ylabel('Index')
legend('Actual','Predicted','Location','Southwest')

figure(7)
plot(tdat,tfpdat./tfpdat(1)*100,'Linewidth',2,'Color','b','Marker','o')
line(t,tfp/tfp(1)*100,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,95,100])
title('Total Factor Productivity','Fontsize',16);
ylabel('Index')
legend('Actual','Predicted','Location','Southwest')

figure(8)
plot(t,xt/gdp0,'Linewidth',2,'Color','b','Marker','o')
line(t,q.*xi/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,.04,.12])
title('Predicted Business Investments','Fontsize',16);
legend('Tangible','Intangible','Location','Southwest')

figure(9)
plot(tdat,gdpdat/gdpdat(1)*100,'Linewidth',2,'Color','b','Marker','o')
line(t,toty/toty(1)*100,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,85,105])
title('US GDP and Predicted Total Output','Fontsize',16);
legend('GDP','Total Output','Location','Southwest')

%
% Implied intangible investment
%
qxidat   = psi*(1+taucdat).*cdat.*(hbdat)./((1-thet1-alph1)*(1-hdat).* ...
           (1-tauhdat))-ydat;
figure(10)
plot(tdat,qxidat,'Linewidth',2,'Color','b','Marker','o')
line(t,q.*xi,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
axis([2004,2011,.03,.12])
title('Intangible Investment','Fontsize',16);
legend('Residual','Predicted','Location','Southwest')


a1      = z1.^(1-thet1-alph1);
a2      = z2.^(1-thet2-alph2);
usdata  = [tdat,[gdpdat/gdpdat(1),cdat/cdat(1),xdat/xdat(1), ...
           hdat/hdat(1)]*100,aggpdat,buspdat];
inputs  = [t,tauc,tauh,a1,a2,[yn,xn,hn]/gdp0]; % Assume others held constant
outputs = [t,[gdp/gdp0,c/c(1),x/x(1),h/h(1)]*100,aggp,busp];
other   = [t,[xt/xt(1),q.*xi/q(1)/xi(1),(gdp+q.*xi)/(gdp0+q(1)*xi(1))]*100];


disp('               RESULTS FOR MODEL WITH INTANGIBLE CAPITAL              ')
disp('======================================================================')
disp(' ')
disp('  Parameters: ')
disp(' ')
disp('     Preferences ')
disp(sprintf('       Discount factor                       %17.4f',beta))
disp(sprintf('       Weight on log of Leisure              %17.4f',psi))
disp(' ')
disp('     Final G&S Sector Production')
disp(sprintf('       Tangible capital share                %17.4f',thet1))
disp(sprintf('       Intangible capital share              %17.4f',alph1))
disp(' ')
disp('     Intangible Sector Production')
disp(sprintf('       Tangible capital share                %17.4f',thet2))
disp(sprintf('       Intangible capital share              %17.4f',alph2))
disp(' ')
disp('     Depreciation rates')
disp(sprintf('       Tangible capital                      %17.4f',delt))
disp(sprintf('       Intangible capital                    %17.4f',deli))
disp(' ')
disp('     Growth rates')
disp(sprintf('       Technology                            %17.4f',gz))
disp(sprintf('       Population                            %17.4f',gn))
disp(' ')
disp('     Capital tax rates')
disp(sprintf('       Property                              %17.4f',tauk(1)))
disp(sprintf('       Profits                               %17.4f',taup(1)))
disp(sprintf('       Dividends                             %17.4f',taud(1)))
disp(' ')
disp(sprintf('     Share of Expensing by Firms             %17.4f',chi))
disp(' ')
disp('   Exogenous inputs:                                                  ')
disp(' ')
disp(['     Year      tauc      tauh       a1        a2        Yn', ...
     '        Xn        Hn  ']) 
disp(sprintf('%9g%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n',inputs(1:T,:)'))
disp('   Data:                                                              ')
disp('                                                      Productivity    ')
disp('     Year      GDP        C         X         H    Aggregate Business ')
disp(sprintf('%9g%10.4f%10.4f%10.4f%10.4f%10.2f%10.2f\n',usdata(1:T,:)'))
disp('   Comparable Model Outputs:                                          ')
disp('                                                      Productivity    ')
disp('     Year      GDP        C         X         H    Aggregate Business ')
disp(sprintf('%9g%10.4f%10.4f%10.4f%10.4f%10.2f%10.2f\n',outputs(1:T,:)'))
disp('   Model Investments and Total Output:                                ')
disp(' ')
disp('     Year      XT       Q*XI     GDP+Q*XI                             ')
disp(sprintf('%9g%10.4f%10.4f%10.4f\n',other(1:T,:)'))

