clear;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: July 10, 2022
%
% Program: This program integrates the analysis of the aggregates and
% and their standard deviations
% 1. Compute the aggregate growth and aggregate variance decomposition for
% export sales
% 2. Compute analytical standard deviations for the aggregate variance and its compo-
% nents using
% 3. Prepare the data for figure plots (Aggregate variance and components with standard errors) 
% 4. Export output files in a proper dataset
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%% Destination effects yearly frequency

% This part reads in the baseline regression micro results with destination effects only and a yearly frequency.
% Get raw data 
data=csvread('Output/growth_regs_capays_yearly.csv',1,0);

id_fp = data(:,1);
id_f = data(:,2);
id_p = data(:,3);
year = data(:,4);
id_py = data(:,5);
gr_variable = data(:,6);
gr_idio = data(:,7);
gr_idio_com = data(:,8);
gr_idio_dest = data(:,9);
gr_macro = data(:,10);
w_fp = data(:,11);
w_f = data(:,12);
w_p = data(:,13);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: July 10, 2022
%
% The following portion of code
%
%   1) create aggregate variance and components using constant weights
%   2) create aggregate vairance and components using time-varying weights
%   3) compute Herfindahl indices based on the distribution of weights
% 
% References to the following functions: 
%   ag_growth.m (aggregate growth and its variance, based on time-varying 
%       weights)
%   var_growth.m (variance of aggregate growth, based on constant weights
%   var_only_growth.m (variance of aggregate growth, ignoring within 
%       components covariance terms)
%   herfindahl.m (Herfindahl index)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


wgr_variable = w_fp.* gr_variable;
[agr_variable sd_agvariable] = ag_growth(id_fp,year,wgr_variable);
wgr_idio = w_fp.* gr_idio;
[agr_idio sd_agidio] = ag_growth(id_fp,year,wgr_idio);
wgr_macro = w_fp.* gr_macro;
[agr_macro sd_agmacro] = ag_growth(id_fp,year,wgr_macro);

% Create aggregate variance components: time-varying weights

[VAR_variable var_variable sd_variable] = var_growth(id_fp,year,gr_variable,w_fp);
[VAR_idio var_idio sd_idio] = var_growth(id_fp,year,gr_idio,w_fp);
[VAR_idio_dest var_idio_dest sd_idio_dest] = var_growth(id_fp,year,gr_idio_dest,w_fp);
[VAR_idio_com var_idio_com sd_idio_com] = var_growth(id_fp,year,gr_idio_com,w_fp);
[VAR_macro var_macro sd_macro] = var_growth(id_fp,year,gr_macro,w_fp);


% Herfindahl index

herf_fp = herfindahl(id_fp,year,w_fp);
herf_f = herfindahl(id_f,year,w_f);
herf_p = herfindahl(id_p,year,w_p);

t = linspace(min(year),max(year),(max(year)-min(year)+1))';


% Results in Table 5 (Panel III)
% Standard Deviations
actual_vol = mean(sd_variable)
firm_vol = mean(sd_idio)
macro_vol = mean(sd_macro)
% Mean growth rates
actual_mean_gr = mean(agr_variable)
firm_mean_gr = mean(agr_idio)
macro_mean_gr = mean(agr_macro)

statistics_table = table(actual_vol, firm_vol, macro_vol,...
                         actual_mean_gr, firm_mean_gr, macro_mean_gr);



% Calculate variances/covariances

gr_sales = gr_variable;
gr_mac = gr_macro;

[VAR_sales,var_sales,sd_sales] = var_growth(id_fp,year,gr_sales,w_fp);
[VAR_idio,var_idio,sd_idio] = var_growth(id_fp,year,gr_idio,w_fp);
[VAR_mac,var_mac,sd_mac] = var_growth(id_fp,year,gr_mac,w_fp);

% Calculate standard errors

[var_sales_se,sd_sales_se] = std_err(id_fp,year,gr_sales,w_fp,var_sales);
[var_idio_se,sd_idio_se] = std_err(id_fp,year,gr_idio,w_fp,var_idio);
[var_mac_se,sd_mac_se] = std_err(id_fp,year,gr_mac,w_fp,var_mac);


% Plotting

export_vol = sd_sales;
export_vol_se_max = sd_sales+1.96*sd_sales_se;
export_vol_se_min = sd_sales-1.96*sd_sales_se;

export_idio_vol = sd_idio;
export_idio_vol_se_max = sd_idio+1.96*sd_idio_se;
export_idio_vol_se_min = sd_idio-1.96*sd_idio_se;

export_macro_vol = sd_mac;
export_macro_vol_se_max = sd_idio+1.96*sd_mac_se;
export_macro_vol_se_min = sd_idio-1.96*sd_mac_se;

year = t;

res_table = table(year,...
                 export_vol,export_vol_se_max,export_vol_se_min,...
                 export_idio_vol,export_idio_vol_se_max,export_idio_vol_se_min,...
                 export_macro_vol,export_macro_vol_se_max,export_macro_vol_se_min);


writetable(res_table,'Output/ts_outputs_capays_yearly.csv')
writetable(statistics_table,'Output/aggregate_outputs_capays_yearly.csv')

%%%%%%%%%%%%%%%%%% Destination and sector effects yearly frequency

clear;
close all;

% This part reads in the baseline regression micro results with destination and sector effects
% and a yearly frequency.
% Get raw data 
data=csvread('Output/growth_regs_capays_destsec_yearly.csv',1,0);

id_fp = data(:,1);
id_f = data(:,2);
id_ip = data(:,3);
id_i = data(:,4); 
id_p = data(:,5);
year = data(:,6);
id_ipy = data(:,7); 
gr_variable = data(:,8);
gr_idio = data(:,9);
gr_idio_com = data(:,10);
gr_idio_dest = data(:,11);
gr_sec = data(:,12);
gr_macro = data(:,13);
w_fp = data(:,14); 
w_f = data(:,15);
w_ip = data(:,16);
w_i = data(:,17);
w_p = data(:,18);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: July 10, 2022
%
% The following portion of code
%
%   1) create aggregate variance and components using constant weights
%   2) create aggregate vairance and components using time-varying weights
%   3) compute Herfindahl indices based on the distribution of weights
% 
% References to the following functions: 
%   ag_growth.m (aggregate growth and its variance, based on time-varying 
%       weights)
%   var_growth.m (variance of aggregate growth, based on constant weights
%   var_only_growth.m (variance of aggregate growth, ignoring within 
%       components covariance terms)
%   herfindahl.m (Herfindahl index)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


gr_macrosec = gr_macro+gr_sec;

% Create aggregate growth rates and components and compute standard
% deviations
wgr_variable = w_fp.*gr_variable;
[agr_variable sd_agvariable] = ag_growth(id_fp,year,wgr_variable);
wgr_idio = w_fp.*gr_idio;
[agr_idio sd_agidio] = ag_growth(id_fp,year,wgr_idio);
wgr_macrosec = w_fp.*gr_macrosec;
[agr_macrosec sd_agmacrosec] = ag_growth(id_fp,year,wgr_macrosec);

% Create aggregate variance components: time-varying weights

[VAR_variable var_variable sd_variable] = var_growth(id_fp,year,gr_variable,w_fp);
[VAR_idio var_idio sd_idio] = var_growth(id_fp,year,gr_idio,w_fp);
[VAR_idio_dest var_idio_dest sd_idio_dest] = var_growth(id_fp,year,gr_idio_dest,w_fp);
[VAR_idio_com var_idio_com sd_idio_com] = var_growth(id_fp,year,gr_idio_com,w_fp);
[VAR_macrosec var_macrosec sd_macrosec] = var_growth(id_fp,year,gr_macrosec,w_fp);

% Herfindahl index

herf_fp = herfindahl(id_fp,year,w_fp);
herf_f = herfindahl(id_f,year,w_f);
herf_ip = herfindahl(id_ip,year,w_ip);
herf_i = herfindahl(id_i,year,w_i);
herf_p = herfindahl(id_p,year,w_p);

t = linspace(min(year),max(year),(max(year)-min(year)+1))';

% Results in Table 5 (Panel III)
% Standard Deviations
actual_vol = mean(sd_variable)
firm_vol = mean(sd_idio)
macrosec_vol = mean(sd_macrosec)
% Mean growth rates
actual_mean_gr = mean(agr_variable)
firm_mean_gr = mean(agr_idio)
macrosec_mean_gr = mean(agr_macrosec)

statistics_table = table(actual_vol, firm_vol, macrosec_vol,...
                         actual_mean_gr, firm_mean_gr, macrosec_mean_gr);



% Calculate variances/covariances

gr_sales = gr_variable;
gr_macsec = gr_macro+gr_sec;

% Calculate variances/covariances

[VAR_sales,var_sales,sd_sales] = var_growth(id_fp,year,gr_sales,w_fp);
[VAR_idio,var_idio,sd_idio] = var_growth(id_fp,year,gr_idio,w_fp);
[VAR_macsec,var_macsec,sd_macsec] = var_growth(id_fp,year,gr_macsec,w_fp);

% Calculate standard errors

[var_sales_se,sd_sales_se] = std_err(id_fp,year,gr_sales,w_fp,var_sales);
[var_idio_se,sd_idio_se] = std_err(id_fp,year,gr_idio,w_fp,var_idio);
[var_macsec_se,sd_macsec_se] = std_err(id_fp,year,gr_macsec,w_fp,var_macsec);


% Plotting

export_vol = sd_sales;
export_vol_se_max = sd_sales+1.96*sd_sales_se;
export_vol_se_min = sd_sales-1.96*sd_sales_se;

export_idio_vol = sd_idio;
export_idio_vol_se_max = sd_idio+1.96*sd_idio_se;
export_idio_vol_se_min = sd_idio-1.96*sd_idio_se;

export_macrosec_vol = sd_macsec;
export_macrosec_vol_se_max = sd_macsec+1.96*sd_macsec_se;
export_macrosec_vol_se_min = sd_macsec-1.96*sd_macsec_se;

year = t;

res_table = table(year,...
                 export_vol,export_vol_se_max,export_vol_se_min,...
                 export_idio_vol,export_idio_vol_se_max,export_idio_vol_se_min,...
                 export_macrosec_vol,export_macrosec_vol_se_max,export_macrosec_vol_se_min);


writetable(res_table,'Output/ts_outputs_capays_destsec_yearly.csv')
writetable(statistics_table,'Output/aggregate_outputs_capays_destsec_yearly.csv')



%%%%%%%%%%%%%%%%%% Destination effects quarterly frequency


clear;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: July 10, 2022
%
% Program: This program integrates the analysis of the aggregates and
% and their standard deviations
% 1. Compute the aggregate growth and aggregate variance decomposition for
% export sales
% 2. Compute analytical standard deviations for the aggregate variance and its compo-
% nents using
% 3. Prepare the data for figure plots (Aggregate variance and components with standard errors) 
% 4. Export output files in a proper dataset
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% This part reads in the baseline regression micro results with destination effects only and a yearly frequency.
% Get raw data 
data=csvread('Output/growth_regs_capays_quarterly.csv',1,0);


id_fp = data(:,1);
id_f = data(:,2);
id_p = data(:,3);
period = data(:,4);
id_py = data(:,5);
gr_variable = data(:,6);
gr_idio = data(:,7);
gr_idio_com = data(:,8);
gr_idio_dest = data(:,9);
gr_macro = data(:,10);
w_fp = data(:,11);
w_f = data(:,12);
w_p = data(:,13);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: July 10, 2022
%
% The following portion of code
%
%   1) create aggregate variance and components using constant weights
%   2) create aggregate vairance and components using time-varying weights
%   3) compute Herfindahl indices based on the distribution of weights
% 
% References to the following functions: 
%   ag_growth.m (aggregate growth and its variance, based on time-varying 
%       weights)
%   var_growth.m (variance of aggregate growth, based on constant weights
%   var_only_growth.m (variance of aggregate growth, ignoring within 
%       components covariance terms)
%   herfindahl.m (Herfindahl index)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


wgr_variable = w_fp.* gr_variable;
[agr_variable sd_agvariable] = ag_growth(id_fp,period,wgr_variable);
wgr_idio = w_fp.* gr_idio;
[agr_idio sd_agidio] = ag_growth(id_fp,period,wgr_idio);
wgr_macro = w_fp.* gr_macro;
[agr_macro sd_agmacro] = ag_growth(id_fp,period,wgr_macro);

% Create aggregate variance components: time-varying weights

[VAR_variable var_variable sd_variable] = var_growth(id_fp,period,gr_variable,w_fp);
[VAR_idio var_idio sd_idio] = var_growth(id_fp,period,gr_idio,w_fp);
[VAR_idio_dest var_idio_dest sd_idio_dest] = var_growth(id_fp,period,gr_idio_dest,w_fp);
[VAR_idio_com var_idio_com sd_idio_com] = var_growth(id_fp,period,gr_idio_com,w_fp);
[VAR_macro var_macro sd_macro] = var_growth(id_fp,period,gr_macro,w_fp);


% Herfindahl index

herf_fp = herfindahl(id_fp,period,w_fp);
herf_f = herfindahl(id_f,period,w_f);
herf_p = herfindahl(id_p,period,w_p);

t = linspace(min(period),max(period),((max(period)-min(period))/25)+1)';


% Results in Table 5 (Panel III)
% Standard Deviations
actual_vol = mean(sd_variable)
firm_vol = mean(sd_idio)
macro_vol = mean(sd_macro)
% Mean growth rates
actual_mean_gr = mean(agr_variable)
firm_mean_gr = mean(agr_idio)
macro_mean_gr = mean(agr_macro)

statistics_table = table(actual_vol, firm_vol, macro_vol,...
                         actual_mean_gr, firm_mean_gr, macro_mean_gr);



% Calculate variances/covariances

gr_sales = gr_variable;
gr_mac = gr_macro;

[VAR_sales,var_sales,sd_sales] = var_growth(id_fp,period,gr_sales,w_fp);
[VAR_idio,var_idio,sd_idio] = var_growth(id_fp,period,gr_idio,w_fp);
[VAR_mac,var_mac,sd_mac] = var_growth(id_fp,period,gr_mac,w_fp);

% Calculate standard errors

[var_sales_se,sd_sales_se] = std_err(id_fp,period,gr_sales,w_fp,var_sales);
[var_idio_se,sd_idio_se] = std_err(id_fp,period,gr_idio,w_fp,var_idio);
[var_mac_se,sd_mac_se] = std_err(id_fp,period,gr_mac,w_fp,var_mac);


% Plotting

export_vol = sd_sales;
export_vol_se_max = sd_sales+1.96*sd_sales_se;
export_vol_se_min = sd_sales-1.96*sd_sales_se;

export_idio_vol = sd_idio;
export_idio_vol_se_max = sd_idio+1.96*sd_idio_se;
export_idio_vol_se_min = sd_idio-1.96*sd_idio_se;

export_macro_vol = sd_mac;
export_macro_vol_se_max = sd_idio+1.96*sd_mac_se;
export_macro_vol_se_min = sd_idio-1.96*sd_mac_se;

year = floor(t/100);
quarter = ((t/100-floor(t/100))*100)/25 + 1;

res_table = table(year,quarter,...
                 export_vol,export_vol_se_max,export_vol_se_min,...
                 export_idio_vol,export_idio_vol_se_max,export_idio_vol_se_min,...
                 export_macro_vol,export_macro_vol_se_max,export_macro_vol_se_min);


writetable(res_table,'Output/ts_outputs_capays_quarterly.csv')
writetable(statistics_table,'Output/aggregate_outputs_capays_quarterly.csv')



%%%%%%%%%%%%%%%%%% Destination and sector effects quarterly frequency


clear;
close all;

% This part reads in the baseline regression micro results with destination and sector effects
% and a yearly frequency.
% Get raw data 
data=csvread('Output/growth_regs_capays_destsec_quarterly.csv',1,0);

    
id_fp = data(:,1);
id_f = data(:,2);
id_ip = data(:,3);
id_i = data(:,4); 
id_p = data(:,5);
period = data(:,6);
id_ipq = data(:,7); 
gr_variable = data(:,8);
gr_idio = data(:,9);
gr_idio_com = data(:,10);
gr_idio_dest = data(:,11);
gr_sec = data(:,12);
gr_macro = data(:,13);
w_fp = data(:,14); 
w_f = data(:,15);
w_ip = data(:,16);
w_i = data(:,17);
w_p = data(:,18);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: July 10, 2022
%
% The following portion of code
%
%   1) create aggregate variance and components using constant weights
%   2) create aggregate vairance and components using time-varying weights
%   3) compute Herfindahl indices based on the distribution of weights
% 
% References to the following functions: 
%   ag_growth.m (aggregate growth and its variance, based on time-varying 
%       weights)
%   var_growth.m (variance of aggregate growth, based on constant weights
%   var_only_growth.m (variance of aggregate growth, ignoring within 
%       components covariance terms)
%   herfindahl.m (Herfindahl index)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


gr_macrosec = gr_macro+gr_sec;

% Create aggregate growth rates and components and compute standard
% deviations
wgr_variable = w_fp.*gr_variable;
[agr_variable sd_agvariable] = ag_growth(id_fp,period,wgr_variable);
wgr_idio = w_fp.*gr_idio;
[agr_idio sd_agidio] = ag_growth(id_fp,period,wgr_idio);
wgr_macrosec = w_fp.*gr_macrosec;
[agr_macrosec sd_agmacrosec] = ag_growth(id_fp,period,wgr_macrosec);

% Create aggregate variance components: time-varying weights

[VAR_variable var_variable sd_variable] = var_growth(id_fp,period,gr_variable,w_fp);
[VAR_idio var_idio sd_idio] = var_growth(id_fp,period,gr_idio,w_fp);
[VAR_idio_dest var_idio_dest sd_idio_dest] = var_growth(id_fp,period,gr_idio_dest,w_fp);
[VAR_idio_com var_idio_com sd_idio_com] = var_growth(id_fp,period,gr_idio_com,w_fp);
[VAR_macrosec var_macrosec sd_macrosec] = var_growth(id_fp,period,gr_macrosec,w_fp);

% Herfindahl index

herf_fp = herfindahl(id_fp,period,w_fp);
herf_f = herfindahl(id_f,period,w_f);
herf_ip = herfindahl(id_ip,period,w_ip);
herf_i = herfindahl(id_i,period,w_i);
herf_p = herfindahl(id_p,period,w_p);

t = linspace(min(period),max(period),((max(period)-min(period))/25)+1)';


% Results in Table 5 (Panel III)
% Standard Deviations
actual_vol = mean(sd_variable)
firm_vol = mean(sd_idio)
macrosec_vol = mean(sd_macrosec)
% Mean growth rates
actual_mean_gr = mean(agr_variable)
firm_mean_gr = mean(agr_idio)
macrosec_mean_gr = mean(agr_macrosec)

statistics_table = table(actual_vol, firm_vol, macrosec_vol,...
                         actual_mean_gr, firm_mean_gr, macrosec_mean_gr);



% Calculate variances/covariances

gr_sales = gr_variable;
gr_macsec = gr_macro+gr_sec;

% Calculate variances/covariances

[VAR_sales,var_sales,sd_sales] = var_growth(id_fp,period,gr_sales,w_fp);
[VAR_idio,var_idio,sd_idio] = var_growth(id_fp,period,gr_idio,w_fp);
[VAR_macsec,var_macsec,sd_macsec] = var_growth(id_fp,period,gr_macsec,w_fp);

% Calculate standard errors

[var_sales_se,sd_sales_se] = std_err(id_fp,period,gr_sales,w_fp,var_sales);
[var_idio_se,sd_idio_se] = std_err(id_fp,period,gr_idio,w_fp,var_idio);
[var_macsec_se,sd_macsec_se] = std_err(id_fp,period,gr_macsec,w_fp,var_macsec);


% Plotting

export_vol = sd_sales;
export_vol_se_max = sd_sales+1.96*sd_sales_se;
export_vol_se_min = sd_sales-1.96*sd_sales_se;

export_idio_vol = sd_idio;
export_idio_vol_se_max = sd_idio+1.96*sd_idio_se;
export_idio_vol_se_min = sd_idio-1.96*sd_idio_se;

export_macrosec_vol = sd_macsec;
export_macrosec_vol_se_max = sd_macsec+1.96*sd_macsec_se;
export_macrosec_vol_se_min = sd_macsec-1.96*sd_macsec_se;

year = floor(t/100);
quarter = ((t/100-floor(t/100))*100)/25 + 1;

res_table = table(year,quarter,...
                 export_vol,export_vol_se_max,export_vol_se_min,...
                 export_idio_vol,export_idio_vol_se_max,export_idio_vol_se_min,...
                 export_macrosec_vol,export_macrosec_vol_se_max,export_macrosec_vol_se_min);


writetable(res_table,'Output/ts_outputs_capays_destsec_quarterly.csv')
writetable(statistics_table,'Output/aggregate_outputs_capays_destsec_quarterly.csv')



