%COMPARE   plot data and model after running benchsim.
%

%          Ellen McGrattan
%          Revised, erm, 6-16-08

clear,clc

%
% Run the benchmark model
%
disp('Which model to simulate?')
disp('   1 = benchmark')
imodel = input('   2 = capacity utilization');
iskip  = 1;  % skip additional model simulations
if imodel==1;
  cd bench
  benchsim39
else
  cd caputil
  caputsim39
end;
cd ..

%
% Load data compiled by Braun and McGrattan for their
% paper in the NBER Macro annual 1993.
%
t         = [t0:t1]';
T         = t-1788;
t1941     = find(t==1941);
t1946     = find(t==1946);
t1947     = find(t==1947);
t1948     = find(t==1948);
tt0       = t0;    % First year plotted
tt1       = 1946;  % Last year plotted

load data/nipr_us.dat                   % billions of 1982 dollars
load data/nipn_us.dat                   % billions of dollars
load data/ibtx_us.dat                   % billions of dollars
load data/cons_us.dat                   % billions of 1982 dollars
load data/capr_us.dat                   % billions of 1982 dollars
load data/ginv_us.dat                   % millions of 1982 dollars
load data/fpur_us.dat                   % billions of dollars or 1982$
load data/spur_us.dat                   % billions of dollars or 1982$
load data/popu_us.dat                   % thousands
load data/hour_us.dat                   % millions (for series used)
load data/empl_us.dat                   % thousands (for series used)
load data/full_us.dat                   % thousands
load data/mhrs_us.dat                   % millions
load data/nfhr_us.dat                   % millions
load data/comp_us.dat                   % millions of dollars
load data/wind_us.dat                   % millions of dollars
load data/prop_us.dat                   % millions of dollars
load data/ivan_us.dat                   % millions of dollars
load data/ninc_us.dat                   % billions of dollars
load data/kinv_us.dat                   % billions of dollars or 1982$
load data/mtax_us.dat                   % percent

%
% Pull relevant columns from raw data files.
%
pop16     = popu_us(T,5)/1000000;       % Population over 16
gnpdef82  = nipr_us(T,8)/100;           % GNP deflator 1982$
gnp82     = nipr_us(T,7);               % Real GNP, 1982$
c82       = cons_us(T,2);               % Real personal consumption exp., 1982$
cd82      = cons_us(T,3);               % Real durable expenditures, 1982$
cnd82     = cons_us(T,4);               % Real nondurable expenditures, 1982$
cs82      = cons_us(T,5);               % Real service expenditures, 1982$
g82       = nipr_us(T,6);               % Real government purchases, 1982$
if82      = nipr_us(T,3);               % Real fixed investment, 1982$
ii82      = nipr_us(T,4);               % Real change in inventories, 1982$
nx82      = nipr_us(T,5);               % Real net exports, 1982$
fedgi82   = ginv_us(T,2)/1000;          % Real federal govt. investment, 1982$
stlgi82   = ginv_us(T,9)/1000;          % Real state&local govt. inv., 1982$
fedme82   = ginv_us(T,4)/1000;          % Real federal mil. equip. inv., 1982$
fedms82   = ginv_us(T,8)/1000;          % Real federal mil. bldgs. inv., 1982$
fedcomp82 = fpur_us(T,10);              % Real compensation, federal employees
stlcomp82 = spur_us(T,3);               % Real compensation, s&l employees
thrs      = mhrs_us(T,2)/1000;          % Manhours, national economy
mhrs      = mhrs_us(T,5)/1000;          % Manhours, military
fhrs      = mhrs_us(T,8)/1000;          % Manhours, farm
mfhrs     = nfhr_us(T,4)/1000;          % Manhours, manufacturing
temp      = empl_us(T,2)/1000000;       % Persons engaged, national economy
memp      = empl_us(T,4)/1000000;       % Persons engaged, military
tfte      = full_us(T,2)/1000000;       % Full-time equiv. emp., national ec.
mfte      = full_us(T,17)/1000000;      % Full-time equiv. emp., military
kpnre82   = capr_us(T-1,2);             % Real stock, priv. nonres. eq., 1982$
kpnrs82   = capr_us(T-1,3);             % Real stock, priv. nonres. st., 1982$
kpr82     = capr_us(T-1,4);             % Real stock, priv. residential, 1982$
kgnre82   = capr_us(T-1,5);             % Real stock, govt. nonres. eq., 1982$
kgnrs82   = capr_us(T-1,6);             % Real stock, govt. nonres. st., 1982$
kgr82     = capr_us(T-1,7);             % Real stock, govt. residential, 1982$
kcd82     = capr_us(T-1,8);             % Real stock, consumer durables, 1982$
kmile82   = capr_us(T-1,9);             % Real stock, military equip., 1982$
kmils82   = capr_us(T-1,10);            % Real stock, military struc., 1982$
ngnp      = nipn_us(T,7);               % Nominal GNP 
nsaletax  = ibtx_us(T,3)+ibtx_us(T,4);  % Nominal excise plus sales tax
ncomp     = comp_us(T,2)/1000;          % Nominal compensation, total
nagcomp   = comp_us(T,5)/1000;          % Nominal compensation, agriculture
nmfcomp   = comp_us(T,8)/1000;          % Nominal compensation, manufacturing
nfedcomp  = fpur_us(T,5);               % Nominal compensation, fed. employees
nstlcomp  = spur_us(T,2);               % Nominal compensation, s&l employees
nmilws    = fpur_us(T,12);              % Nominal military wages 
nmilcomp  = [wind_us(1:159,17);
             comp_us(160:210,17)]/1000; % ---Nominal military compensation---
