%% 1. PRELIMINARIES
% =======================================================================
clear all; clear session; close all; clc
warning off

addpath './VAR_Toolbox/VAR'
addpath './VAR_Toolbox/Utils'
addpath './VAR_Toolbox/Auxiliary'
addpath './VAR_Toolbox/ExportFig'
addpath './VAR_Toolbox/Figure'
addpath './VAR_Toolbox/Stats'

country_list = ["Australia", "Austria", "Belgium", "Brazil", "Canada", "Chile", "China", ...
    "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 Kingdom","United States", ];


country_list_unc =["Australia", "Austria", "Belgium", "Brazil", "Canada", "Chile", "China", ...
    "Colombia", "Czech", "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", ...
    "US", "UK"];


%TOTAL TRADE in 2019
GDP_2019_weight =[0.019110889
    0.013367298
    0.024006263
    0.014958951
    0.032464302
    0.005012019
    0.127046765
    0.003566241
    0.009981783
    0.010788194
    0.00122197
    0.005889525
    0.049603466
    0.096123799
    0.004569749
    0.007674108
    0.000563922
    0.028553074
    0.006363796
    0.033863684
    0.049897997
    0.037664166
    0.001169872
    0.002253929
    0.00707264
    0.028078254
    0.003215592
    0.040060953
    0.008370712
    0.017681273
    0.005785454
    0.005407645
    0.002410105
    0.025829652
    0.014393346
    0.024745618
    0.014511805
    0.055078352
    0.161642835];

%equal weight
%GDP_2019_weight = 1/39*ones(39,1);

%load US and control variables
[xlsdata,text] = xlsread('./Data/main_data_2022update.xlsx','United States');

dates_l = 2021+5/12;

nn = size(xlsdata,1);
dates_end = text(nn+1,1);

temp = char(dates_end(1));
year_end = str2double(temp(1:4));
month_end  = str2double(temp(6:end));
date_end = year_end +month_end/12;
index_end = ceil((date_end - dates_l)*12);

xlsdata = xlsdata(1:end-index_end,:);

ip_us = (xlsdata(:,2));
reu_us = (xlsdata(:,1));
ff  = (xlsdata(:,3));
emp  = (xlsdata(:,4));
sp500  = (xlsdata(:,5));
cpi  = (xlsdata(:,6));
vxo = (xlsdata(:,7));

for jj = 1:size(country_list,2)
    
     disp(jj);
    % Load
    [xlsdata,text] = xlsread('./Data/main_data_2022update.xlsx',country_list{jj});
    dates_l = 2021+5/12;

    nn = size(xlsdata,1);
    dates_end = text(nn+1,1);
    
    temp = char(dates_end(1));
    year_end = str2double(temp(1:4));
    month_end  = str2double(temp(6:end));
    date_end = year_end +month_end/12;
    index_end = ceil((date_end - dates_l)*12);
    
    xlsdata = xlsdata(1:end-index_end,:);
    
    [xlsdata_unc,text_unc] = xlsread('./Data/uncertainty_index_wps.xlsx',country_list_unc{jj});
    
    REU = [xlsdata_unc(1:end,1);NaN(0,1)];
    dates_1 = text(2,1);
    
    
    temp = char(dates_1(1));
    year = str2double(temp(1:4));
    month = str2double(temp(6:end));
    
    dates_first = year + month/12;
    
    dates_f = 1960+7/12;
    
    
    index = floor((dates_first - dates_f)*12);
    
    
    
    ip = (xlsdata(:,2));
    
    
    REU_index(:,jj) = [NaN(index,1);REU];
    ip_index(:,jj) = [NaN(index,1);ip];
    
end


T = length(ip);

