% -------------------------------------------------------------------------
%
% Project: Uneven Growth
%
% Purpose: Create Figure 9
%
% -------------------------------------------------------------------------

clear variables
close all

cd('C:/Users/Pascual Restrepo/Dropbox/Research/wealth_inequality/AutomationInequality/replication_mrr_january_2021/additional_evidence/investment_data')
addpath 'C:/Users/Pascual Restrepo/Documents/MATLAB/scripts'
addpath 'raw_data'
addpath 'output'

% -------------------------------------------------------------------------
%% - Preamble -
% -------------------------------------------------------------------------

% Set latex as default interpreter
set(groot, 'defaulttextinterpreter','latex'); 
set(groot, 'defaultAxesTickLabelInterpreter','latex'); 
set(groot, 'defaultLegendInterpreter','latex');

% -------------------------------------------------------------------------
%% - Figure 9 -
% -------------------------------------------------------------------------

% - Get data --------------------------------------------------------------

data_KY = csvread('bls_capital_series.csv',1,0); %Notes: Series from Gourio and Klein
 
% Private sector capital to output ratio % changes since 1950
KY_year              = data_KY(:,1);
KY_BLS_PS            = data_KY(:,2);
KY_BLS_PS_nonfarm    = data_KY(:,3);

%Nipa: investment
data_pnfi_value = csvread('nipa-t5-2-5-millions_investment.csv',6,2, [6 2 29 90]);
pnfi_value      = data_pnfi_value(10,:); %Table 5.2.5, Line 10: Private non-residential fixed investment.

data_pnfi_qty   = csvread('nipa-t5-2-3-quantity_investment.csv',6,2, [6 2 29 90]);
pnfi_qty        = data_pnfi_qty(10,:); %Table 5.2.3, Line 10: Private non-residential fixed investment.

%NIPA: gross domestic product
data_gdp_value = csvread('nipa-t1-1-5-millions_gdp.csv',6,2, [6 2 18 90]);
gdp_value      = data_gdp_value(1,:); %Table 1.1.5, Line 1: Gross domestic product.
pnfi2_value    = data_gdp_value(9,:); %Table 1.1.5, Line 9: Private non-residential fixed investment / should equal series above.

data_gdp_qty   = csvread('nipa-t1-1-3-quantity_gdp.csv',6,2, [6 2 18 90]);
gdp_qty        = data_gdp_qty(1,:); %Table 1.1.3, Line 1: Gross domestic product.
pnfi2_qty      = data_gdp_qty(9,:); %Table 1.1.3, Line 9: Private non-residential fixed investment / should equal series above.

assert(abs(max(pnfi_qty-pnfi2_qty))<0.001);
assert(abs(max(pnfi_value-pnfi2_value))<0.001);

%FIXED ASSET TABLES: capital net of depreciation
data_pnfk_value = csvread('fixed-assets-t4-1-billions_nonresidential-assets.csv',6,6, [6 6 9 94]);
data_pnfd_value = csvread('fixed-assets-t4-4-billions_nonresidential-depreciation.csv',6,6, [6 6 9 94]);
pnfk_value      = data_pnfk_value(1,:)*1000+data_pnfd_value(1,:)*1000; %Table 4.1, Line 1 + Table 4.4, Line 1: (gross) Private nonresidential fixed assets

data_pnfk_qty   = csvread('fixed-assets-t4-2-quantity_nonresidential-assets.csv',6,6, [6 6 9 94]);
data_pnfd_qty   = csvread('fixed-assets-t4-5-quantity_nonresidential-depreciation.csv',6,6, [6 6 9 94]);

%Compute chain quantity index of gross capital

%Get price indices:
knet_qty        = data_pnfk_qty(1,:);
delta_qty       = data_pnfd_qty(1,:);
knet_price      = data_pnfk_value(1,:)./data_pnfk_qty(1,:);
delta_price     = data_pnfd_value(1,:)./data_pnfd_qty(1,:);

growth_rate     = sqrt( (knet_price(2:end).*knet_qty(2:end)+delta_price(2:end).*delta_qty(2:end))./(knet_price(2:end).*knet_qty(1:end-1)+delta_price(2:end).*delta_qty(1:end-1)) ...
                       .*(knet_price(1:end-1).*knet_qty(2:end)+delta_price(1:end-1).*delta_qty(2:end))./(knet_price(1:end-1).*knet_qty(1:end-1)+delta_price(1:end-1).*delta_qty(1:end-1)) ...
                       );

pnfk_qty        = pnfk_value(1)*exp(cumsum([0, log(growth_rate)]));

% Compute relative growth since 1950 for all series
nipa_year=1950:1:2017;

rel_pnfi_value = 100*(log(pnfi_value(22:end)/pnfi_value(22))-log(gdp_value(22:end)/gdp_value(22)));
rel_pnfk_value = 100*(log(pnfk_value(22:end)/pnfk_value(22))-log(gdp_value(22:end)/gdp_value(22)));

rel_pnfi_qty   = 100*(log(pnfi_qty(22:end)/pnfi_qty(22))-log(gdp_qty(22:end)/gdp_qty(22)));
rel_pnfk_qty   = 100*(log(pnfk_qty(22:end)/pnfk_qty(22))-log(gdp_qty(22:end)/gdp_qty(22)));

% - Plot figures --------------------------------------------------------------

figure_ky = figure('units','centimeters','position',[10 10 7*(16/9)*2 10]);
 
subplot(1,2,1)
hold on
plot(KY_year,KY_BLS_PS, 'k-', 'LineWidth',1.2);
message = sprintf('private \nsector (BLS)');
text(2017,1.05*KY_BLS_PS(end),message, 'color', 'k');

plot(KY_year,KY_BLS_PS_nonfarm, 'b-', 'LineWidth', 1.2);
message = sprintf('nonfarm private \nsector (BLS)');
text(2017,KY_BLS_PS_nonfarm(end),message, 'color', 'b');

plot(nipa_year,rel_pnfk_value, 'r-', 'LineWidth', 1.2);
message = sprintf('private fixed \nassets (BEA)');
text(2017,0.85*rel_pnfk_value(end),message,'color', 'r');

ylabel('percent change relative to 1950')
title('Panel A. capital-output ratios')
xlim([1950, 2029])
set(gca,'TickLabelInterpreter','latex')
set(gcf, 'Color', 'w');
 
 
subplot(1,2,2)
hold on
plot(nipa_year, rel_pnfi_qty, 'k-',  'LineWidth', 1.2)
message = sprintf('quantity \ninvestment to \nquantity output \n(NIPA)');
text(2018,0.92*rel_pnfi_qty(end),message, 'color', 'k');

plot(nipa_year, rel_pnfi_value, 'r-', 'LineWidth', 1.2)
message = sprintf('value \ninvestment to \nvalue of \noutput (NIPA)');
text(2018,0.9*rel_pnfi_value(end),message,'color', 'r');

ylabel('percent change relative to 1950')
title('Panel B. investment-output ratios')
xlim([1950, 2029])
 
export_fig output/Figure_paper_investment.pdf
