
% This code reads the epidemic dataset deposited in Zenodo and described in
% Marani et al. (2021)

clear all
close all


data=csvread('Epidemic16001945March2021.csv');  % exclude post-WWII data to eliminate vaccines/treatment/policy effects

anno0=min(data(:,1));   % epidemic start year
annomax=max(data(:,1));  % last start year on record

deaths0=zeros(length(data),1);      
duration=zeros(length(data),1); 
population=zeros(length(data),1); 

deaths0=data(:,3); % absolute number of deaths in thousands
duration=data(:,2)-data(:,1)+1; % duration of epidemics in years
population=data(:,4); % global population at start of epidemic in thousands
anno=zeros(annomax-anno0+1,1);


wlength=20; % length of window in years for the estimation of the MEVD
wlength0=1; % length of window in years to estimate the number of events: keep wlength0=1 year
wgev=10; % length of window, in years, for computing block maxima and GEV fitting


thresh=0.001; % threshold for GPD fit determined after several tests

% count the number of all epidemics falling in sliding windows with length
% wlength0

nevent0=zeros(floor((annomax-anno0+1)/wlength0),1);
windowj0=zeros(floor((annomax-anno0+1)/wlength0),1);
avgpop=zeros(floor((annomax-anno0+1)/wlength0),1);

time=data(:,1);

for j=1:floor((annomax-anno0+1)/wlength0) % loop over all events in each year (if wlength0=1)
    windowj0(j,1)=anno0+(j-1)*wlength0;
    nevent0(j,1) = sum(time(:,1)>=anno0+(j-1)*wlength0 & time(:,1)<=anno0+j*wlength0);
    avgpop(j,1) = mean(data(data(:,1)>=anno0+(j-1)*wlength0 & data(:,1)<=anno0+j*wlength0,4));
end


figure
subplot(211)  % plot number of events/year
bar(windowj0(:,1),nevent0(:,1)./wlength0)
xlabel('Year')
ylabel('number of events/year')
axis([1500 2050 0 ceil(max(nevent0(:,1))./wlength0)])
subplot(212)                            % plot global population
plot(windowj0(:,1),avgpop(:,1))
xlabel('Year')
ylabel('Global Population (thousands)')
axis([1600 2050 0 ceil(max(avgpop(:,1)))])


% Compute exceedance probability

data1=data;
deaths=zeros(length(data1),1);
deaths=data1(:,3)./data1(:,4)*1000./(data(:,2)-data(:,1)+1); % relative intensity: number of deaths/global population/duration*1000


% traditional GEV analysis on all magnitudes > 0
% compute block maxima

indices = find(deaths(:,1)<=0);
data1(indices,:) = [];
deaths(indices,:)=[]; 
duration(indices,:)=[]; 
population(indices,:)=[]; 

deathsc=deaths; % copy to keep original values


windowgev=zeros(floor((annomax-anno0+1)/wgev),1);
blockmaxima=zeros(floor((annomax-anno0+1)/wgev),1);


for j=1:floor((annomax-anno0+1)/wgev) % loop over all windows indexed with start year
    windowgev(j,1)=anno0+(j-1)*wgev;
    clear indices
    indices = find(data1(:,2)>=windowgev(j,1) & data1(:,2)<windowgev(j,1)+wgev);
    
    if isempty(indices)
       blockmaxima(j,1)=0;
    else
       blockmaxima(j,1)=max(deaths(indices,1));
    end
        
end

figure   % plot yearly maxima
stem(windowgev(:,1), blockmaxima(:,1),'.')
 xlabel('Year')
 ylabel('Max Epidemic intensity in block (per mil/year)')
 set(gca,'YScale','log')
axis([1600 2050 0 50]) 


sortmax=sort(blockmaxima(:,1),'ascend');
Sdmax=zeros(length(sortmax),1);

for i=1:length(sortmax)
    Sdmax(i,1)=1-i/(length(sortmax)+1);
end

figure   % plot block maxima frequencies and GEV fit

Slowmax=zeros(length(sortmax),1);  % confidence intervals
Supmax=zeros(length(sortmax),1);
nn=length(sortmax);
sigf=1.96; % significance level. 1 = one standard deviation, 1.96= 95% confidence interval
   