%% Global REU
Global_reu = func_gdp_weight(REU_index,GDP_2019_weight');
dates = 1960+7/12:1/12:2021+5/12;

plot(dates, Global_reu,'linewidth',2)
grid on
xlim([min(dates),max(dates)])
title('Global Real Economic Uncertainty')
xlabel('Year')
set(gca,'fontsize',12)
dim = [8,6];
set(gcf,'papersize',dim,'paperposition',[0,0,dim]);
figure(gcf)
print(gcf,'-dpdf','./figure/ts_global_tradew_REU_2022') %Change Figure Name Here
close all

%% Figure 4
for jj = 1:size(country_list,2)
    disp(jj);
    
    domestic_reu = REU_index(:,jj);
    domestic_ip = ip_index(:,jj);
    
    %GDP weighted
    weight = GDP_2019_weight([1:(jj-1),(jj+1):size(country_list,2)])./sum(GDP_2019_weight([1:(jj-1),(jj+1):size(country_list,2)]));
    foreign_reu = func_gdp_weight(REU_index(:,[1:(jj-1),(jj+1):size(country_list,2)]),weight');
    foreign_ip = func_gdp_weight(ip_index(:,[1:(jj-1),(jj+1):size(country_list,2)]),weight');
    
    X = 100*[vxo, domestic_reu,foreign_reu,log(domestic_ip),log(foreign_ip)];
    ind1 =isnan(domestic_reu) ==0;
    ind2 =isnan(domestic_ip) ==0;
    ind3 = isnan(vxo) == 0;
    ind = ind1&ind2&ind3;
    
    
    X = X(ind,:);
    vnames = {'VIX','Domestic REU','Foreign REU','Domestic IP','Foreign IP'};
    
    % VAR ESTIMATION
    % =======================================================================
    % Set the case for the VARout (0, 1, or 2)
    det = 2;
    % Set number of nlags
    nlags = 6;
    % Estimate
    [VAR, VARopt] = VARmodel(X,nlags,det);
    % Print at screen and create table
    VARopt.vnames = vnames;
    [beta, tstat, TABLE] = VARprint(VAR,VARopt);
    
    
    % IMPULSE RESPONSE
    % =======================================================================
    % Set options some options for IRF calculation
    VARopt.nsteps = 12*5;
    VARopt.ident = 'oir';
    VARopt.quality = 0;
    VARopt.pick = 1;
    VARopt.pctg = 68;
    % Compute IRF
    [IRF, VAR] = VARir(VAR,VARopt);
    % Compute error bands
    [IRFINF(:,:,:,jj),IRFSUP(:,:,:,jj),IRFMED(:,:,:,jj)] = VARirband(VAR,VARopt);
    % Plot
    % VARirplot(IRFMED,VARopt,IRFINF,IRFSUP,'test',country_list{jj},2);
    
    if  jj <= 20
        subplot(4,5,jj)
        VARirplot_spillover(IRFMED(:,:,:,jj),VARopt,IRFINF(:,:,:,jj),IRFSUP(:,:,:,jj),'test',country_list{jj},4);
    else
        %figure
        subplot(4,5,jj-20)
        VARirplot_spillover(IRFMED(:,:,:,jj),VARopt,IRFINF(:,:,:,jj),IRFSUP(:,:,:,jj),'test',country_list{jj},4);
    end
    
    if jj == 20
        
        dim = [8,6];
        set(gcf,'papersize',dim,'paperposition',[0,0,dim]);
        figure(gcf)
        print(gcf,'-dpdf','./figure/IRF_spillover_REU_tradew_first20') %Change Figure Name Here
        close all
    elseif jj == 39
        
        dim = [8,6];
        set(gcf,'papersize',dim,'paperposition',[0,0,dim]);
        figure(gcf)
        print(gcf,'-dpdf','./figure/IRF_spillover_REU_tradew_next19') %Change Figure Name Here
        close all
        
    end
    
    % FORECAST ERROR VARIANCE DECOMPOSITION
    % % =======================================================================
    % % Compute FEVD
    VARopt.pick = 4;
    [FEVD, VAR] = VARfevd(VAR,VARopt);
    
    FEVD_domesticREU_avg(jj) = mean(FEVD(1:24,4,2));
    FEVD_foreignREU_avg(jj) = mean(FEVD(1:24,4,3));
    FEVD_domesticREU_max(jj) = max(FEVD(1:24,4,2));
    FEVD_foreignREU_max(jj) = max(FEVD(1:24,4,3));
    
end
