
clear all;
close all;
clc;

% Add path of data and functions
addpath('functions');
addpath('data');


p = 3;              % p is number of lags in the VAR part

% Impulse response setting                               
nhor = 40;          % Identification horizon
ident = 3;          % 1: Long+Short; 2: Barsky-Sims (2011), : baseline Max Share identification 4: Kurmann Sims identification (set nhor=80)
nirf = 40;          % Impulse response horizon
ndraws = 500;       % Number of draws for simulation

stand_data = 0;       % 0: leave data untouched, 1: standardize data

plot_news_ser = 0;    % Plot estimated TFP News Shock time series
plot_fevd = 0;        % Plot FEVD for all variables
plot_fevd_tfp = 0;    % Plot FEVD only for TFP
plot_news_irf = 1;    % Plot TFP News Shock IRFs
plot_surpr_irf = 0;   % Plot TFP Surprise shock IRFs
plot_confidence = 0;  % Plot Confidence shock IRFs
plot_cross_corr = 0;  % Plot Cross Correlations
plot_counterf = 0;    % Plot counterfactual impulse responses
monte_carlo = 0;      % Do Monte carlo on model data

%----------------------------- END OF PRELIMINARIES --------------------------------
tic;
%----------------------------------LOAD DATA----------------------------------------   



% Inventory accumulation paper
Ydata = xlsread('data_invt.xlsx','quarterly_long','DK160:FW280'); % baseline horizon (280 for ECC/DCC)
% Ydata = xlsread('data_invt_updateSep18.xlsx','quarterly_long','DU160:GP293'); % longer horizon for appendix version of figure 1
% Ydata = xlsread('data_invt.xlsx','quarterly_long','DK188:FS280'); % when using buiness inventories
Ynames = {'TFP'; 'TFP_I'; 'TFP_C'; 'GDP'; 'Investment'; 'Consumption'; 'Hours'; 'BAA'; 'E5Y - Confidence Indicator'; 'SP500';...
    'UMSENT'; 'USSLIND'; 'TB3MS'; 'Business Inventories'; 'Inventories: Retail Sector'; 'Inventories: Wholesale Sector'; 'Inventories: Manufacturing Sector'; 'Finished Goods Inventories'; 'Work in Process Inventories'; 'Input Inventories';...
    'Inventory to Sales_alternative'; 'IS_r'; 'IS_w'; 'IS_m'; 'invtp'; 'Private Non-Farm Inventories'; 'Invtp_m'; 'Invtp_r'; 'Invtp_w'; 'FFR';...
    'Inflation'; 'Captial Utilization'; 'EBP'; 'RP'; 'DR'; 'IBES'; 'Inventory to Sales'; 'ICC3'; 'Implied Cost of Capital: GORDON'; 'Implied Cost of Capital: OJ';...
    'Implied Cost of Capital: MPEG'; 'Implied Cost of Capital: CT'; 'Implied Cost of Capital: GLS'; 'ICC5'; 'Equity Cost of Capital'; 'Debt Cost of Capital'; 'LPROD'; 'RRF'; 'INVTPaS'; 'MC1';...
    'Marginal Cost CES: Nekarda-Ramey 1'; 'Marginal Cost CES: Gali et al.'; 'Marginal Cost CES: Nekarda-Ramey 2'; 'MC_CD1'; 'Marginal Cost CD: Nekarda-Ramey 1'; 'Marginal Cost CD: Gali et al.'; 'Marginal Cost CD: Nekarda-Ramey 2'; 'Real Wage'; 'ECC';...
    'DCC'; 'Inventory-to-Sales'; 'BAA yield'; 'Real Rate'; 'Inflation Expectations'; 'Real BAA Bond Yield';
};

    % private non-farm inventories is 'INVTpa'
    % business inventories is INVT
    % business inventories: retail is INVT_r
    % business inventories: wholesale is INVT_w
    % business inventories: manufacturing is INVT_m
    % inventories in manufacturing sector: Finished Goods Inventories is invt_mf
    % inventories in manufacturing sector: Work in Process Inventories is invt_mwip
    % inventories in manufacturing sector: Input Inventories is invt_mrm
    % equity cost of capital is ECC
    % debt cost of capital is DCC
    % marginal cost CES - Nekarda Ramey 1 is MC2
    % Marginal Cost CES: Gali et al. is MC3
    % Marginal Cost CES: Nekarda-Ramey 2 is MC_PW4
    % Marginal Cost CD - Nekarda Ramey 1 is MC_CD2
    % Marginal Cost CD - Gali et al. is MC_CD3
    % Marginal Cost CD - Nekarda Ramey 2 is MC_CD_PW4
    

