%BCA_Simul   Generate time series for a prototype business cycle model 
%            with variations in three wedges along with observed
%            fiscal variables. Before running this code make sure to 
%            first run bca_steady and bca_wedges. See details in the 
%            Appendix to:
%
%               Unmeasured Investment and the 1990s US Hours Boom
%               by Ellen McGrattan and Ed Prescott
                                                                                
%            Ellen McGrattan, 11-1-05
%            Revised, ERM, 8-6-07
%---------------------------------------------------------------------
% 1. Default parameters 
%
clear
bca_params

T0         = 1990;
bhat       = beta*(1+gamma)^(-sigma);
thet1      = 1-theta;
delt1      = 1-delta;
grate      = (1+eta)*(1+gamma);
param      = [beta;delta;eta;gamma;psi;sigma;theta];

%---------------------------------------------------------------------
% 2.  Compute equilibrium
%---------------------------------------------------------------------

load inputs.dat
T          = length(inputs);
clc
disp('Business Cycle Accounting for 1-Sector Growth Model')
disp('---------------------------------------------------')
iexog = input('  Set any exogenous variables to 1990 level? (1=yes,0=no) ');
if iexog==1;
  disp('    This variable set to 1990 level? (1=yes,0=no)')
  labels     = ['    TFP?                              ';
                '    Labor wedge?                      ';
                '    Investment wedge?                 ';
                '    Government spending?              ';
                '    Government tax --Consumption?     ';
                '    Government tax --Labor?           ';
                '    Government tax --Investment?      ';
                '    Government tax --Property?        ';
                '    Government tax --Profits?         ';
                '    Government tax --Distributions?   '];
  for i=1:10;
    if any(abs(inputs(:,i)-inputs(1,i))>1e-3);
      ifix     = input(labels(i,:));
      if ifix==1; 
         inputs(:,i) = inputs(1,i)*ones(T,1); 
      end;
    end;
  end;
end;
At         = inputs(:,1);
lwedget    = inputs(:,2);
xwedget    = inputs(:,3);
gt         = inputs(:,4);
tauct      = inputs(:,5);
gtauht     = inputs(:,6);
gtauxt     = inputs(:,7);
taukt      = inputs(:,8);
taupt      = inputs(:,9);
taudt      = inputs(:,10);
tauht      = 1-(1-gtauht).*lwedget;
tauxt      = (1+gtauxt)./xwedget-1;

load vinitial.dat
V0         = vinitial;
deriv      = V0(T,1)/V0(T-1,1);
res        = bca_resid(V0(:),[param;ks;T;inputs(:);deriv]);
V          = secant('bca_resid',V0(:),[param;ks;T;inputs(:);deriv]);
V1         = zeros(T,2);
V1(:)      = V;

ktp        = V1(:,1);
ht         = V1(:,2);
kt         = [ks;ktp(1:T-1)];
xt         = grate*ktp-(1-delta)*kt;
yt         = At.*kt.^theta.*ht.^thet1;
ct         = yt-xt-gt;
lt         = 1-ht;
uct        = ct.^(-sigma).*lt.^(psi*(1-sigma));
ult        = psi*ct.^(1-sigma).*lt.^(psi*(1-sigma)-1);
prt        = yt./ht;
rt         = theta*yt./kt;
proft      = rt.*kt-delta*kt-tauks*kt;
val        = (1-tauds)*(1+tauxs)*ktp;
dval       = [NaN;val(2:T)*grate-val(1:T-1)];
                                                                                
load ../../DATA/observables.dat
ydat       = observables(:,2)*gdps;
cdat       = observables(:,4)*gdps;
hdat       = observables(:,7);
xdat       = observables(:,5)*gdps;
tauhdat    = observables(:,11);
kdat(1)    = ks;
for i=1:T-1;
  kdat1(i,1)  = (delt1*kdat(i)+xdat(i))/grate;
  kdat(i+1,1) = kdat1(i);
end;
kdat1(T)   = (delt1*kdat(T)+xdat(T))/grate;
mtfpt      = gdps*ydat./(kdat.^theta.*hdat.^thet1);
prdat      = ydat./hdat;

T1         = T0+T-1;

figure(1)
set(gca,'Fontsize',14)
plot(T0:T1,hdat/hdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,ht/ht(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Per Capita Hours and Model Prediction')
legend('US Hours','Model Hours','Location','Southwest')

figure(2)
set(gca,'Fontsize',14)
plot(T0:T1,ydat/ydat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,yt/yt(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Per Capita Real GDP and Model Prediction (/1.02^t)')
legend('US GDP','Model GDP','Location','Southwest')

figure(3)
set(gca,'Fontsize',14)
plot(T0:T1,prdat/prdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,prt/prt(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Labor Productivity and Model Prediction (/1.02^t)')
legend('US GDP Per Hour','Model GDP Per Hour','Location','Southwest')

figure(4)
set(gca,'Fontsize',14)
plot(T0:T1,xdat/xdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,xt/xt(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Per Capita Real Investment and Model Prediction (/1.02^t)')
legend('US Investment','Model Investment', ...
       'Location','Southwest')

figure(5)
set(gca,'Fontsize',14)
plot(T0:T1,cdat/cdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,ct/ct(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Per Capita Real Consumption and Model Prediction (/1.02^t)')
legend('US Consumption','Model Consumption','Location','Southeast')

figure(6)
set(gca,'Fontsize',14)
plot(T0:T1,hdat/hdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,At/At(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,lwedget/lwedget(1)*100,'Linewidth',2,'Color','b','Marker','*')
line(T0:T1,xwedget/xwedget(1)*100,'Linewidth',2,'Color','c','Marker','+')
title( 'US Hours and the BCA Wedges','Fontsize',14)
legend('US Hours','TFP','Labor Wedge','Investment Wedge', ...
       'Location','Northwest')