for i=1:nn
    Supmax(i,1)=1-((i+0.5*sigf^2)/(nn+sigf+1)-sigf/(nn+sigf^2)*((i*(nn-i)/nn+sigf^2/4))^.5); % Wilson Score estimate of confidence interval
    Slowmax(i,1)=max(10^(-10),1-((i+0.5*sigf^2)/(nn+sigf+1)+sigf/(nn+sigf^2)*((i*(nn-i)/nn+sigf^2/4))^.5)); % https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Wilson_score_interval   
end

loglog(sortmax(:,1),wgev./Sdmax(:,1),'o')
    hold on
loglog(sortmax(:,1),wgev./Supmax,'color',[1,0,0])
loglog(sortmax(:,1),wgev./Slowmax,'color',[0,1,0]) 

xlabel('Epidemic intensity (per mil/year)')
ylabel('Average recurrence time (years)')

% Fit GEV Using ML to relative size in per mil
parmgev = gevfit(sortmax)
pgev = gevcdf(sortmax,parmgev(1),parmgev(2),parmgev(3));


hold on
loglog(sortmax(:,1),wgev./(1-pgev(:,1)))
title('GEV Fit of Epidemic Intensity Block Maxima')

indices = find(deaths(:,1)<=thresh);
data1(indices,:) = [];
deaths(indices,:)=[]; 
duration(indices,:)=[]; 
population(indices,:)=[]; 


length(deaths)


time00=data1(:,1);


nbelow=length(deaths0)-length(deaths); % number of events below threshold thresh, also counting '-999's
pb=nbelow/length(deaths0); % frequency of events below threshold

 
 figure  % plotting all epidemic intensities
 stem(data1(:,1), deaths(:,1),'.')
 xlabel('Year')
 ylabel('Epidemic intensity (per mil/year)')
 set(gca,'YScale','log')
axis([1600 2050 0 50]) 

sorted0=sort(deaths(:,1),'ascend');
Sd0=zeros(length(sorted0),1);

for i=1:length(sorted0)
    Sd0(i,1)=1-i/(length(sorted0)+1);
end

figure

Slow0=zeros(length(sorted0),1);
Sup0=zeros(length(sorted0),1);
nn=length(sorted0);
sigf=1.96; % significance level. 1 = one standard deviation, 1.96= 95% confidence interval
   
for i=1:nn
    Sup0(i,1)=1-((i+0.5*sigf^2)/(nn+sigf+1)-sigf/(nn+sigf^2)*((i*(nn-i)/nn+sigf^2/4))^.5); % Wilson Score estimate of confidence interval
    Slow0(i,1)=max(10^(-4),1-((i+0.5*sigf^2)/(nn+sigf+1)+sigf/(nn+sigf^2)*((i*(nn-i)/nn+sigf^2/4))^.5)); % https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Wilson_score_interval   
end


% perturb observations by applying noise to original data
% generate npert perturbations 
tic
npert=10000;  % number of realizations
amplitude=1.; % amplitude of perturbation as a fraction
parm1i=zeros(npert,1);
parm2i=zeros(npert,1);

for i=1:npert
    deathi=zeros(length(deathsc),1);
    deathi(:,1)=((rand(length(deathsc(:,1)),1)).*amplitude+1-amplitude/2).*deathsc(:,1);  % apply perturbation
    
    indices = find(deathi(:,1)<=thresh);
    deathi(indices,:)=[]; 
    
      sorted0i=sort(deathi(:,1),'ascend');
  
 % Fit GPD Using ML to relative size in per mil for each perturbed
 % realization
parmhat = gpfit(sorted0i-thresh);

p = gpcdf(sorted0i-thresh,parmhat(1),parmhat(2),0);

parm1i(i,1)=parmhat(1);
parm2i(i,1)=parmhat(2);

loglog(sorted0i(:,1),1-(pb+(1-pb)*p(:,1)),'color',[0.5,.5,.5,0.006])  % draw GPD fit for each realization in grey
 

    hold on
    
end


loglog(sorted0(:,1),1-(pb+(1-pb)*(1-Sd0(:,1))),'ok','MarkerSize',10,'MarkerFaceColor', 'w')  % plot empirical frequencies of unperturbed dataset
    hold on