savefig = 1; % figure saved if ==1

% Figures in paper:

% Figure 1 --- baseline sample of line 38
% model = [1 4 5 6 7 26 10];   % figure 1 (also for surprise shock in appendix)
% figurename = 'May19_85Q1_15Q1_fig1';




% Figure with inventories by sector --- business inventories: change sample in line 38 and 40
% model = [1 4 5 6 7 15 10]; 
% figurename = 'May19_92Q1_15Q1_businvt_retail';

% model = [1 4 5 6 7 16 10];
% figurename = 'May19_92Q1_15Q1_businvt_wholes';

% model = [1 4 5 6 7 17 10];
% figurename = 'May19_92Q1_15Q1_businvt_manuf';




% figure inventory types --- business inventories: change sample in line 38 and 40
% model = [1 4 5 6 7 18 10]; 
% figurename = 'May19_92Q1_15Q1_businvt_manuf_fg';

% model = [1 4 5 6 7 20 10];
% figurename = 'May19_92Q1_15Q1_businvt_manuf_rm';




% Figure with ECC and DCC measures --- baseline sample of line 38
% model = [1 4 45 6 7 26 10];  % figure ECC+DCC
% figurename = 'May19_85Q1_15Q1_ecc';

% model = [1 4 46 6 7 26 10];  % figure ECC+DCC
% figurename = 'May19_85Q1_15Q1_dcc';




% Figure with ICC measures --- baseline sample of line 38
% model = [1 4 39 6 7 26 10];  % figure ICC: Gordon
% figurename = 'May19_85Q1_15Q1_iccGordon';

% model = [1 4 40 6 7 26 10];  % figure ICC: OJ
% figurename = 'May19_85Q1_15Q1_iccOJ';

% model = [1 4 41 6 7 26 10];  % figure ICC: MPEG
% figurename = 'May19_85Q1_15Q1_iccMPEG';

% model = [1 4 43 6 7 26 10];  % figure ICC: GLS
% figurename = 'May19_85Q1_15Q1_iccGLS';



% Figure with marginal cost measures --- baseline sample of line 38
% model = [1 4 55 6 7 26 10];  % figure MC: CD NR1
% figurename = 'May19_85Q1_15Q1_mcCD_NR1';

% model = [1 4 51 6 7 26 10];  % figure MC: CES NR1
% figurename = 'May19_85Q1_15Q1_mcCES_NR1';

% model = [1 4 57 6 7 26 10];  % figure MC: CD NR2
% figurename = 'May19_85Q1_15Q1_mcCD_NR2';

% model = [1 4 53 6 7 26 10];  % figure MC: CES NR2
% figurename = 'May19_85Q1_15Q1_mcCES_NR2';




% Figure with real baa bond yield --- baseline sample of line 38
% model = [1 4 65 6 7 26 10]; % figure with real BAA bond yield
% figurename = 'May19_85Q1_15Q1_baa';




%% CODE
ydata = [Ydata(:,model) ];
names = Ynames(model);

% Check if all observations are in order (some series might begin later)
f = zeros(size(ydata,2),1);
for i = 1:size(ydata,2)
    f(i,:) = min(find(1-isnan(ydata(:,i))));
end
chop = max(f);
ydata = ydata(chop:end,:);
Ydata = Ydata(chop:end,:);

%Standardize data
if stand_data == 1
    Y = standardize(ydata);
else    
    Y = ydata;
end

% Number of observations and dimension of X and Y
t = size(Y,1); % t is the time-series observations of Y
M = size(Y,2); % M is the dimensionality of Y
t2 = size(Y,1);

% Estimate VAR in full sample, calculate IRFs, FEVDs and NEWS shock
disp(['Running VAR on ' num2str(t2) ' observations....'])
[impulses_news2,D2,ND2, news1, surp1, impulses_sur2] = bvarnews(Y,p,ndraws,nhor,ident);
disp('Success!')

