 
% wavelet cohesion  
% distances  
 
% 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'); 
 
ta = 1;        te = 148;      
dd41 = (log(dat1(ta+4:te,2:end)) - log(dat1(ta:te-4,2:end))) * 100; 
dd42 = (log(dat2(ta+4:te,2:end)) - log(dat2(ta:te-4,2:end))) * 100; 
dd43 = (log(dat3(ta+4:te,2:end)) - log(dat3(ta:te-4,2:end))) * 100; 

indd = 1;       % selects country, 1 is DE 
dd4a = [ dd41(1:end,indd) dd42(1:end,indd) dd43(1:end,indd) ]; 
ttime = dat1(ta+4:te,1); 
ntt = size(ttime,1); 
 
mnames = char('LNF', 'LHH', 'RPP'); 

% 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 = 0;     % 100  
sur_type = 'ARMAEcon'; 
p = 4;    q = 1; 
  
% ######################################################################## 
% cohesion 
ncc = size(dd4a,2);    sncc = sum(1:ncc-1); 
nwave = 90; 
COHES = NaN(nwave,ntt,sncc); 
cind = 1; 
nSur = 0; 
for cc1 = 1 : ncc-1                
    for cc2 = cc1+1 : ncc         
        ts1 = dd4a(:,cc1); 
        ts2 = dd4a(:,cc2); 
        
        % wavelet coherency  
        [ WCO, d1, d2, period, d3, coi, pvCO, WCoh ] = ...
            AWCOh(ts1,ts2,dt,dj,lowPeriod,upPeriod,[],[],[],[],... 
            wT_type,wT_size,wS_type,wS_size,nSur,sur_type,p,q); 
        COHES(:,:,cind) = WCoh;
        cind = cind + 1; 
    end 
end 

% unweighted cohesion 
CCO1 = mean(COHES,3); 
 
% ######################################################################## 
% bootstrapping cohesion 
% simulating series  
nSurB = 1000; 
surXX = NaN(ntt,nSurB,ncc);  
for vvi = 1 : ncc 
    xi = dd4a(:,vvi); 
    surxi = SurrogateARMAEcon15(xi,p,q,nSurB); 
    surXX(:,:,vvi) = surxi; 
end 
 
% calculating cohesions 
rowsWCO = 90;      
colsWCO = ntt;     
COHESB = NaN(rowsWCO,colsWCO,nSurB); 
for bbi = 1 : nSurB  
    COHESi = NaN(rowsWCO,colsWCO,sncc); 
    cind = 1; 
    for cc1 = 1 : ncc-1               
        for cc2 = cc1+1 : ncc  
            ts1 = surXX(:,bbi,cc1); 
            ts2 = surXX(:,bbi,cc2); 
         
            % wavelet coherency  
            [ WCO, d1, d2, period, d3, coi, pvCO, WCoh ] = ...
                AWCOh(ts1,ts2,dt,dj,lowPeriod,upPeriod,[],[],[],[],... 
                wT_type,wT_size,wS_type,wS_size,nSur,sur_type,p,q); 
            COHESi(:,:,cind) = WCoh;
            cind = cind + 1; 
        end    
    end  
    % equal weights 
    CCO1i(:,:,bbi) = mean(COHESi,3);  
end   
  
% calculating confidence regions 
nRows = size(CCO1,1); 
nCols = size(CCO1,2);  
pvCOB = zeros(nRows,nCols); 
for iSurB = 1 : nSurB  
    CCO1ib = CCO1i(:,:,iSurB); 
    for jj = 1 : nRows 
        for kk = 1 : nCols 
            if CCO1ib(jj,kk) >= CCO1(jj,kk) 
                pvCOB(jj,kk) = pvCOB(jj,kk) + 1; 
            end  
        end  
    end  
end  
  
% save('Coh_DE_FS3.mat','CCO1')
% save('pvCO_DE_FS3.mat','pvCOB') 

% ########################################################################
% graph 
xlim = [ min(ttime) max(ttime) ]; 
xticks = [ 1990 2000 2010 ]; 
ylim = log2( [ min(period),max(period) ] );
yticks = [ 3 6 10 16 24 ]; 

figure(2) 
imagesc(ttime,log2(period),CCO1); 
title('Cohesion DE (FS3)','FontSize',10);
ylabel('Period (years)','FontName','times','FontSize',10);
set(gca,'XGrid','on','XTick',xticks,'FontSize',8);
set(gca,'YLim',log2([min(period),max(period)]),'YDir','reverse', ...
        'YTick',log2(yticks), ...
        'YTickLabel',yticks,'YGrid','on','FontSize',8);     
caxis manual 
caxis([ -1.0 1.0 ]); 
colormap(jet) 
% colorbar('location','eastoutside')
hold on
        
% COI
plot(ttime,log2(coi),'-r','LineWidth',2.5); 
 
% levels of significance
if nSurB > 0
    contour(ttime,log2(period),(pvCOB),[5,5],'k-','LineWidth',1.5);
end 
hold off 
 




