% Calculation of wavelet power spectra    
%    using CWT-Toolbox written by Aguiar-Conraria, L. and Soares, M.J. (2012),
%    "The Continuous Wavelet Transform: Moving Beyond Uni-and Bivariate
%    Analysis"
 
% 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_EA4.xlsx'); 

[ nT, nvar ] = size(dat1);   
  
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; 
dd44 = (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); 
 
% indices     
dd4a(:,:,1) = dd41(1:end,:); 
dd4a(:,:,2) = dd42(1:end,:); 
dd4a(:,:,3) = dd43(1:end,:); 
dd4a(:,:,4) = dd45(1:end,:); 
 
cnames = char('DE', 'ES', 'FR', 'IT'); 
fnames = char('LNF', 'LHH', 'RPP', 'GDPR'); 

% 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 = 1000;           
sur_type = 'ARMAEcon'; 
p = 0;    q = 0; 
  
% wavelet power spectra  
figure(1) 
% loop over financial series     
ncc = 4;    nfs = 4; 
cind = 1; 
for ff = 1 : nfs 
    for ii = 1 : ncc    
        d4llfsr = dd4a(:,ii,ff);  
        
        % wavelet power spectrum 
        [ wave, period, scale, coi, power, pvPower ] = ...
                AWTh(d4llfsr,dt,dj,lowPeriod,upPeriod,[],[],[],[],nSur,sur_type,p,q);
 
        % local maxima of WPS   
        maxPower = MatrixMax(power,5,.14);
   
        xlim = [ min(ttime) max(ttime) ]; 
        xticks = [ 1990 2000 2010 ]; 
        ylim = log2( [ min(period),max(period) ] );
        yticks = [ 3 6 10 16 24 ]; 
  
        pictEnh = 1/5;      
        perc5 = 5/100;      
  
        subplot(nfs,ncc,cind) 
        imagesc(ttime,log2(period),(power).^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   
        if ff == 4 && ii == 1 
            ylabel('GDPR','FontName','times','FontWeight','bold','FontSize',8);
        end   
     
        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);     
        hold on
          
        % ridges
        contour(ttime,log2(period),maxPower,[1,1],'w-','LineWidth',1);
        
        % COI
        plot(ttime,log2(coi),'-r','LineWidth',2.5);
        
        % levels of significance
        if nSur>0
             contour(ttime,log2(period),(pvPower),[perc5,perc5],...
                    'k-','LineWidth',1.5);
        end
        hold off 
        cind = cind + 1; 
    end 
end
 
  
  

 
 