% Get quantiles of news shocks
newshock2 = squeeze(quantile(impulses_news2,.5,3));   
newshocklowalt2 = squeeze(quantile(impulses_news2,.16,3));
newshockhighalt2 = squeeze(quantile(impulses_news2,.84,3));
surshock2 = squeeze(quantile(impulses_sur2,.5,3));   
surshocklowalt2 = squeeze(quantile(impulses_sur2,.16,3));
surshockhighalt2 = squeeze(quantile(impulses_sur2,.84,3));
D_F2 = squeeze(quantile(D2,[0.16, 0.5, 0.84]));
N_F2 = squeeze(quantile(ND2,[0.16, 0.5, 0.84]));
news1r = squeeze(quantile(news1,.5,2));
surp1r = squeeze(quantile(surp1,.5,2));

% Plot figures

figure1 = figure('Color',[1 1 1]);
set(gcf,'units','normalized','outerposition',[0 0 1 1]);
for i=1:(M)     
    subplot(round((M)/2),3,i)  
    shadedplot(1:nirf,100*newshocklowalt2(i,1:nirf),100*newshockhighalt2(i,1:nirf),[0.4,0.4,0.4],[0.4,0.4,0.4]);
    hold all
    plot(squeeze(100*newshock2(i,1:nirf))','k','Linewidth',2);
     %hold on
     %plot(squeeze(100*newshock2mfev(i,:))','--k','Linewidth',2)
%      hold on
%      plot(squeeze(100*newshockiMF(i,:))','-.k','Linewidth',3)
    xlim([1 nirf])
    grid on
    xlabel(['quarters'])
    title([cell2mat(names(i))], 'FontSize', 16)
end

if savefig ==1 
    % save figure
    set(gcf, 'Renderer', 'opengl')
    set(figure1,'Units','Inches');
    pos = get(figure1,'Position');
    set(figure1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
    fname = ['Fig_' figurename '_news'];
    %print(figure1,fname,'-dpdf','-r0')
    % save data
    newshock_median = 100*newshock2;
    newshock_high = 100*newshockhighalt2;
    newshock_low = 100*newshocklowalt2;
    irf.newshock_median = newshock_median;
    irf.newshock_high = newshock_high;
    irf.newshock_low = newshock_low;
    irf.names = names;
    irf.surp_median = 100*surshock2;
    irf.surp_high = 100*surshockhighalt2;
    irf.surp_low = 100*surshocklowalt2;
    eval(['save ' figurename '_MaxSnews_irf irf N_F2'])
    %eval(['save ' figurename '_KS80news_irf irf N_F2'])
end



figure2 = figure('Color',[1 1 1]);  
set(gcf,'units','normalized','outerposition',[0.5 0 .5 1]);
for i=1:M
    subplot(round((M)/2),3,i)
    shadedplot(1:nirf,squeeze(N_F2(1,i,1:nirf))',squeeze(N_F2(3,i,1:nirf))',[0.4,0.4,0.4],[0.4,0.4,0.4]);
    hold all
    plot(squeeze(N_F2(2,i,1:nirf)),'black','Linewidth',2);
%     hold on
%     plot(squeeze(N_F2c(1,i,:)), 'k-.','Linewidth',2)
%     hold on
%     plot(squeeze(N_F2i(1,i,:)),'k--','Linewidth',3)
    xlim([1 nirf])
    grid on
%     title(['%fevd of news shock on the forecast error of ' cell2mat(names(i))])
    title(['FEVD: ' cell2mat(names(i))], 'FontSize', 16)
end
%saveas(gcf,'fevd_second_sample','fig')
close(figure2)

figure3 = figure('Color',[1 1 1]);
set(gcf,'units','normalized','outerposition',[0 0 1 1]);
for i=1:(M)     
    subplot(round((M)/2),3,i)  
    shadedplot(1:nirf,100*surshocklowalt2(i,1:nirf),100*surshockhighalt2(i,1:nirf),[0.4,0.4,0.4],[0.4,0.4,0.4]);
    hold all
    plot(squeeze(100*surshock2(i,1:nirf))','k','Linewidth',2);
     %hold on
     %plot(squeeze(100*newshock2mfev(i,:))','--k','Linewidth',2)
%      hold on
%      plot(squeeze(100*newshockiMF(i,:))','-.k','Linewidth',3)
    xlim([1 nirf])
    grid on
    xlabel(['quarters'])
    title([cell2mat(names(i))], 'FontSize', 16)
end

if savefig ==1 
    % save figure
    set(gcf, 'Renderer', 'opengl')
    set(figure3,'Units','Inches');
    pos = get(figure3,'Position');
    set(figure3,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
    fname = ['Fig_' figurename '_surp'];
    %print(figure3,fname,'-dpdf','-r0')
end
close(figure3)
disp('Finito')

