%RIntang     Generate time series for the 2-sector business cycle
%  _Simul    model with intangible and variations in three wedges,
%            fiscal variables, and nonbusiness activity. Before 
%            running this code make sure to first run rintang_steady
%            and rintang_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-13-07
%---------------------------------------------------------------------
% 1. Default parameters 

rintang_params;

T0         = 1990;
bhat       = beta*(1+gamma)^(-sigma);
lshare     = 1-theta-phi;
grate      = (1+eta)*(1+gamma);
param      = [beta;deltI;deltT;eta;gamma;phi;psi;sigma;theta];
%---------------------------------------------------------------------
% 2.  Compute equilibrium
%---------------------------------------------------------------------

load inputs.dat
[T,jnk]    = size(inputs);
clc
disp('BC Accounting for Model with Intangibles and Neutral Technologies')
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?   ';
                '    Government tax --Capital gains?   ';
                '    %Intangibles financed by profits? '];
  for i=1:12;
    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);
taugt      = inputs(:,11);
chit       = inputs(:,12);
ynt        = inputs(:,13);
xnt        = inputs(:,14);
hnt        = inputs(:,15);

load vinitial.dat
V0         = vinitial;
bderiv     = V0(T,1)/V0(T-1,1);
uderiv     = V0(T,2)/V0(T-1,2);
res        = rintang_resid(V0(:),[param;kTs;kIs;T;inputs(:);bderiv;uderiv]);
V          = secant('rintang_resid',V0(:),[param;kTs;kIs;T;inputs(:);
                     bderiv;uderiv]);
V1         = zeros(T,3);
V1(:)      = V;

kTtp       = V1(:,1);
kItp       = V1(:,2);
hbt        = V1(:,3);
kTt        = [kTs;kTtp(1:T-1)];
kIt        = [kIs;kItp(1:T-1)];
xTt        = grate*kTtp-(1-deltT)*kTt;
xmt        = xTt+xnt;
xIt        = grate*kItp-(1-deltI)*kIt;
ht         = hbt+hnt;
tybt       = At.*kTt.^theta.*kIt.^phi.*hbt.^lshare;
ybt        = tybt-xIt;
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;
rTt        = theta*tybt./kTt;
rIt        = phi*tybt./kIt;
bcompt     = lshare*tybt;
proft      = rTt.*kTt+rIt.*kIt-deltT*kTt-taukt.*kTt-xIt;
gain       = xIt;
val        = (1-taudt).*(1+gtauxt).*kTtp+(taugt+(1-taudt).*(1-taupt)).*kItp;
dval       = [NaN;val(2:T)*grate-val(1:T-1)];

load ../../DATA/observables.dat
kbdat(1)   = kTs;
cdat       = observables(:,4)*gdps;
xTdat      = (observables(:,5)-observables(:,6))*gdps;
ybdat      = (observables(:,2)-observables(:,3))*gdps;
gdpdat     = observables(:,2)*gdps;
hdat       = observables(:,7);
hbdat      = hdat-hnt;
mprdat     = gdpdat./hdat;
bprdat     = ybdat./hbdat;
for i=1:T-1;
  kbdat1(i,1)  = ((1-deltT)*kbdat(i)+xTdat(i))/grate;
  kbdat(i+1,1) = kbdat1(i);
end;
kbdat1(T)  = ((1-deltT)*kbdat(T)+xTdat(T))/grate;
mtfpdat    = gdpdat./(kbdat.^(.33).*hbdat.^(.67));

T1         = T0+T-1;
figure(1)
set(gca,'Fontsize',14)
plot(T0:T1,hdat,'Linewidth',2,'Color','k')
line(T0:T1,ht,'Linewidth',2,'Color','r','Marker','^')
title( 'US Hours and Model Prediction')
legend('US Hours','Model Hours','Location','Southwest')

figure(2)
set(gca,'Fontsize',14)
plot(T0:T1,gdpdat,'Linewidth',2,'Color','k')
line(T0:T1,ymt,'Linewidth',2,'Color','r','Marker','^')
title( 'US GDP and Model Prediction')
legend('US GDP','Model GDP','Location','Southwest')

figure(3)
set(gca,'Fontsize',14)
plot(T0:T1,ybdat,'Linewidth',2,'Color','k')
line(T0:T1,ybt,'Linewidth',2,'Color','r','Marker','^')
title( 'US Business Value-added and Model Prediction')
legend('US Business Value-added','Model Business Value-added', ...
       'Location','Southwest')

figure(4)
set(gca,'Fontsize',14)
plot(T0:T1,mprdat,'Linewidth',2,'Color','k')
line(T0:T1,mprt,'Linewidth',2,'Color','r','Marker','^')
title( 'US Measured Labor Productivity and Model Prediction')
legend('US GDP Per Hour','Model GDP Per Hour','Location','Southwest')

figure(5)
set(gca,'Fontsize',14)
plot(T0:T1,bprdat,'Linewidth',2,'Color','k')
line(T0:T1,bprt,'Linewidth',2,'Color','r','Marker','^')
title( 'US Measured Business Labor Productivity and Model Prediction')
legend('US Business Value Added Per Hour', ...
       'Model Business Value Added Per Hour','Location','Southwest')

figure(6)
set(gca,'Fontsize',14)
plot(T0:T1,xTdat,'Linewidth',2,'Color','k')
line(T0:T1,xTt,'Linewidth',2,'Color','r','Marker','^')
title( 'US Business Tangible Investment and Model Prediction')
legend('US Business Tangible Investment', ...
       'Model Business Tangible Investment', ...
       'Location','Southwest')

figure(7)
set(gca,'Fontsize',14)
plot(T0:T1,cdat,'Linewidth',2,'Color','k')
line(T0:T1,ct,'Linewidth',2,'Color','r','Marker','^')
title( 'US Consumption and Model Prediction')
legend('US Consumption','Model Consumption','Location','Southwest')

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','*')
legend('FOF Holding Gains/GDP','Model (with foreign gains)', ...
       'Model (without foreign gains)','Location','Southwest')
axis([1990,2004,-.05,.3])

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')

figure(10)
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 Hours and Model Prediction (1990=100)')
legend('US Hours','Model Hours','Location','Southeast')