nmilcomp  = nmilcomp(T);                 
wgt       = nmilcomp(t1948)./nmilws(t1948);
nmilcomp(t1946) = nmilcomp(t1946).* wgt;
nmilcomp(t1947) = nmilcomp(t1947).* wgt;% ---Nominal military compensation---
nnfprop   = (prop_us(T,2)- ...
             prop_us(T,3)+ ...          
             ivan_us(T,11))/1000;       % Nominal nonfarm propreitors' income
nagprop   = (prop_us(T,12)+ ...
             prop_us(T,3))/1000;        % Nominal farm propreitors' income
nmfprop   = (prop_us(T,6)+ ...           
             ivan_us(T,14))/1000;       % Nominal manufacturing prop. income
nprop     = ninc_us(T,4);               % Nominal proprietors' income
taukdat   = mtax_us(T,6)/100;           % Tax rate on capital
tauldat   = mtax_us(T,10)/100;          % Tax rate on labor
%tauldat  = mtax_us(T,5)/100;           % Tax rate on labor (Barro & Sahasakul)
kinv82    = kinv_us(T-1,2);
j         = find(~isnan(kinv82));
for i=j(1):-1:2; 
  kinv82(i-1) = kinv82(i)-ii82(i-1);
end;

% Use data to construct variables to be compared with data.
%

milws82   = nmilws.* ...            
            fedcomp82./nfedcomp;        % Real military wages, 1982$
saletax82 = nsaletax./gnpdef82;         % Real sales taxes, 1982$
nmvadef82 = (ngnp-nmilws)./(gnp82-fedcomp82.* ...
            nmilws./nfedcomp);          % Non-military deflator
pvadef82  = (ngnp-nfedcomp-nstlcomp)./(gnp82- ...
            fedcomp82-stlcomp82);       % Private value-added deflator
chrs      = thrs-mhrs;                  % Manhours, civilian
nfhrs     = chrs-fhrs;                  % Manhours, civilian non-farm
cemp      = temp-memp;                  % Persons engaged, civilian
cfte      = tfte-mfte;                  % Full-time equiv. emp., civilian
fehrs     = chrs.* ( cfte./cemp );      % Manhours, civilian FTEs
rgnppc    = (gnp82-saletax82+ ...
             .04*kcd82)./pop16;         % Real per capita GNP
rcppc     = (cnd82+cs82-saletax82.* ...
            (cnd82+cs82)./(cd82+cnd82+ ...
             cs82)+.04*kcd82)./pop16;   % Real per capita private consumption
