%% 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.021914715	0.0164005	0.02965338	0.01764604	0.038664162	0.012567189	0.007332389	0.06156293	0.116631872	0.005643817	0.000713125	0.031903312	0.007598036	0.041332126	0.061980673	0.044132257	0.008305101	0.033781825	0.048166481	0.003886935	0.009913463	0.021590137	0.007153702	0.032029123	0.014444005	0.017147829	0.016540975	0.016278826	0.062381765	0.192703306];


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];



%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));
    if jj < 39
        ip_manu =  (xlsdata(:,3));
        ip_service =  (xlsdata(:,4));
    else %US
        ip_manu =  (xlsdata(:,8));
        ip_service =  (xlsdata(:,9));
    end
    
    REU_index(:,jj) = [NaN(index,1);REU];
    
    
    ip_index(:,jj) = [NaN(index,1);ip];
    ip_manu_index(:,jj) = [NaN(index,1);ip_manu];
    ip_service_index(:,jj) = [NaN(index,1);ip_service];
    
end

T = length(ip);
%%

for jj = 1:size(country_list,2)
    disp(jj);
    
    for kkk = 1:2
        domestic_reu = REU_index(:,jj);
        
        if kkk == 1
            domestic_ip = ip_manu_index(:,jj);
        else
            domestic_ip = ip_service_index(:,jj);
        end
        
        %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');
        
        
        %         foreign_reu = nanmean(REU_index(:,[1:(jj-1),(jj+1):size(country_list,2)])')';
        %         foreign_ip = nanmean(ip_index(:,[1:(jj-1),(jj+1):size(country_list,2)])')';
        %
        %
        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 = 2;
        VARopt.pctg = 68;
        % Compute IRF
        [IRF, VAR] = VARir(VAR,VARopt);
        % Compute error bands
        
        if kkk == 1
            [IRFINF_manu(:,:,:,jj),IRFSUP_manu(:,:,:,jj),IRFMED_manu(:,:,:,jj)] = VARirband(VAR,VARopt);
        else
            
            [IRFINF_service(:,:,:,jj),IRFSUP_service(:,:,:,jj),IRFMED_service(:,:,:,jj)] = VARirband(VAR,VARopt);
        end
    end
end


[nsteps, nvars, nshocks] = size(IRF);
x_axis = zeros(1,nsteps);
report_out = [];
for jj = 1:39
    
    
    subplot(5,8,jj)
    % VARirplot_spillover(IRFMED(:,:,:,jj),VARopt,IRFINF(:,:,:,jj),IRFSUP(:,:,:,jj),'test',country_list{jj},4);
    if jj == 11
        plot(IRFMED_manu(:,4,3,jj),'LineStyle','-','Color',[0.01 0.09 0.44],'LineWidth',2);
    else
        plot(IRFMED_manu(:,4,3,jj),'LineStyle','-','Color',[0.01 0.09 0.44],'LineWidth',2);
    end
    hold on
    plot(IRFMED_service(:,4,3,jj),'LineStyle','-','Color','r','LineWidth',2);
    hold on
    ciplot_h(IRFINF_manu(:,4,3,jj),IRFSUP_manu(:,4,3,jj),1:60,[0.01 0.09 0.44],0.1);
    hold on
    ciplot_h(IRFINF_service(:,4,3,jj),IRFSUP_service(:,4,3,jj),1:60,'r',0.1);
    
    hold on
    plot(x_axis,'k','LineWidth',0.5)
    title(country_list{jj})
    
    if jj == 1
        leg = legend('Manu.','Retail Sales');
        set(leg','orientation','horizontal')
    end
    
    report_temp = [IRFMED_manu(:,4,3,jj),IRFMED_service(:,4,3,jj)];
    report_out = [report_out, report_temp];
    
    
end
dim = [12,6];
set(gcf,'papersize',dim,'paperposition',[0,0,dim]);
figure(gcf)
print(gcf,'-dpdf','./figure/IRF_spillover_manu_service_2022') %Change Figure Name Here
%close all
