%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;
xn      = model3(:,28);
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);
x       = xt+xn;
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;
prof    = (rt-delt-tauk).*kt+ri.*ki-chi*q.*xi;
div     = rt.*kt-xt+ri.*ki-chi*q.*xi-txk-txx-txp-txu;
gdp     = y+yn;
toty    = gdp+q.*xi;
t       = 1929:1928+length(gdp);
gdp0    = gdp(1);


%
% Load data from ../../accounts, where t=1929:1940 
%

load usdata.dat
tdat    = [1929:1940]';
gdpdat  = usdata(:,1);
cdat    = usdata(:,2);
xdat    = usdata(:,3);
gdat    = usdata(:,4);
ktdat   = usdata(:,5);
hdat    = usdata(:,6);
hndat   = usdata(:,7);
hbdat   = hdat-hndat;
ydat    = usdata(:,8);
yndat   = usdata(:,9);
xtdat   = usdata(:,10);
xndat   = usdata(:,11);
Wbdat   = usdata(:,12);
Wndat   = usdata(:,13);
Wdat    = Wbdat+Wndat;
Rbdat   = usdata(:,14);
Rndat   = usdata(:,15);
Rdat    = Rbdat+Rndat;
Depbdat = usdata(:,16);
Depndat = usdata(:,17);
Depdat  = Depbdat+Depndat;
vdat    = usdata(:,18);
taucdat = usdata(:,19);
tauhdat = usdata(:,20);
taukdat = usdata(:,21);
taupdat = usdata(:,22);
tauddat = usdata(:,23);
kndat   = usdata(:,25);
kdat    = ktdat+kndat;
kn(1)   = kndat(1);
deln    = .05;
for i=2:length(xn);
  kn(i,1) = ((1-deln)*kn(i-1)+xn(i-1))/((1+gz)*(1+gn));
end;
kdat2   = kdat;
ktdat2  = ktdat;
kndat2  = kndat;
delta   = delt;
for i=2:length(xndat);
  kdat2(i,1)  = ((1-delta)*kdat2(i-1)+xdat(i-1))/((1+gz)*(1+gn));
  ktdat2(i,1) = ((1-delta)*ktdat2(i-1)+xtdat(i-1))/((1+gz)*(1+gn));
  kndat2(i,1) = ((1-delta)*kndat2(i-1)+xndat(i-1))/((1+gz)*(1+gn));
end;
lwdgdat = cdat.*hdat./(gdpdat.*(1-hdat));
ewdgdat = gdpdat./(kdat.^.33.*hdat.^.67);
lwedge  = c.*h./(gdp.*(1-h));
ewedge  = gdp./((kt+kn).^.33.*h.^.67);

newtauh = 1-psi*(1+taucdat).*cdat.*hbdat./((1-alph1-thet1)*(1-hdat).*ydat) * ...
          h1(1)/(h1(1)+h2(1));
snewtauh= hptrend(newtauh,1);


%
% Plot results
%

figure(1)
plot(tdat,gdpdat,'Linewidth',2,'Color','b','Marker','o')
line(t,gdp/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
a=axis; axis([1929,1940,a(3),a(4)])
title('Actual and Predicted GDP','Fontsize',16);
legend('Actual','Predicted','Location','Southeast')

figure(2)
plot(tdat,cdat,'Linewidth',2,'Color','b','Marker','o')
line(t,c/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
a=axis; axis([1929,1940,a(3),a(4)])
title('Actual and Predicted Consumption','Fontsize',16);
legend('Actual','Predicted','Location','Northeast')

figure(3)
plot(tdat,xdat,'Linewidth',2,'Color','b','Marker','o')
line(t,x/gdp0,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
a=axis; axis([1929,1940,a(3),a(4)])
title('Actual and Predicted Investment','Fontsize',16);
legend('Actual','Predicted','Location','Southeast')

figure(4)
plot(tdat,hdat,'Linewidth',2,'Color','b','Marker','o')
line(t,h,'Linewidth',2,'Color','r','Marker','^')
set(gca,'Fontsize',16)
a=axis; axis([1929,1940,a(3),a(4)])
title('Actual and Predicted Hours','Fontsize',16);
legend('Actual','Predicted','Location','Northeast')

figure(5)
plot(t,xt,'Linewidth',2,'Color','r','Marker','^')
line(t,xi,'Linewidth',2,'Color','g','Marker','*')
set(gca,'Fontsize',16)
a=axis; axis([1929,1940,a(3),a(4)])
title('Predicted Business investments','Fontsize',16);
legend('Tangible','Intangible','Location','Southeast')
