%Intang      Generate time series for the 3-sector business cycle
%  _Simul    model with variations in technology parameters, fiscal
%            variables, an investment wedge, and nonbusiness activity.
%            Before running this code make sure to first run intang_steady 
%            and intang_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

intang_params

T0         = 1990;
bhat       = beta*(1+gamma)^(-sigma);
lshare1    = 1-thet1-phi1;
lshare2    = 1-thet2-phi2;
grate      = (1+eta)*(1+gamma);
param      = [beta;deltI;deltT;eta;gamma;phi1;phi2;psi;sigma;thet1;thet2];
%---------------------------------------------------------------------
% 2.  Compute equilibrium
%---------------------------------------------------------------------

load inputs.dat
[T,jnk]    = size(inputs);
disp('BC Accounting for Model with Intangibles and Non-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 final g&s sector?          ';
                '    TFP in intangible sector?         ';
                '    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?   ';
                '    Intangible 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;

z1t        = inputs(:,1);
z2t        = 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        = intang_resid(V0(:),[param;kTs;kIs;T;inputs(:);bderiv;uderiv]);
V          = secant('intang_resid',V0(:), ...
                    [param;kTs;kIs;T;inputs(:);bderiv;uderiv]);
V          = secant('intang_resid',real(V(:)), ...
                    [param;kTs;kIs;T;inputs(:);bderiv;uderiv]);
V1         = zeros(T,4);
V1(:)      = V;

kTtp       = V1(:,1);
kItp       = V1(:,2);
h1t        = V1(:,3);
h2t        = V1(:,4);
kTt        = [kTs;kTtp(1:T-1)];
kIt        = [kIs;kItp(1:T-1)];
kTtp2      = [kTtp(2:T);kTtp(T)];
kItp2      = [kItp(2:T);kItp(T)];
xTt        = grate*kTtp-(1-deltT)*kTt;
xmt        = xTt+xnt;
xIt        = grate*kItp-(1-deltI)*kIt;
hbt        = h1t+h2t;
ht         = hbt+hnt;
qxIybt     = lshare1/lshare2 *h2t./h1t;
k1Tt       = thet1*kTt./(thet1+thet2*qxIybt);
k2Tt       = kTt-k1Tt;
ybt        = k1Tt.^thet1.*kIt.^phi1.*(z1t.*h1t).^lshare1;
ymt        = ybt+ynt;
qxIt       = qxIybt.*ybt;
tybt       = ybt+qxIt;
qt         = qxIt./xIt;
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;
tprt       = (ymt+qxIt)./ht;
tbprt      = (ybt+qxIt)./hbt;
A1t        = z1t.^lshare1;
A2t        = z2t.^lshare2;
rTt        = thet1*ybt./k1Tt;
rIt        = (phi1*ybt+phi2*qxIt)./kIt;
bcompt     = lshare1*ybt+lshare2*qt.*xIt-(1-chit).*qt.*xIt;
proft      = rTt.*kTt+rIt.*kIt-deltT*kTt-taukt.*kTt-chit.*qt.*xIt;
gain       = qxIt;
val        = (1-taudt).*(1+gtauxt).*kTtp+(taugt+chit.*(1-taudt).*(1-taupt)+ ...
               (1-chit).*(1-gtauht)).*qt.*kItp;
dval       = [NaN;val(2:T)*grate-val(1:T-1)];

load ../../DATA/observables.dat
kbdat(1)   = kTs;
kdat(1)    = kTs+xns/(grate-1+deltT);
cdat       = observables(:,4)*gdps;
xmdat      = observables(:,5)*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);
  kdat1(i,1)   = ((1-deltT)*kdat(i)+xmdat(i))/grate;
  kdat(i+1,1)  = kdat1(i);
