% Calculating wavelet coherencies and time differences relative to GDPR    
%  using CWT-Toolbox written by Aguiar-Conraria, L. and Soares, M.J. (2012) 

% clear all variables and screen  
clear all; 
clc;

% update path   
addpath('.\ASToolbox\Functions\Auxiliary');
addpath('.\ASToolbox\Functions\WaveletTransforms');
 
% seed 
defaultStream = RandStream.getGlobalStream;
savedState = defaultStream.State;  
 
% data
dat1 = xlsread('Data_LNF_EA4.xlsx');  
dat2 = xlsread('Data_LHH_EA4.xlsx'); 
dat3 = xlsread('Data_RPP_EA4.xlsx'); 
dat4 = xlsread('Data_YER_EA6.xlsx'); 
 
[ nT, nvar ] = size(dat1);  
 
ta = 1;        te = 148; 
dd4a(:,:,1) = (log(dat1(ta+4:te,2:end)) - log(dat1(ta:te-4,2:end))) * 100; 
dd4a(:,:,2) = (log(dat2(ta+4:te,2:end)) - log(dat2(ta:te-4,2:end))) * 100; 
dd4a(:,:,3) = (log(dat3(ta+4:te,2:end)) - log(dat3(ta:te-4,2:end))) * 100; 
dd4a(:,:,4) = (log(dat4(ta+4:te,2:end)) - log(dat4(ta:te-4,2:end))) * 100; 
ttime = dat1(ta+4:te,1); 
ntt = size(ttime,1); 
 
cnames = char('DE', 'ES', 'FR', 'IT'); 
fnames = char('LNF', 'LHH', 'RPP', 'YR'); 
   
% specification for wavelet analysis 
dt = 1 / 4;    dj = 1 / 25; 
lowPeriod = 2.0;      upPeriod = 24.0; 
pad = 0; 
mother = 'Morlet'; 
be = 6.0; 
ga = []; 

wT_type = 'Ham'; 
wT_size = 2; 
wS_type = 'Ham'; 
wS_size = 2; 
 
nSur = 100; 
sur_type = 'ARMAEcon'; 
p = 4; 
q = 1; 

figure(4) 
% loop over countries    
ncc = 4;   nfs = 3; 
cind = 1; 
timeDiff1 = NaN(ncc,ntt);  
timeDiff2 = NaN(ncc,ntt);  
lowPF1 = 6;     upPF1 = 10; 
lowPF2 = 10;    upPF2 = 16; 
for ff = 1 : nfs 
    for ii = 1 : ncc 
        d4lfs = dd4a(:,ii,ff);    
        d4lyr = dd4a(:,ii,4);  
        
        % wavelet coherency  
        [ WCO, d1, d2, period, d3, coi, pvCO ] = ...
            AWCOh(d4lfs,d4lyr,dt,dj,lowPeriod,upPeriod,[],[],[],[],... 
            wT_type,wT_size,wS_type,wS_size,nSur,sur_type,p,q);
   
        % time differences 
        [ d1, d2, d3, dphaseXY1, timeLag1 ] = ... 
              AWCOOutput(WCO,(WCO),period,lowPF1,upPF1);
        timeDiff1(cind,:) = timeLag1; 
           
        [ d1, d2, d3, dphaseXY2, timeLag2 ] = ... 
              AWCOOutput(WCO,(WCO),period,lowPF2,upPF2);
        timeDiff2(cind,:) = timeLag2; 

        
        xlim = [ min(ttime) max(ttime) ]; 
        xticks = [ 1990 2000 2010 ]; 
        ylim = log2( [ min(period),max(period) ] );
        yticks = [ 3 6 10 16 ]; 
        pictEnh = 7.5;   
        perc5 = 5/100;        
 
        subplot(nfs,ncc,cind); 
        imagesc(ttime,log2(period),abs(WCO).^pictEnh); 
        colormap(jet) 
        if ff == 1 
            title(cnames(ii,:),'FontSize',8); 
            if ii == 1 
                ylabel('LNF','FontName','times','FontWeight','bold','FontSize',8);
            end  
        end 
        if ff == 2 && ii == 1 
            ylabel('LHH','FontName','times','FontWeight','bold','FontSize',8);  
        end 
        if ff == 3 && ii == 1   
            ylabel('RPP','FontName','times','FontWeight','bold','FontSize',8); 
        end   

        set(gca,'XGrid','on','XTick',xticks,'FontSize',6);
        set(gca,'YLim',log2([min(period),max(period)]),'YDir','reverse', ...
                'YTick',log2(yticks), ...
                'YTickLabel',yticks,'YGrid','on','FontSize',6);     
        hold on 
       
        % COI
        plot(ttime,log2(coi),'-r','LineWidth',2.5);
       
        % levels of significance
        if nSur>0
             contour(ttime,log2(period),(pvCO),[perc5,perc5],...
                    'k-','LineWidth',1.5);
        end 
        hold off 
        cind = cind + 1; 
    end 
end 
 
 
% ######################################################################## 
% time differences 
%       relative to GDPR      
   
timeDiff1B = NaN(12,ntt); 
timeDiff1B(:,23:122) = timeDiff1(:,23:122); 
timeDiff2B = NaN(12,ntt); 
timeDiff2B(:,36:109) = timeDiff2(:,36:109); 
  
% 6_10 
figure(5) 
plot(ttime,timeDiff1(1,:),'--r',ttime,timeDiff1(2,:),'--b', ... 
     ttime,timeDiff1(3,:),'--g',ttime,timeDiff1(4,:),'--k', ... 
     'LineWidth',2.0);
