% Code for section 3.3. of "The Global Transmission of Real Economic
% Uncertainty", JMCB (forthcoming)

% Juan-Miguel Londono

clear all;

% Add librarias
addpath(genpath('./Data'));
%addpath(genpath('U:\fm\matlab_libs'));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load list of countries.
% First column are the names of the files that contain trade matrices
% Second column are the names of the countries used by the IMF in the
% matrices
filename=char('./Data/country_names.xlsx');
[NUM,c_names,RAW]=xlsread(filename,'Sheet1','A2:B41');
file_names=c_names(:,1);        
c_names=c_names(:,2);  


% Control variables
% 1. VIX and VXO from CBOE
load('mvix.mat');
load('vxo.mat');

% Load transmission variables: financial and trade integration/exposure,
% and the control variables: size and distance between capitals. 
% See file "trade_indexes.m"
load('info.mat');
% Create unconditional estimates of all transmission indexes
% For all matrices, the elements are organized as follows: 
% [ a_(1,1) a_(1,2) .....a_(1,n)
% [ a_(2,1) a_2,3) .........
% [ ....
% [ a_(n,1)...............a_(n,n)]
%
% element a_(i,j) corresponds to any measure of trade or assets between i
% and j, divided by country's i GDP


imp_gdp=[];
exp_gdp=[];
trade_gdp=[];
ca_gdp=[];
relsize=[];
asset_gdp=[];
liab_gdp=[];
flow_gdp=[];
nd_gdp=[];
for j=1:length(c_names)-1
    imp_gdp=[imp_gdp;nanmean(info(j).imp_gdp(1:end-1,2:end-1))];
    exp_gdp=[exp_gdp;nanmean(info(j).exp_gdp(1:end-1,2:end-1))];
    trade_gdp=[trade_gdp;nanmean(info(j).trade_gdp(1:end-1,2:end-1))];
    ca_gdp=[ca_gdp;nanmean(info(j).ca_gdp(1:end-1,2:end-1))];
    relsize=[relsize;nanmean(info(j).relsize(1:end-1,2:end-1))];
    asset_gdp=[asset_gdp;nanmean(info(j).asset_gdp(1:end-1,2:end-1))];
    liab_gdp=[liab_gdp;nanmean(info(j).liab_gdp(1:end-1,2:end-1))];
    flow_gdp=[flow_gdp;nanmean(info(j).flow_gdp(1:end-1,2:end-1))];
    nd_gdp=[nd_gdp;nanmean(info(j).nd_gdp(1:end-1,2:end-1))];
end
distance=info(1).distance;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Exploration of the relation between IRFs and the transmission variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LOad average IRFs of effect of shocks to foreign REU on domestic IP
% Averaged over 24 months or 60 monhts
% Order of countries
%["Australia","Austria","Belgium","Canada","Chile","Colombia","Czech Republic","Denmark",
%"Estonia","Finland","France","Germany","Greece","Hungary","Iceland","Ireland","Israel",
%"Italy","Japan","Korea","Latvia","Lithuania","Luxembourg","Mexico","Netherlands","New Zealand",
%"Norway","Poland","Portugal","Slovak Republic","Slovenia","Spain","Sweden","Switzerland","Turkey",
%"United States","United Kingdom","Brazil","China"]]

% Calculate average IRFs
%load('U:\IFS\REU\JMCB\output\out_files_all_ji_newChinaChileChe.mat ')
%newout=out_IRF_fREUip;
load('newout.mat')

irf24=nan(length(c_names)-1,length(c_names)-1);
for i=1:length(c_names)-1
    for j=1:length(c_names)-1
        newnewout=newout(i,j,:);
        irf24(i,j)=mean(newnewout(2:24));
    end
end

% Match to make sure the integration variables are in the same order as the
% IRF countries
sel=[1	2	3	5	6	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32	33	34	35	36	37	39	38	4	7];
newirf24=[];
%newirf60=[];

newtrade=[];
newca=[];
newrelsize=[];
newflow=[];
newdistance=[];
for j=1:length(sel)
    for jj=1:length(sel)
        if j==jj
            continue
        else   
        newtrade=[newtrade;trade_gdp(sel(j),sel(jj))];
        newca=[newca;ca_gdp(sel(j),sel(jj))];
        newrelsize=[newrelsize;relsize(sel(j),sel(jj))];
        newflow=[newflow;flow_gdp(sel(j),sel(jj))];
        newdistance=[newdistance;1./distance(sel(j),sel(jj))];
        newirf24=[newirf24;[irf24(j,jj) j jj]];
        %newirf60=[newirf60;[irf60(j,jj) j jj]];
        end
    end
end


aa=[newirf24(:,1) normalize(newtrade) normalize(newca) normalize(newrelsize) normalize(newflow) normalize(newdistance) newirf24(:,2:3)];
aa=aa(~any(isnan(aa),2),:);


n=1;
yvar(1).var=aa(:,2);
n=n+1;
yvar(n).var=aa(:,3);
n=n+1;
yvar(n).var=aa(:,4);
n=n+1;
yvar(n).var=aa(:,5);
n=n+1;
yvar(n).var=aa(:,6);


xvar(1).var=aa(:,1);
%xvar(2).var=aa(:,2);


% Scatter plots and tables for unconditional setting


% Tables
% Univariate
table1=[];
%table1_60=[];
% Univariate with very high very low
table2=[];
%table2_60=[];
clear prlevel
prlevel(1).p=1;
prlevel(2).p=2;
prlevel(3).p=5;
prlevel(4).p=10;


for j=[2 4 5 6]
newout=nwest(aa(:,1),[ones(length(aa),1) aa(:,j)],0);
table1=[table1 [newout.beta(2);newout.tstat(2)]];
newout=nwest(aa(:,2),[ones(length(aa),1) aa(:,j)],0);
%table1_60=[table1_60 [newout.beta(2);newout.tstat(2)]];
% Dummies for extreme values
newtable2=[];
%newtable2_60=[];
for jj=1:size(prlevel,2)
    dlow=zeros(length(aa),1);
    dlow(aa(:,j)<=prctile(aa(:,j),prlevel(jj).p),1)=1;
    dhigh=zeros(length(aa),1);
    dhigh(aa(:,j)>=prctile(aa(:,j),100-prlevel(jj).p),1)=1;
    newout=nwest(aa(:,1),[ones(length(aa),1) aa(:,j) dlow dhigh],0);
    newtable2=[newtable2;[newout.beta(2);newout.tstat(2)]];
    newout=nwest(aa(:,2),[ones(length(aa),1) aa(:,j) dlow dhigh],0);
    %newtable2_60=[newtable2_60;[newout.beta(2);newout.tstat(2)]];
end
table2=[table2 newtable2];
%table2_60=[table2_60 newtable2_60];
end

% Multivariate: 
% [trade size flow proximity]
out_multi=nwest(aa(:,1),[ones(length(aa),1) aa(:,[2 4 5 6])],0);
table3=[[out_multi.beta(1);out_multi.tstat(1)] [out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];


%out_multi=nwest(aa(:,2),[ones(length(aa),1) aa(:,[3 5 6 7])],0);
%table3_60=[[out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];


table4=[];
%table4_60=[];
for jj=1:size(prlevel,2)
    newtable4=[];
 %   newtable4_60=[];
    dlow=[];
    dhigh=[];
    %for j=[3 5 6 7]
    for j=[2 5]
    newdlow=zeros(length(aa),1);
    newdlow(aa(:,j)<=prctile(aa(:,j),prlevel(jj).p),1)=1;
    newdhigh=zeros(length(aa),1);
    newdhigh(aa(:,j)>=prctile(aa(:,j),100-prlevel(jj).p),1)=1;
    dlow=[dlow newdlow];
    dhigh=[dhigh newdhigh];
    end
    out_multi=nwest(aa(:,1),[ones(length(aa),1) aa(:,[2 4 5 6]) dlow dhigh],0);
    table4=[table4;[out_multi.beta(1);out_multi.tstat(1)] [out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];
    %out_multi=nwest(aa(:,2),[ones(length(aa),1) aa(:,[3 5 6 7]) dlow dhigh],0);
    %table4_60=[table4_60;[out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];
end

% remove offshore financial centers: Luxembourg and Switzerland
newaa=aa;
newaa(newaa(:,end-1)==25,:)=[];
newaa(newaa(:,end)==25,:)=[];
newaa(newaa(:,end-1)==36,:)=[];
newaa(newaa(:,end)==36,:)=[];

out_multi=nwest(newaa(:,1),[ones(length(newaa),1) newaa(:,[2 4 5 6])],0);
table4_ofc=[[out_multi.beta(1);out_multi.tstat(1)] [out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];

% Remove same countries as in the conditional version 
% Remove countries for which we do not have enough IP data:

% Benchmark starting in 1990
% China
% Estonia
% Iceland
% Latvia
% Lithuania
% Slovenia


newaa=aa;
%newaa(newaa(:,end-1)==25,:)=[];
%newaa(newaa(:,end)==25,:)=[];
newaa(newaa(:,end-1)==7,:)=[];
newaa(newaa(:,end)==7,:)=[];
newaa(newaa(:,end-1)==11,:)=[];
newaa(newaa(:,end)==11,:)=[];
newaa(newaa(:,end-1)==17,:)=[];
newaa(newaa(:,end)==17,:)=[];
newaa(newaa(:,end-1)==23,:)=[];
newaa(newaa(:,end)==23,:)=[];
newaa(newaa(:,end-1)==24,:)=[];
newaa(newaa(:,end)==24,:)=[];
newaa(newaa(:,end-1)==33,:)=[];
newaa(newaa(:,end)==33,:)=[];



out_multi=nwest(newaa(:,1),[ones(length(newaa),1) newaa(:,[2 4 5 6])],0);
table4_remove=[[out_multi.beta(1);out_multi.tstat(1)] [out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];



table4_removeoutliers=[];
%table4_60=[];
for jj=1:size(prlevel,2)
    newtable4=[];
 %   newtable4_60=[];
    dlow=[];
    dhigh=[];
    %for j=[3 5 6 7]
    for j=[2 5]
    newdlow=zeros(length(newaa),1);
    newdlow(newaa(:,j)<=prctile(newaa(:,j),prlevel(jj).p),1)=1;
    newdhigh=zeros(length(newaa),1);
    newdhigh(newaa(:,j)>=prctile(newaa(:,j),100-prlevel(jj).p),1)=1;
    dlow=[dlow newdlow];
    dhigh=[dhigh newdhigh];
    end
    out_multi=nwest(newaa(:,1),[ones(length(newaa),1) newaa(:,[2 4 5 6]) dlow dhigh],0);
    table4_removeoutliers=[table4_removeoutliers;[out_multi.beta(1);out_multi.tstat(1)] [out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];
    %out_multi=nwest(aa(:,2),[ones(length(aa),1) aa(:,[3 5 6 7]) dlow dhigh],0);
    %table4_60=[table4_60;[out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];
end



% remove offshore financial centers: Luxembourg and Switzerland
newaa=aa;
newaa(newaa(:,end-1)==7,:)=[];
newaa(newaa(:,end)==7,:)=[];
newaa(newaa(:,end-1)==11,:)=[];
newaa(newaa(:,end)==11,:)=[];
newaa(newaa(:,end-1)==17,:)=[];
newaa(newaa(:,end)==17,:)=[];
newaa(newaa(:,end-1)==23,:)=[];
newaa(newaa(:,end)==23,:)=[];
newaa(newaa(:,end-1)==24,:)=[];
newaa(newaa(:,end)==24,:)=[];
newaa(newaa(:,end-1)==33,:)=[];
newaa(newaa(:,end)==33,:)=[];
newaa(newaa(:,end-1)==25,:)=[];
newaa(newaa(:,end)==25,:)=[];
newaa(newaa(:,end-1)==36,:)=[];
newaa(newaa(:,end)==36,:)=[];

out_multi=nwest(newaa(:,1),[ones(length(newaa),1) newaa(:,[2 4 5 6])],0);
table4_removeofc=[[out_multi.beta(1);out_multi.tstat(1)] [out_multi.beta(2);out_multi.tstat(2)] [out_multi.beta(3);out_multi.tstat(3)] [out_multi.beta(4);out_multi.tstat(4)] [out_multi.beta(5);out_multi.tstat(5)]];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATE SUMMARY OF IRFs and Variance decomposition similar to DY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%load('U:\IFS\REU\JMCB\output\out_files_all_ji_newChinaChileChe.mat ')
load('out_IRF_fREUdREU')
ncountries=length(c_names)-1;

% Average IRF foreign REU to domestic REU


new=100*out_IRF_fREUdREU;


table=nan(ncountries,ncountries);
for j=1:ncountries
    for jj=1:ncountries
        table(j,jj)=mean(new(j,jj,2:24));
    end
end
% Let's do the same as in DY and average the "to" and "from"
aa=diag(table);
%table=[[table (sum(table,2)-aa)./ncountries];[(sum(table,1)-aa')./ncountries nan]];
%table=[table;[table(end,1:end-1)-table(1:end-1,end)' nan]];
% [OWN from_other to_others]
table=[aa (sum(table,2)-aa)./ncountries (sum(table,1)-aa')'./ncountries];
tablereu=table;

% Average IRF foreign REU to IP
load('out_IRF_fREUip.mat');
new=out_IRF_fREUip;

ncountries=length(c_names)-1;
table=nan(ncountries,ncountries);
for j=1:ncountries
    for jj=1:ncountries
        table(j,jj)=mean(new(j,jj,2:24));
    end
end
% Let's do the same as in DY and average the "to" and "from"
aa=diag(table);
%table=[[table (sum(table,2)-aa)./ncountries];[(sum(table,1)-aa')./ncountries nan]];
%table=[table;[table(end,1:end-1)-table(1:end-1,end)' nan]];
% [OWN from_other to_others]
table=[aa (sum(table,2)-aa)./ncountries (sum(table,1)-aa')'./ncountries];
table=[tablereu table];
table=[table;mean(table)];

table_irf_24=table;




% NOw let us do the variance decomposition for the relevance

load('out_VARd_fREUdREU.mat')
new=100*out_VARd_fREUdREU;

table=new;
% Let's do the same as in DY and average the "to" and "from"
aa=diag(table);
%table=[[table (sum(table,2)-aa)./ncountries];[(sum(table,1)-aa')./ncountries nan]];
%table=[table;[table(end,1:end-1)-table(1:end-1,end)' nan]];
% [OWN from_other to_others]
table=[aa (sum(table,2)-aa)./ncountries (sum(table,1)-aa')'./ncountries];
tablereu=table;

% Average IRF foreign REU to IP

load('out_VARd_fREUip.mat')
new=100*out_VARd_fREUip;
table=new;
% Let's do the same as in DY and average the "to" and "from"
aa=diag(table);
%table=[[table (sum(table,2)-aa)./ncountries];[(sum(table,1)-aa')./ncountries nan]];
%table=[table;[table(end,1:end-1)-table(1:end-1,end)' nan]];
% [OWN from_other to_others]
table=[aa (sum(table,2)-aa)./ncountries (sum(table,1)-aa')'./ncountries];
table=[tablereu table];
table=[table;mean(table)];

table_vard_24=table;