loglog(sorted0(:,1),1-(pb+(1-pb)*(1-Sup0)),'color',[0,0,0])   % plot uncertainty interval around empirical frequencies
loglog(sorted0(:,1),1-(pb+(1-pb)*(1-Slow0)),'color',[0,0,0]) 


xlabel('Epidemic intensity (per mil/year)')
ylabel('Exceedance Probability')


% Fit GPD Using ML to relative intensity in per mil/year to the unperturbed
% dataset

parmhat = gpfit(sorted0-thresh)

p = gpcdf(sorted0-thresh,parmhat(1),parmhat(2),0);

hold on
loglog(sorted0(:,1),1-(pb+(1-pb)*p(:,1)),'LineWidth',2,'color',[1,0,0])  %plot GPD fit to unperturbed dataset in red

% compute some statistics on GPD probabilities obtained from perturbed
% realizations. these statistics are not plotted but printed out

meanp=zeros(length(sorted0),1); % initialize mean of GPD probabilities
stdevp=zeros(length(sorted0),1); % initialize std dev of GPD probabilities
meanparm1=zeros(1,1); % initialize mean of GPD shape parameter
meanparm2=zeros(1,1); % initialize mean of GPD scale parameter
stdparm1=zeros(1,1); % initialize std dev of GPD shape parameter
stdparm2=zeros(1,1); % initialize std dev of GPD scale parameter
pipert=zeros(length(sorted0),npert);
pup=zeros(length(sorted0),1);
pdown=zeros(length(sorted0),1);

for i=1:npert  % now compute GPD values for mean observations and for all perturbations and compute their mean and standard dev
   
   pipert(:,i)=gpcdf(sorted0-thresh,parm1i(i,1),parm2i(i,1),0);
   meanp=meanp+gpcdf(sorted0-thresh,parm1i(i,1),parm2i(i,1),0);
   stdevp=stdevp+gpcdf(sorted0-thresh,parm1i(i,1),parm2i(i,1),0).^2;
   
end

for i=1:length(sorted0) % compute 2.5 and 97.5 percentiles across perturbed samples
pup(i,1)=prctile(pipert(i,:),2.5);
pdown(i,1)=prctile(pipert(i,:),97.5);
end

meanp=meanp./npert;  % mean GPD probabilities
stdevp=sqrt((stdevp./npert-meanp.^2)); % standard deviation of above GPD probability
meanparm1=mean(parm1i)
meanparm2=mean(parm2i)
stdparm1=std(parm1i)
stdparm2=std(parm2i)


axis([.001 10 0.001 1]) 


% Build the Metastatistical Extreme Value Distribution

wj=zeros(floor((annomax-anno0+1)/wlength),1);


mevd=zeros(length(sorted0),floor((annomax-anno0+1)/wlength));
mevdmin=zeros(length(sorted0),floor((annomax-anno0+1)/wlength));
mevdmax=zeros(length(sorted0),floor((annomax-anno0+1)/wlength));

nwindow=floor((annomax-anno0+1)/wlength); % number of windows that fit in the record
j0=length(nevent0)-wlength*nwindow  % start from here to have an integer number of windows in the record

for i=1:length(sorted0) % compute MEVD values for all observed quantiles
for k=1:j0-1 %  works properly only if wlength0=1 year

     mevd(i,1)=mevd(i,1)+(pb+(1-pb)*meanp(i,1))^nevent0(k,1); % use average GPD computed across perturbed epidemic intensities

     mevdmax(i,1)=mevdmax(i,1)+(pb+(1-pb)*(pup(i,1)))^nevent0(k,1); % 97.5 percentile for ordinary distribution
     mevdmin(i,1)=mevdmin(i,1)+(pb+(1-pb)*(pdown(i,1)))^nevent0(k,1); % 2.5 percentile of ordinary distribution

end
mevd(i,1)=mevd(i,1)/(j0-1);
mevdmin(i,1)=mevdmin(i,1)/(j0-1);
mevdmax(i,1)=mevdmax(i,1)/(j0-1);


end
wj(1,1)=anno0+j0+wlength;


for j=2:nwindow % loop over all windows indexed with start year
    wj(j,1)=anno0+j0+j*wlength;
    
    % MEVD-GPD
    for i=1:length(sorted0) % compute MEVD values for all observed quantiles
        for k=1:wlength % compute a MEVD as the mean over all years in each window of wlength years
