%NoIntang    Generate time series for a prototype 2-sector business cycle 
%  _Simul    model with variations in three wedges along with observed
%            fiscal variables and nonbusiness activity. Before running 
%            this code make sure to first run nointang_steady and 
%            nointang_wedges. See details in the Appendix to:
%
%               Unmeasured Investment and the 1990s US Hours Boom
%               by Ellen McGrattan and Ed Prescott
                                                                                
%            Ellen McGrattan, 7-26-05
%            Revised, ERM, 8-6-07
%---------------------------------------------------------------------
% 1. Default parameters 
%

clear
nointang_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 2-Sector Growth Model, No Intangible')
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 in business sector            ';
                '    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);
ynt        = inputs(:,11);
xnt        = inputs(:,12);
hnt        = inputs(:,13);
tauht      = 1-(1-gtauht).*lwedget;
tauxt      = (1+gtauxt)./xwedget-1;

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

kbtp       = V1(:,1);
hbt        = V1(:,2);
kbt        = [kbs;kbtp(1:T-1)];
xbt        = grate*kbtp-delt1*kbt;
xmt        = xbt+xnt;
ht         = hbt+hnt;
ybt        = At.*kbt.^theta.*hbt.^thet1;
ymt        = ybt+ynt;
ct         = ymt-xmt-gt;
lt         = 1-ht;
uct        = ct.^(-sigma).*lt.^(psi*(1-sigma));
ult        = psi*ct.^(1-sigma).*lt.^(psi*(1-sigma)-1);
mprt       = ymt./ht;
bprt       = ybt./hbt;
rbt        = theta*ybt./kbt;
bcompt     = thet1*ybt;
proft      = rbt.*kbt-delta*kbt-tauks*kbt;
val        = (1-tauds)*(1+tauxs)*kbtp;
dval       = [NaN;val(2:T)*grate-val(1:T-1)];

load ../../DATA/observables.dat
gdpdat     = observables(:,2)*gdps;
cdat       = observables(:,4)*gdps;
xmdat      = observables(:,5)*gdps;
xbdat      = (observables(:,5)-observables(:,6))*gdps;
ybdat      = (observables(:,2)-observables(:,3))*gdps;
hdat       = observables(:,7);
hbdat      = hdat-hnt;
mprdat     = gdpdat./hdat;
bprdat     = ybdat./hbdat;
kdat(1)    = kbs;
for i=1:T-1;
  kdat1(i,1)  = (delt1*kdat(i)+xbdat(i))/grate;
  kdat(i+1,1) = kdat1(i);
end;
kdat1(T)   = (delt1*kdat(T)+xbdat(T))/grate;
mtfpdat    = ybdat./(kdat.^theta.*hbdat.^thet1);

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,gdpdat/gdpdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,ymt/ymt(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,ybdat/ybdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,ybt/ybt(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Per Capita Real Business Value-Added and Model Prediction (/1.02^t)')
legend('US Business Value Added','Model Business Value Added', ...
       'Location','Southwest')

figure(4)
set(gca,'Fontsize',14)
plot(T0:T1,mprdat/mprdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,mprt/mprt(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(5)
set(gca,'Fontsize',14)
plot(T0:T1,bprdat/bprdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,bprt/bprt(1)*100,'Linewidth',2,'Color','r','Marker','^')
ylabel('Index')
title( 'US Business Labor Productivity and Model Prediction (/1.02^t)')
legend('US Business Labor Productivity', ...
       'Model Business Labor Productivity','Location','Southwest')

figure(6)
set(gca,'Fontsize',14)
plot(T0:T1,xmdat/xmdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,xmt/xmt(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(7)
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(8)
set(gca,'Fontsize',14)
profratio=[ ...
    1990    0.3518
    1991    0.3240
    1992    0.2839
    1993    0.2616
    1994    0.2340
    1995    0.2409
    1996    0.2251
    1997    0.2171
    1998    0.2285
    1999    0.2579
    2000    0.3580
    2001    0.4313
    2002    0.2896
    2003    0.2822]; % =(natprof-prof)./(prof-ctax);
                     % = adjustment so that foreign/domestic gains
                     %   is equal to foreign/domestic after-tax profits
fofavg = [[1953:2003]',[5.98*ones(42,1);12.14*ones(9,1)]/100];
plot(fofavg(:,1),fofavg(:,2),'Linewidth',2,'Color','k')
line(T0:T1,dval./ymt.*(1+profratio(:,2)),'Linewidth',2,'Color','r', ...
     'Marker','^')
line(T0:T1,dval./ymt,'Linewidth',2,'Color','b','Marker','*')
title( 'US Average Real Holding Gains and Model Prediction')
legend('FOF Holding Gains/GDP', ...
       'Model Gains/GDP (with foreign)', ...
       'Model Gains/GDP (without foreign)','Location','Southwest')
axis([1990,2003,-.02,.15001])

figure(9)
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')