end;
kbdat1(T)= ((1-deltT)*kbdat(T)+xTdat(T))/grate;
kdat1(T) = ((1-deltT)*kdat(T)+xmdat(T))/grate;
mtfpbdat = ybdat./(kbdat.^(.33).*hbdat.^(.67));
mtfpdat  = gdpdat./(kdat.^(.33).*hdat.^(.67));

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 Hours and Model Prediction')
legend('US Hours','Model Hours','Location','Northwest')

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 GDP and Model Prediction')
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 Business Value Added and Model Prediction')
legend('US Business Value Added','Model Business Value Added', ...
       'Location','Northwest')

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 Measured Labor Productivity and Model Prediction')
legend('US GDP Per Hour','Model GDP Per Hour','Location','Northeast')

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 Measured Business Labor Productivity and Model Prediction')
legend('US Business Value Added Per Hour', ...
       'Model Business Value Added Per Hour','Location','Northwest')

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 Tangible Investment and Model Prediction')
legend('US Tangible Investment', ...
       'Model Tangible 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 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)]];
plot(fofavg(:,1),fofavg(:,2),'Linewidth',2,'Color','k')
line(T0:T1,dval./ymt.*(1+profratio(:,2))*100,'Linewidth',2, ...
     'Color','r','Marker','^')
line(T0:T1,dval./ymt*100,'Linewidth',2,'Color','b','Marker','*')
title( 'US Capital Gains Relative to GDP and Model Prediction')
legend('US Holding Gains/GDP', ...
       'Model Gains/GDP (with foreign)', ...
       'Model Gains/GDP (without foreign)','Location','Southwest')
axis([1990,2003,0,15])

NIPAbwt = [ ....
    1990  100.0000
    1991   95.6084
    1992   96.5825
    1993   96.2953
    1994   96.5439
    1995   96.7738
    1996   97.5143
    1997   99.2501
    1998  103.5733
    1999  106.2740
    2000  107.7714
    2001  104.2886
    2002  100.0662
    2003   97.9144]; % = (ssectcomp+.7*ssectpinc)./(dgdp.*grate.*pcps)/gdp0
                     %   from ./DATA/accounts.m or accounts2.m

figure(9)
set(gca,'Fontsize',14)
plot(NIPAbwt(:,1),NIPAbwt(:,2),'Linewidth',2,'Color','k')
line(T0:T1,bcompt/bcompt(1)*100,'Linewidth',2,'Color','r','Marker','^')
title( 'US Measured Compensation and Model Prediction')
legend('US NIPA Compensation','Model NIPA Compensation', ...
       'Location','Northwest')

figure(10)
set(gca,'Fontsize',14)
plot(T0:T1,hdat/hdat(1)*100,'Linewidth',2,'Color','k')
line(T0:T1,A1t/A1t(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,A2t/A2t(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, Final G&S','TFP, Intangibles', ...
       'Investment Wedge', ...
       'Location','Northwest')

figure(11)
set(gca,'Fontsize',14)
plot(T0:T1,ymt/ymt(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,(ymt+qxIt)/(ymt(1)+qxIt(1))*100,'Linewidth',2,'Color','b', ...
     'Marker','o')
ylabel('Index')
title( 'Model Outputs','Fontsize',14)
legend('Model GDP', ...
       'Model True Output','Location','Southwest')

figure(12)
set(gca,'Fontsize',14)
plot(T0:T1,ybt/ybt(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,tybt/tybt(1)*100,'Linewidth',2,'Color','b', ...
     'Marker','o')
ylabel('Index')
title( 'Model Business Value Added','Fontsize',14)
legend('Model Measured Business Value Added', ...
       'Model True Business Value Added','Location','Northwest')

figure(13)
set(gca,'Fontsize',14)
plot(T0:T1,mprt/mprt(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,tprt/tprt(1)*100,'Linewidth',2,'Color','b','Marker','o')
title( 'Model Labor Productivities','Fontsize',14)
legend('Model GDP Per Hour', ...
       'Model True Output Per Hour','Location','Southwest')

figure(14)
set(gca,'Fontsize',14)
plot(T0:T1,bprt/bprt(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,tbprt/tbprt(1)*100,'Linewidth',2,'Color','b','Marker','o')
ylabel('Index')
title( 'Model Business Labor Productivities','Fontsize',14)
legend('Model Measured Business Labor Productivity', ...
       'Model True Business Labor Productivity','Location','Northwest')

figure(15)
set(gca,'Fontsize',14)
plot(T0:T1,xmt/xmt(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,(xmt+qxIt)/(xmt(1)+qxIt(1))*100,'Linewidth',2,'Color','b', ...
       'Marker','o')
title( 'Model Investments','Fontsize',14)
legend('Model Measured Investment', ...
       'Model True Investment','Location','Southwest')

figure(16)
set(gca,'Fontsize',14)
plot(T0:T1,xTt/xTt(1)*100,'Linewidth',2,'Color','r','Marker','^')
line(T0:T1,(xTt+qxIt)/(xTt(1)+qxIt(1))*100,'Linewidth',2,'Color','b', ...
       'Marker','o')
title( 'Model Business Investments','Fontsize',14)
legend('Model Measured Business Investment', ...
       'Model True Business Investment','Location','Northwest')

figure(17)
set(gca,'Fontsize',14)
plot(T0:T1,qxIt./(ymt+qxIt)*100,'Linewidth',2,'Color','r','Marker','^')
title( 'Model Intangible Investment Share of Output (%)','Fontsize',14)