hold on 
plot(ttime,timeDiff1B(1,:),'-r',ttime,timeDiff1B(2,:),'-b', ... 
     ttime,timeDiff1B(3,:),'-g',ttime,timeDiff1B(4,:),'-k', ... 
     'LineWidth',2.0);
title('LNF: Time difference (relative to GDPR)',... 
                          'FontName','arial','FontSize',10);
ylabel('Time difference','FontName','times','FontSize',8);
set(gca,'XLim',xlim,'FontSize',7,'XGrid','off','XTick',xticks); 
ord = (lowPF1 + upPF1) / 4; 
set(gca,'YLim',[-ord ord],'FontSize',7,'YGrid','on');
set(gca,'YTick',-ord:ord/2:ord,'YTickLabel',{-ord,-ord/2,0,ord/2,ord}) 
legend('DE', 'ES', 'FR', 'IT', ... 
       'Location','SouthOutside','Orientation','horizontal')
 

figure(7) 
plot(ttime,timeDiff1(5,:),'--r',ttime,timeDiff1(6,:),'--b', ... 
     ttime,timeDiff1(7,:),'--g',ttime,timeDiff1(8,:),'--k', ... 
     'LineWidth',2.0);
hold on 
plot(ttime,timeDiff1B(5,:),'-r',ttime,timeDiff1B(6,:),'-b', ... 
     ttime,timeDiff1B(7,:),'-g',ttime,timeDiff1B(8,:),'-k', ...
     'LineWidth',2.0);
title('LHH: Time difference (relative to GDPR)',... 
                          'FontName','arial','FontSize',10);
ylabel('Time difference','FontName','times','FontSize',8);
set(gca,'XLim',xlim,'FontSize',7,'XGrid','off','XTick',xticks); 
ord = (lowPF1 + upPF1) / 4; 
set(gca,'YLim',[-ord ord],'FontSize',7,'YGrid','on');
set(gca,'YTick',-ord:ord/2:ord,'YTickLabel',{-ord,-ord/2,0,ord/2,ord}) 
legend('DE', 'ES', 'FR', 'IT', ... 
       'Location','SouthOutside','Orientation','horizontal')
 
  
% ######################################################################## 
% 10_16 
figure(6) 
plot(ttime,timeDiff2(1,:),'--r',ttime,timeDiff2(2,:),'--b', ... 
     ttime,timeDiff2(3,:),'--g',ttime,timeDiff2(4,:),'--k', ... 
     'LineWidth',2.0);
hold on 
plot(ttime,timeDiff2B(1,:),'-r',ttime,timeDiff2B(2,:),'-b', ... 
     ttime,timeDiff2B(3,:),'-g',ttime,timeDiff2B(4,:),'-k', ... 
     'LineWidth',2.0);
title('LNF: Time difference (relative to GDPR)',... 
                          'FontName','arial','FontSize',10);
ylabel('Time difference','FontName','times','FontSize',8);
set(gca,'XLim',xlim,'FontSize',7,'XGrid','off','XTick',xticks); 
ord = (lowPF2 + upPF2) / 4; 
set(gca,'YLim',[-ord ord],'FontSize',7,'YGrid','on');
set(gca,'YTick',-ord:ord/2:ord,'YTickLabel',{-ord,-ord/2,0,ord/2,ord}) 
legend('DE', 'ES', 'FR', 'IT', ... 
       'Location','SouthOutside','Orientation','horizontal')
 

figure(8) 
plot(ttime,timeDiff2(5,:),'--r',ttime,timeDiff2(6,:),'--b', ... 
     ttime,timeDiff2(7,:),'--g',ttime,timeDiff2(8,:),'--k', ... 
     'LineWidth',2.0);
hold on 
plot(ttime,timeDiff2B(5,:),'-r',ttime,timeDiff2B(6,:),'-b', ... 
     ttime,timeDiff2B(7,:),'-g',ttime,timeDiff2B(8,:),'-k', ...
     'LineWidth',2.0);
title('LHH: Time difference (relative to GDPR)',... 
                          'FontName','arial','FontSize',10);
ylabel('Time difference','FontName','times','FontSize',8);
set(gca,'XLim',xlim,'FontSize',7,'XGrid','off','XTick',xticks); 
ord = (lowPF2 + upPF2) / 4; 
set(gca,'YLim',[-ord ord],'FontSize',7,'YGrid','on');
set(gca,'YTick',-ord:ord/2:ord,'YTickLabel',{-ord,-ord/2,0,ord/2,ord}) 
legend('DE', 'ES', 'FR', 'IT', ... 
       'Location','SouthOutside','Orientation','horizontal')
 
 
figure(9) 
plot(ttime,timeDiff2(9,:),'--r',ttime,timeDiff2(10,:),'--b', ... 
     ttime,timeDiff2(11,:),'--g',ttime,timeDiff2(12,:),'--k', ... 
     'LineWidth',2.0);
hold on 
plot(ttime,timeDiff2B(9,:),'-r',ttime,timeDiff2B(10,:),'-b', ... 
     ttime,timeDiff2B(11,:),'-g',ttime,timeDiff2B(12,:),'-k', ... 
     'LineWidth',2.0);
title('RPP: Time difference (relative to GDPR)',... 
                          'FontName','arial','FontSize',10);
ylabel('Time difference','FontName','times','FontSize',8);
set(gca,'XLim',xlim,'FontSize',7,'XGrid','off','XTick',xticks); 
ord = (lowPF2 + upPF2) / 4; 
set(gca,'YLim',[-ord ord],'FontSize',7,'YGrid','on');
set(gca,'YTick',-ord:ord/2:ord,'YTickLabel',{-ord,-ord/2,0,ord/2,ord}) 
legend('DE', 'ES', 'FR', 'IT', ... 
       'Location','SouthOutside','Orientation','horizontal')
 
 