rippc     = (cd82-saletax82.* ...
             cd82./(cd82+cnd82+cs82)+ ...
             if82+ii82)./pop16;         % Real per capita priv. dom. investment
rgpc      = (g82+nx82)./pop16;          % Real per capita govt spending+net ex.
rigpc     = (fedgi82+stlgi82-fedme82- ...
             fedms82)./pop16;           % Real per capita nonmil. govt. inv.
rmwspc    = milws82./pop16;             % Real per capita military wages
rcgpc     = rgpc-rigpc-rmwspc;          % Real per capita govt. consumption
rnxpc     = nx82./pop16;                % Real per capita net exports
rcdpc     = cd82./pop16;                % Real per capita cons. durables
twhrpc    = thrs./pop16 /52;            % Weekly hours per capita, total
mwhrpc    = mhrs./pop16 /52;            % Weekly hours per capita, military
cwhrpc    = chrs./pop16 /52;            % Weekly hours per capita, civilian
nfwhrpc   = nfhrs./pop16 /52;           % Weekly hours per capita, nonfarm 
mfwhrpc   = mfhrs./pop16 /52;           % Weekly hours per capita, manufacturing
temppc    = temp./pop16;                % Total Employees/population 16+
memppc    = memp./pop16;                % Military employees/ population 16+
cemppc    = cemp./pop16;                % Civilian employees/ population 16+
twhrwk    = thrs./temp /52;             % Weekly hours per worker, total
mwhrwk    = mhrs./memp /52;             % Weekly hours per worker, military
cwhrwk    = chrs./cemp /52;             % Weekly hours per worker, civilian
rkppc_bea = (kcd82+kpnre82+kpnrs82+ ...
             kpr82+kinv82)./pop16;      % Real stock, private capital (BEA)
rkgpc_bea = (kgnre82+kgnrs82+kgr82- ...
             kmile82-kmils82)./ ...
             pop16;                     % Real stock, nonmil. govt. cap. (BEA)
ncivcomp  = ncomp-nmilcomp;             % Nominal compensation, civilian
nnfcomp   = ncomp-nmilcomp-nagcomp;     % Nominal compensation, civ. nonfarm

%
% Detrend and normalize variables to make comparable to model
%