%            mevd(i,j)=mevd(i,j)+(pb+(1-pb)*p(i,1))^nevent0(j0+(j-1)*wlength+k,1); % use GPD from avg epidemic intensities 
             mevd(i,j)=mevd(i,j)+(pb+(1-pb)*meanp(i,1))^nevent0(j0+(j-1)*wlength+k,1); % use average GPD computed across perturbed epidemic intensities
%             mevdmax(i,j)=mevdmax(i,j)+(pb+(1-pb)*(meanp(i,1)+stdevp(i,1)))^nevent0(j0+(j-1)*wlength+k,1);
%             mevdmin(i,j)=mevdmin(i,j)+(pb+(1-pb)*(meanp(i,1)-stdevp(i,1)))^nevent0(j0+(j-1)*wlength+k,1); % subtract a standard deviation
             mevdmax(i,j)=mevdmax(i,j)+(pb+(1-pb)*(pup(i,1)))^nevent0(j0+(j-1)*wlength+k,1); % 2.5 percentile
             mevdmin(i,j)=mevdmin(i,j)+(pb+(1-pb)*(pdown(i,1)))^nevent0(j0+(j-1)*wlength+k,1); % 97.5 percentile

        end
        mevd(i,j)=mevd(i,j)/wlength;
        mevdmin(i,j)=mevdmin(i,j)/wlength;
        mevdmax(i,j)=mevdmax(i,j)/wlength;

     end
end

figure  % plot Extreme value distribution using number of events in the last time slice on record
plot(sorted0(:,1),1-mevd(:,floor((annomax-anno0+1)/wlength)),'r','LineWidth',2)
hold on
xxx=[0 sorted0(1,1)];
yyy=[1 1-mevd(1,floor((annomax-anno0+1)/wlength))];
plot(xxx,yyy,'r','LineWidth',2)


sort2 = [sorted0', fliplr(sorted0')];        % repeat x values
yym = [(1-mevdmin(:,floor((annomax-anno0+1)/wlength)))', fliplr((1-mevdmax(:,floor((annomax-anno0+1)/wlength)))')];   % vector of upper & lower boundaries
patch(sort2,yym,[0.5 0.5 0.5],'FaceAlpha',0.3,'LineStyle','none')    % fill area defined by x & yy in grey
set(gca, 'XScale','linear', 'YScale','log')


xlabel('Epidemic intensity (per mil year^{-1})')
ylabel('Yearly probability of exceedance')

    figure
    stem(time00(:,1),deaths(:,1).*data1(:,4),'.')
 
    xlabel('Year')
    ylabel('number of deaths/year')
    
    

mevd0=zeros(floor((annomax-anno0+1)/wlength),1);
mevd0plus=zeros(floor((annomax-anno0+1)/wlength),1);
mevd0minus=zeros(floor((annomax-anno0+1)/wlength),1);

windj=zeros(floor((annomax-anno0+1)/wlength),1);
parmhat1=1.3981; % values of GPD parameters from 1600-1945
parmhat2=0.0113;

meanp0=zeros(1,1); % initialize mean of GPD probability of prescribed event
stdevp0=zeros(1,1); % initialize std dev of GPD probability of prescribed event
pi=zeros(npert,1);

for i=1:npert  % now compute GPD value for mean observations and for all perturbations 
    
   pi(i,1)=gpcdf(intensity-thresh,parm1i(i,1),parm2i(i,1),0); % store all values of probab. for all perturbations
   meanp0=meanp0+gpcdf(intensity-thresh,parm1i(i,1),parm2i(i,1),0); % compute avg
   stdevp0=stdevp0+gpcdf(intensity-thresh,parm1i(i,1),parm2i(i,1),0).^2;  % compute stdev
   
end

meanp0=meanp0/npert % mean GPD probability of prescribed event across random perturbations
stdevp0=sqrt((stdevp0/npert-meanp0^2))  % standard deviation of the above

% compute 2.5 and 97.5 percentiles of the probab from all perturbations
p25=prctile(pi(:,1),2.5)
p975=prctile(pi(:,1),97.5)

toc                                  