GNP46     = rgnppc(t1946);
Y46       = rcppc(t1946)+rippc(t1946)+rcgpc(t1946)+rigpc(t1946);
cpdat     = rcppc ./((1+gz).^([t0:t1]'-1946))/Y46;
cgdat     = rcgpc ./((1+gz).^([t0:t1]'-1946))/Y46;
ipdat     = rippc ./((1+gz).^([t0:t1]'-1946))/Y46;
igdat     = rigpc ./((1+gz).^([t0:t1]'-1946))/Y46;
gdat      = rgpc  ./((1+gz).^([t0:t1]'-1946))/Y46;
ydat      = (rcppc+rippc+rcgpc+rigpc)./((1+gz).^([t0:t1]'-1946))/Y46;
mwsdat    = rmwspc./((1+gz).^([t0:t1]'-1946))/Y46;
nxdat     = rnxpc ./((1+gz).^([t0:t1]'-1946))/Y46;
cddat     = rcdpc ./((1+gz).^([t0:t1]'-1946))/Y46;
gnpdat    = rgnppc./((1+gz).^([t0:t1]'-1946))/Y46;
adat      = memppc;
ndat      = cemppc./(1-adat);
hdat      = cwhrwk/84;
lcdat     = ndat.*hdat;
lpdat     = (1-adat).*ndat.*hdat;
ldat      = lpdat+adat*50/84;
kpdat_bea = rkppc_bea ./((1+gz).^([t0:t1]'-1946))/Y46;
kgdat_bea = rkgpc_bea ./((1+gz).^([t0:t1]'-1946))/Y46;
kpdat     = kpdat_bea;
kgdat     = kgdat_bea;
for i=2:length(ipdat);
  kpdat(i) = ((1-delta)*kpdat(i-1)+ipdat(i-1))/((1+gp)*(1+gz));
  kgdat(i) = ((1-delta)*kgdat(i-1)+igdat(i-1))/((1+gp)*(1+gz));
end;

z1dat_bea = (ydat./((kpdat_bea+kgdat_bea).^(theta).*lpdat.^(1-theta))).^ ...
            (1/(1-theta));
z1dat     = (ydat./((kpdat+kgdat).^(theta).*lpdat.^(1-theta))).^(1/(1-theta));
z2dat_bea = (ydat./((kpdat_bea+kgdat_bea).^(theta).*(ndat.*(1-adat)).^ ...
            (1-theta).*hdat)).^(1/(1-theta));
z2dat     = (ydat./((kpdat+kgdat).^(theta).*(ndat.*(1-adat)).^ ...
            (1-theta).*hdat)).^(1/(1-theta));

pshare    = 1-theta;
wcivdat   = ((ncivcomp+pshare*nprop)./ ...
             (chrs.*nmvadef82))./((1+gz).^([t0:t1]'-1946));
wnfdat    = ((nnfcomp+pshare*nnfprop)./ ...
             (nfhrs.*nmvadef82))./((1+gz).^([t0:t1]'-1946));
wmfdat    = ((nmfcomp+pshare*nmfprop)./ ...
             (mfhrs.*nmvadef82))./((1+gz).^([t0:t1]'-1946));
wfedat    = (ncivcomp./ ... 
             (fehrs.*nmvadef82))./((1+gz).^([t0:t1]'-1946));
wagdat    = ((nagcomp+pshare*nagprop)./ ...
             (fhrs.*nmvadef82))./((1+gz).^([t0:t1]'-1946));

mpkdat    = theta.*ydat./(kpdat_bea+kgdat_bea);
mpldat    = (1-theta)*ydat./lpdat;  % =w in benchmark model
                                    % =w(h)/h in capacity utilization model
atwcivdat = (1-tauldat).*wcivdat;
atwnfdat  = (1-tauldat).*wnfdat;
atwfedat  = (1-tauldat).*wfedat;
atwmfdat  = (1-tauldat).*wmfdat;
atr2ldat  = (1-tauldat).*mpldat;
atr2kdat  = (1-taukdat).*(mpkdat-delta)*100;
t46       = find(t==1946);
t57       = find(t==1957);
t60       = find(t==1960);
lpscale   = mean(lpdat(t46:t60));
atr2lsc   = mean(atr2ldat(t46:t60));
atwcivsc  = mean(atwcivdat(t46:t60));
atwnfsc   = mean(atwnfdat(t46:t60));
atwfesc   = mean(atwfedat(t46:t60));
atwmfsc   = mean(atwmfdat(t46:t57));
r2lsc     = mean(mpldat(t46:t60));
wcivsc    = mean(wcivdat(t46:t60));
wnfsc     = mean(wnfdat(t46:t60));
wfesc     = mean(wfedat(t46:t60));
wmfsc     = mean(wmfdat(t46:t57));
lprddat   = ydat./lpdat *lpscale;
lprdmod   = ymod./lpmod *lpscale;
mplmod    = (1-theta)*ymod./lpmod;  
atr2lmod  = (1-tauldat).*mplmod;
htable    = [t,lpmod/lpscale,lpdat/lpscale];
wtable    = [t,atr2lmod/atr2lsc,atr2ldat/atr2lsc,atwcivdat/atwcivsc, ...
             atwnfdat/atwnfsc,atwfedat/atwfesc,atwmfdat/atwmfsc];
%
% Fixed point experiments where we adjust taul until either
%     (1)  pre-tax wages match: w_model(taul) = w_data
% or  (2)  hours match: h_model(taul) = h_data
%
iexp = 0; 
if iexp;
disp(' ')
texp  = [t,htable(:,2:3),mplmod/r2lsc,mpldat/r2lsc,wnfdat/wnfsc, ...
         wmfdat/wmfsc,wcivdat/wcivsc];
texp2 = [t,lpmod/mean(lpmod(t46:t60)),lpdat/lpscale, ...
         mplmod/mean(mplmod(t46:t60)),mpldat/r2lsc,texp(:,6:8)];
disp('Requesting information about the fixed point ...')
disp(' ')
disp(' ==================================================================')
disp('                                       Measures of Compensation    ')
disp('       Nonmilitary                  -------------------------------')
disp('          Hours         MP Labor    Nonfarm  Manufacturing Civilian')
disp('      -------------- -------------- -------------------------------')
disp(' Year  Model  Data    Model  Data      Data       Data       Data  ')
disp(sprintf(' %4g%7.4f%7.4f%8.4f%7.4f%10.4f%11.4f%11.4f\n',texp(1:8,:)'));
disp(' ==================================================================')
disp(' NOTE: all scaling done using data averages for 1946-1960')
disp(' ')
disp(' ==================================================================')
disp('                                       Measures of Compensation    ')
disp('       Nonmilitary                  -------------------------------')
disp('          Hours         MP Labor    Nonfarm  Manufacturing Civilian')
disp('      -------------- -------------- -------------------------------')
disp(' Year  Model  Data    Model  Data      Data       Data       Data  ')
disp(sprintf(' %4g%7.4f%7.4f%8.4f%7.4f%10.4f%11.4f%11.4f\n',texp2(1:8,:)'));
disp(' ==================================================================')
disp(' NOTE: Hours and MP Labor scaled by own averages for 1946-1960')
end;

if (iskip==1); 
  nfig=1; 
else;
  nfig=8;
end;
figure(nfig)
subplot(321)
plot(t,gdat,t,cgdat,t,igdat), axis([tt0,tt1,0,1])
title('INPUTS:Spending')
subplot(322)
plot(t,adat*100),axis([tt0,tt1,0,12])
title('Draft')
subplot(323)
plot(t,tauldat*100),axis([tt0,tt1,5,25])
title('Labor tax')
subplot(324)
plot(t,taukdat*100),axis([tt0,tt1,40,70])
title('Capital tax')
subplot(325)
plot(t,z1dat_bea.^(1-theta)),axis([tt0,tt1,1.25,1.75])
title('TFP--Bench')
subplot(326)
plot(t,z2dat_bea.^(1-theta)),axis([tt0,tt1,1.6,2.25])
title('TFP--CU Model')

figure(nfig+1)
subplot(321)
%plot(t,gnpdat,'r-x',t,gnpmod,'b-o'),axis([tt0,tt1,1,1.5])
plot(t,gnpdat,'r-x',t,gnpmod,'b-o'),axis([tt0,tt1,.8,1.5])
title('OUTPUTS: GNP')
legend('Data=Red','Model=Blue')
subplot(322) 
plot(t,cpdat,'r-x',t,cpmod,'b-o'),axis([tt0,tt1,.45,.65,])
title('Consumption')
subplot(323)
plot(t,ipdat,'r-x',t,max(ipmod,0),'b-o'),axis([tt0,tt1,0,.4])
title('Investment')
subplot(324)
plot(t,lpdat/lpscale,'r-x',t,lpmod/lpscale,'b-o'),axis([tt0,tt1,.9,1.4])
title('Nonmilitary hours')
subplot(325)
plot(t,atr2kdat,'r-x',t,atr2kmod,'b-o'),axis([tt0,tt1,2.5,4.5])
title('After-tax return to capital (%)')
subplot(326)
plot(t,atr2ldat/atr2lsc,'r-x',t,atr2lmod/atr2lsc,'b-o'),axis([tt0,tt1,.92,1.06])
title('After-tax return to labor')


inputs1=[t,gdat,cgdat,igdat,adat,taukdat,tauldat,z1dat_bea];
inputs2=[t,gdat,cgdat,igdat,adat,taukdat,tauldat,z2dat_bea];
outputs=[t,gnpdat,gnpmod,cpdat,cpmod,ipdat,max(ipmod,0), ...
           lpdat/lpscale,lpmod/lpscale, ...
           atr2kdat,atr2kmod,atr2ldat/atr2lsc,atr2lmod/atr2lsc];
           %lprddat,lprdmod];
