% This is version 2 of the code that analyzes global epidemic data as in
% Marani et al., PNAS, 2021
% Changes with respect to the initial version of the code are as follows:
% 1. An error is corrected in the way the code computed the number of epidemics/year.
%    This error attributed to each year the number of events occurred in the
%    following year as well. In turn, this resulted in inflated values of
%    the probability of occurrence (or underestimated values of the associated return
%    times);
% 2. The new code version now analyzes a full dataset 1600-2019 in the file
%    "Epidemic16002020March2021". Data in the period 1600-1945 are used to
%    estimate the GPD distribution of epidemic intensities and are the same
%    as those contained in the file "Epidemic16001945March2021.csv" .
%    The more recent period 2000-2019 is used to estimate
%    current and projected probabilities of epidemic occurrence. 
%    
% The code produces all figures published in Marani et al., PNAS (2021),
% and uses the methods described therein. An erratum is, at the time of
% this publication, under consideration by PNAS.

clear all
close all

data=csvread('Epidemic16002020March2021.csv'); % all data since 1600

%anno0=min(data(:,1));   % epidemic start year
anno0=1600; % start analyses from 1600
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. MEVD is, by definition,
            % the probability distribution of the maximum intensity in a window
            % of wlength years. This is adjusted for at the end, to
            % calculate the MEVD for yearly maxima.

wlength0=1; % length of window in years to record the number of events. This code is now tested
            % for wlength0=1 year. The idea of potentially changing wlength0 is that one may compute the number
            % of events, say, for 10-year windows, to possibly reduce
            % uncertainty/variability in these estimates. This, however, has not been tested, and the results
            % in Marani et al. (2021) assume wlength0=1.

wgev=15; % length of window, in years, for computing block maxima for GEV fitting


thresh=0.001;   % this threshold on relative intensity, needed in the GPD underlying the MEVD,
                % was set after experimenting 
                % with several threshold values and is the value used in
                % the analyses described in Marani et al. (2021)

% count the number of epidemics from all epidemics, including events for which we have no data
% ('-999') and those <10k deaths ('0'), falling in each sliding window with length
% wlength0. This is now only tested for wlength0=1 year.

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)   % compute number of events and average global population
                                            % in each sliding window of length wlength0

    windowj0(j,1)=anno0+(j-1)*wlength0; % keep track of window start year

    nevent0(j,1) = sum(time(:,1)>=anno0+(j-1)*wlength0 & time(:,1)<anno0+j*wlength0); % compute # of events/year
                                                                                     % it is used in the following to compute MEVD

    avgpop(j,1) = mean(data(data(:,1)>=anno0+(j-1)*wlength0 & data(:,1)<anno0+j*wlength0,4)); % compute avg 
                                                                                            % global population
end

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

% Compute exceedance probability
% only keep epidemics for which there is quantitative information, by
% eliminating those for which # of deaths is unknown ('-999') and those for
% which we only know # of deaths was < 10k ('0')

% The exceedance probability, as described in Marani et al 2021, is
% computed only on 1600-1945 data

data1=data; % copy data to data1, to keep original record. events with unknown 
            % # deaths will be removed from data1

deaths=zeros(length(data1),1);

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


% First perform traditional GEV analysis on all magnitudes > 0

indices = find(deaths(:,1)<=0); % find events for which # of deaths is not known ('-999') 
                                % and those with <10k deaths (marked as '0')
data1(indices,:) = [];  % eliminate, for the purpose of GEV and, later, GPD parameter estimation,
                        % all epidemics for which relative intensity cannot
                        % be computed.
deaths(indices,:)=[]; 
duration(indices,:)=[]; 
population(indices,:)=[]; 

deathsc=deaths; % copy to keep original values of relative intensity: values will be needed for MEVD-GPD analysis 

% now only select data in 1600-1945

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

% find block maxima on blocks of length wgev for GEV analysis
for j=1:floor((annomax1-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 2 - Plotting block maxima for wgev-long blocks
figure
stem(windowgev(:,1), blockmaxima(:,1),'.')
 xlabel('Year')
 ylabel('15-year block maxima for GEV analysis (per mil/year)')
 set(gca,'YScale','log')
axis([1600 2050 0 50]) 

% sort block maxima to compute empirical exceedance frequency
sortmax=sort(blockmaxima(:,1),'ascend');
Sdmax=zeros(length(sortmax),1);

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


% Compute confidence interval limits using Wilson's score, e.g. see
% https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Wilson_score_interval
Slowmax=zeros(length(sortmax),1);
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


% Figure 3 - Plot empirical estimates of return periods for GEV visual goodness of fit
figure
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)); % get GEV cdf values for all observed intensities

% Plot GEV return periods to compare with empirical values
hold on
loglog(sortmax(:,1),wgev./(1-pgev(:,1)))
title('GEV Fit of Epidemic Intensity Block Maxima')

% GEV analysis concluded
% GPD fit starts
% Limit GPD fit to data in 1600-1945
indices1 = find(data1(:,1)>1945);
data1(indices1,:) = [];  % these are the original "raw" data, with years, number of deaths, global population
deaths(indices1,:)=[]; % remove values below threshold in the relative intensity data 
duration(indices1,:)=[]; % same for durations
population(indices1,:)=[]; % same for global population number

ntot=length(deaths0); % number of events above and below threshold (events with unknown # of deaths - '-999' - 
                     % and those with <10 k deaths ('0') are counted as
                     % below-threshold. See comments elsewhere in the code.

% Apply now the threshold to the data before fitting the GPD distribution
indices = find(deaths(:,1)<=thresh);
data1(indices,:) = [];  % these are the original "raw" data, with years, number of deaths, global population
deaths(indices,:)=[]; % remove values below threshold in the relative intensity data 
duration(indices,:)=[]; % same for durations
population(indices,:)=[]; % same for global population number


nbelow=ntot-length(deaths); % number of events below threshold thresh,
                                       % assuming '-999's and events below
                                       % 10k (i.e. '0') are all below the
                                       % selected threshold. Since we do
                                       % not know the actual # of deaths we
                                       % cannot determine if they are
                                       % <threshold. One can argue
                                       % that these lesser-known events
                                       % were not major ones in terms of
                                       % deaths.

pb=nbelow/ntot; % frequency of events below threshold


 % Figure 4 - Plot the sequence of epidemic intensities just for the above-threshold
 % events

 figure
 stem(data1(:,1), deaths(:,1),'.')
 xlabel('Year')
 ylabel('Epidemic intensity (per mil/year)')
 title('Sequence of above-threshold events')
 set(gca,'YScale','log')
axis([1600 2050 0 50]) 

% Sort relative intensities for thresholded data
sorted0=sort(deaths(:,1),'ascend');
Sd0=zeros(length(sorted0),1);

% Compute exceedance frequency of above-threshold events
for i=1:length(sorted0)
    Sd0(i,1)=1-i/(length(sorted0)+1);
end


% Figure 5 - Prepare empirical frequencies of all above-threshold events to be later
% compared with GPD exceedance probability distribution
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 Wilson Score method
   
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


% Sensitivity to uncertainty in the data: perturb observations by applying noise to original data
% generate npert perturbations 

npert=10000; % number of perturbation realizations
amplitude=1.; % amplitude of perturbation. If equal to 1 we are generating +/- 50% random perturbations
              % to original number relative intensities

parm1i=zeros(npert,1); % we will fit GPDs to perturbed datasets
parm2i=zeros(npert,1); % these vectors will store the resulting GPD parameters

for i=1:npert
    deathi=zeros(length(deathsc),1); % Here original data, before thresholding, are considered again to
    deathi(:,1)=((rand(length(deathsc(:,1)),1)).*amplitude+1-amplitude/2).*deathsc(:,1); % apply 
                                                                              % the random perturbation.
    indices = find(deathi(:,1)<=thresh);  % The threshold is reapplied after the perturbation. Some events
    deathi(indices,:)=[];                 % that were below threshold may now be above, and viceversa
    
     sorted0i=sort(deathi(:,1),'ascend'); % sorting the perturbed and thresholded dataset
  
% Fit GPD Using ML to relative intensities from the i-th perturbed dataset 
parmhat = gpfit(sorted0i-thresh);
p = gpcdf(sorted0i-thresh,parmhat(1),parmhat(2),0);

parm1i(i,1)=parmhat(1); % store i-th parameter values
parm2i(i,1)=parmhat(2);

% plot in light grey resulting GPD exceedance probability in graph with empirical exceedance frequencies
loglog(sorted0i(:,1),1-(pb+(1-pb)*p(:,1)),'color',[0.5,.5,.5,0.006])
 
hold on
    
end

% Plot empirical exceedance frequencies with closed circles. They are
% plotted here, after the previous for loop, so that the circles are on top
% of the grey lines.
loglog(sorted0(:,1),1-(pb+(1-pb)*(1-Sd0(:,1))),'ok','MarkerSize',10,'MarkerFaceColor', 'w')

loglog(sorted0(:,1),1-(pb+(1-pb)*(1-Sup0)),'color',[0,0,0])
loglog(sorted0(:,1),1-(pb+(1-pb)*(1-Slow0)),'color',[0,0,0]) 

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


% Fit GPD to the un-perturbed data using ML to relative intensity in per mil/year
parmhat = gpfit(sorted0-thresh);
p = gpcdf(sorted0-thresh,parmhat(1),parmhat(2),0);


% Do not plot the GPD exceedance probability fitted to the original data in a thick red line, uncomment if useful 
%loglog(sorted0(:,1),1-(pb+(1-pb)*p(:,1)),'LineWidth',2,'color',[1,0,0])


title('GPD Fit of Epidemic Intensity')
axis([.001 10 0.001 1]) 


% Now compute average values of exceedance probabilities from all
% dataset perturbations
meanp=zeros(length(sorted0),1); % initialize mean of GPD probabilities
stdevp=zeros(length(sorted0),1); % initialize std dev of GPD probabilities

pipert=zeros(length(sorted0),npert);
pup=zeros(length(sorted0),1);
pdown=zeros(length(sorted0),1);

for i=1:npert  % now compute GPD cdf values for all perturbations and compute their mean and standard dev
   
   pipert(:,i)=gpcdf(sorted0-thresh,parm1i(i,1),parm2i(i,1),0); % cdf value from i-th perturbation
   meanp=meanp+gpcdf(sorted0-thresh,parm1i(i,1),parm2i(i,1),0); % compute the mean cdf value across perturbations
   stdevp=stdevp+gpcdf(sorted0-thresh,parm1i(i,1),parm2i(i,1),0).^2; % compute std dev of cdf across perturbations
   
end

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

meanp=meanp./npert;  % mean GPD probabilities across all perturbations
stdevp=sqrt((stdevp./npert-meanp.^2)); % standard deviation of above GPD probability


% Plot, in red, the GPD exceedance probability from averaging across all perturbed datasets 
% Note: this is 'centered' within the cloud of cdfs from all perturbed
% realizations: this is the line that currently appears in Figure 1 in the
% paper.
% The GPD cdf obtained from fitting the original data does not fall in the
% middle of the cloud of cdfs, due to the nonlinearities in the
% randomization and GPD fit process.
loglog(sorted0(:,1),1-(pb+(1-pb)*meanp(:,1)),'LineWidth',1,'color',[1,0,0])

% In the lines above the mean and standard deviation of cdf values are used to
% quantify uncertainty. If uncertainty should be quantified through the mean and standard deviations of the
% GPD parameters uncomment and use lines below.
%meanparm1=mean(parm1i); % Compute mean GPD parameters and standard deviations 
                        % from all parameter values from all perturbations
%meanparm2=mean(parm2i); 
%stdparm1=std(parm1i);
%stdparm2=std(parm2i);


% GPD Estimation completed. Now build the MEVD using the GPD and the number of events/year

% First compute MEVD exceedance probabilities for observed events
% Note that the GPD was fitted to data between 1600 and 1945. The MEVD will
% next be constructed using the number of events/year for all years on
% record, through 2018.

wj=zeros(floor((annomax-anno0+1)/wlength)+1,1); % Initialize vector that will contain start year
                                              % for each window, of length wlength, needed to
                                              % compute MEVD

mevd=zeros(length(sorted0),floor((annomax-anno0+1)/wlength)+1); % initialize vectors that will contain MEVD Exceedance Probability
mevdmin=zeros(length(sorted0),floor((annomax-anno0+1)/wlength)+1); % and min and max from perturbed data
mevdmax=zeros(length(sorted0),floor((annomax-anno0+1)/wlength)+1);

nwindow=floor((annomax-anno0+1)/wlength); % number of whole windows that fit in the record
j0=length(nevent0)-wlength*nwindow;  % start at j0 + 1, rather than at the beginning,
                                    % leaving out a shorter first window, with length=j0 

for i=1:length(sorted0) % compute MEVD values for all observed quantiles

for k=1:j0   % runs through all years within the first time window of length j0
             % has only been tested in the case in which nevent0 contains yearly numbers of events
             % MEVD is computed as the average of the values of cdf^n_events from each year 
             % in the first time slice, or window, of length j0 

%     The next line uses the GPD fitted to the unperturbed data, i.e. the
%     data as they are, to build the MEVD
%     mevd(i,1)=mevd(i,1)+(pb+(1-pb)*p(i,1))^nevent0(k,1); % use GPD fitted to 1600-1945 epidemic intensities 

% The next line uses the mean cdf value across all perturbed datasets to
% build the MEVD
     mevd(i,1)=mevd(i,1)+(pb+(1-pb)*meanp(i,1))^nevent0(k,1); % use average GPD computed across perturbed epidemic intensities

% The next two lines use the standard deviation of the GPD cdf around the
% mean, computed across all perturbed dataset realizations, to quantify
% uncertainty in the MEVD
%     mevdmax(i,1)=mevdmax(i,1)+(pb+(1-pb)*(meanp(i,1)+stdevp(i,1)))^nevent0(k,1); % add a standard deviation
%     mevdmin(i,1)=mevdmin(i,1)+(pb+(1-pb)*(meanp(i,1)-stdevp(i,1)))^nevent0(k,1); % subtract a standard deviation

% The next two lines use the 2.5 and 97.5 percentile of the GPD cdf 
% computed across all perturbed dataset realizations, to quantify
% uncertainty in the MEVD
     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

% store in the appropriate vectors the MEVD computed for the first
% shorter window of length=j0-1
mevd(i,1)=mevd(i,1)/j0;
mevdmin(i,1)=mevdmin(i,1)/j0;
mevdmax(i,1)=mevdmax(i,1)/j0;


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


for j=2:nwindow+1 % loop over all windows, after the initial shorter window, indexed with start year
                % has only been tested in the case in which nevent0 contains yearly numbers of events
                % MEVD is computed as the average of the values of cdf^n_events from each year 
                % in the time slice, or window, of length wlength 
    wj(j,1)=anno0+j0+(j-1)*wlength; % last year in the window, used to index it
    
    % 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 the j-th window of wlength years

%     The next line uses the GPD fitted to the unperturbed data, i.e. the
%     data as they are, to build the MEVD
%            mevd(i,j)=mevd(i,j)+(pb+(1-pb)*p(i,1))^nevent0(j0+(j-2)*wlength+k,1); % use GPD fitted to original,
                                                                                  % unperturbed intensities

% The next line uses the mean cdf value across all perturbed datasets to
% build the MEVD                                                                                  
           mevd(i,j)=mevd(i,j)+(pb+(1-pb)*meanp(i,1))^nevent0(j0+(j-2)*wlength+k,1); % use average GPD computed across perturbed epidemic intensities

% The next two lines use the standard deviation of the GPD cdf around the
% mean, computed across all perturbed dataset realizations, to quantify
% uncertainty in the MEVD
%             mevdmax(i,j)=mevdmax(i,j)+(pb+(1-pb)*(meanp(i,1)+stdevp(i,1)))^nevent0(j0+(j-2)*wlength+k,1);
%             mevdmin(i,j)=mevdmin(i,j)+(pb+(1-pb)*(meanp(i,1)-stdevp(i,1)))^nevent0(j0+(j-2)*wlength+k,1); % subtract a standard deviation

% The next two lines use the 2.5 and 97.5 percentile of the GPD cdf 
% computed across all perturbed dataset realizations, to quantify
% uncertainty in the MEVD
             mevdmax(i,j)=mevdmax(i,j)+(pb+(1-pb)*(pup(i,1)))^nevent0(j0+(j-2)*wlength+k,1); % 2.5 percentile
             mevdmin(i,j)=mevdmin(i,j)+(pb+(1-pb)*(pdown(i,1)))^nevent0(j0+(j-2)*wlength+k,1); % 97.5 percentile

        end % end loop over years ('k') within window

 % store in the appropriate vectors the MEVD computed for the j-th window
        mevd(i,j)=mevd(i,j)/wlength;
        mevdmin(i,j)=mevdmin(i,j)/wlength;
        mevdmax(i,j)=mevdmax(i,j)/wlength;

    end % end loop over observed quantile ('i')
end % end loop over time slices (i.e. windows indexed by 'j')

% Plot exceedance probability from MEVD computed over the last time slice
% of length wlength (2000-2019 if wlength=20 yrs)
figure
plot(sorted0(:,1),1-mevd(:,ceil((annomax-anno0+1)/wlength)),'r','LineWidth',2)
hold on
xxx=[0 sorted0(1,1)];
yyy=[1 1-mevd(1,ceil((annomax-anno0+1)/wlength))];
plot(xxx,yyy,'r','LineWidth',2)
semilogy(sorted0(:,1),1-mevdmin(:,ceil((annomax-anno0+1)/wlength)),'k--','LineWidth',2)
semilogy(sorted0(:,1),1-mevdmax(:,ceil((annomax-anno0+1)/wlength)),'k--','LineWidth',2)

sort2 = [sorted0', fliplr(sorted0')];        % repeat x values
yym = [(1-mevdmin(:,ceil((annomax-anno0+1)/wlength)))', fliplr((1-mevdmax(:,ceil((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')
title('Yearly Probability of Exceedance for the most recent time slice (2000-2019)')

% Computation of MEVD for observed quantiles has ended

% Now, the MEVD for an event of prescribed intensity 'intensity' is
% computed and plotted over time

% initialize new vectors for MEVD computation
mevd0=zeros(floor((annomax-anno0+1)/wlength)+1,1);
mevd0plus=zeros(floor((annomax-anno0+1)/wlength)+1,1);
mevd0minus=zeros(floor((annomax-anno0+1)/wlength)+1,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Choose the intensity of the event to be plotted
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%intensity=0.32439073;  % intensity of COVID-19 as of 3 June 2021
intensity=2.63; % 8 times the COVID-19 intensity to determine the unconstrained avg recurrence interval
%intensity=5.69405149557821 % intensity of Spanish Flu

%intensity=1000;  % intensity that would erase humanity

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


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  % compute GPD value of selected event as a mean across 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 probability of selected event from all perturbations
p25=prctile(pi(:,1),2.5);
p975=prctile(pi(:,1),97.5);


% now plot probability of prescribed event across time

% First, shorter, window, of length j0 as before
for k=1:j0 % compute a MEVD for prescribed event - works properly only if wlength0=1 year
                                                            % use GPD fitted to unperturbed intensities 
%        mevd0(1,1)=mevd0(1,1)+(pb+(1-pb)*gpcdf(intensity-thresh,parmhat(1),parmhat(2),0))^nevent0(k,1); 
        
        mevd0(1,1)=mevd0(1,1)+(pb+(1-pb)*meanp0)^(nevent0(k,1)); % use average GPD computed across perturbed epidemic intensities      

        mevd0minus(1,1)=mevd0minus(1,1)+(pb+(1-pb)*(p25))^(nevent0(k,1)); % this uses 2.5 quantile for lower confidence limit
        mevd0plus(1,1)=mevd0plus(1,1)+(pb+(1-pb)*(p975))^(nevent0(k,1)); % this uses 97.5 quantile for upper confidence limit
end

mevd0(1,1)=mevd0(1,1)/j0;
mevd0minus(1,1)=mevd0minus(1,1)/j0;
mevd0plus(1,1)=mevd0plus(1,1)/j0;
windj(1,1)=anno0+j0-1;

% all other time slices/windows of length wlength
for j=2:nwindow+1 % loop over all windows indexed with end year
    windj(j,1)=anno0+j0+(j-1)*wlength-1;


for k=1:wlength % compute a MEVD for prescribed event
                            % use GPD fitted to unperturbed intensities
%        mevd0(j,1)=mevd0(j,1)+(pb+(1-pb)*gpcdf(intensity-thresh,parmhat(1),parmhat(2),0))^nevent0(j0+(j-2)*wlength+k,1);

         mevd0(j,1)=mevd0(j,1)+(pb+(1-pb)*meanp0)^(nevent0(j0+(j-2)*wlength+k,1)); % use average GPD computed across perturbed epidemic intensities        
         
         mevd0minus(j,1)=mevd0minus(j,1)+(pb+(1-pb)*((p25)))^(nevent0(j0+(j-2)*wlength+k,1));
         mevd0plus(j,1)=mevd0plus(j,1)+(pb+(1-pb)*(p975))^(nevent0(j0+(j-2)*wlength+k,1));
end
mevd0(j,1)=mevd0(j,1)/wlength;
mevd0minus(j,1)=mevd0minus(j,1)/wlength;
mevd0plus(j,1)=mevd0plus(j,1)/wlength;
end

gevtr(1,1)=235.1; % return time from GEV for Spanish Flu event
gevtr(2,1)=235.1;
t(1,1)=min(windj);
t(2,1)=max(windj);


figure
plot(windj(:,1),(1-mevd0).*100,'r-','LineWidth',2) % plot yearly probability of prescribed event according to MEVD-GPD
hold on
plot(windj(:,1),(1-mevd0minus).*100,'k','LineWidth',1) % plot yearly probability of prescribed event according to MEVD-GPD
plot(windj(:,1),(1-mevd0plus).*100,'k','LineWidth',1) % plot yearly probability of prescribed event according to MEVD-GPD

plot(t(:,1),1./gevtr(:,1).*100,'k--','LineWidth',2)

wind2 = [windj', fliplr(windj')];        % repeat x values
yy = [(1-mevd0minus').*100, fliplr((1-mevd0plus').*100)];   % vector of upper & lower boundaries
patch(wind2,yy,[0.5 0.5 0.5],'FaceAlpha',0.3,'LineStyle','none')    % fill area defined by x & yy in grey


xlabel('Year')
ylabel('Yearly probability (%)')

                              
disp('GPD distribution of relative epidemic intensity based on 1600-1945 data:')
X = ['Frequency of events above the threshold=',num2str(1-pb)];
disp(X)

X = ['GPD parameters for unperturbed relative intensity: shape=',num2str(parmhat(1)),', scale=',num2str(parmhat(2))];
disp(X)

X = ['Asymptotic power-law exponent=', num2str(1/parmhat(1))];
disp(X)

disp('Selected Event to be specified in the code: Probability estimates based on 2000-2019 number of occurrences')

X = ['Probability of Exceedance of Selected Event (mean)=',num2str(1-mevd0(nwindow+1,1)),'; corresponding Tr=',num2str(1/(1-mevd0(nwindow+1,1)))];
disp(X)

X = ['Probability of Exceedance of Selected Event (2.5% quantile)=',num2str(1-mevd0minus(nwindow+1,1)),'; corresponding Tr=',num2str(1/(1-mevd0minus(nwindow+1,1)))];
disp(X)

X = ['Probability of Exceedance of Selected Event (97.5% quantile)=',num2str(1-mevd0plus(nwindow+1,1)),'; corresponding Tr=',num2str(1/(1-mevd0plus(nwindow+1,1)))];
disp(X)

X = ['Probability of Occurrence in a 100 yr period=',num2str(1-mevd0(nwindow+1,1)^100),'; 95% confidence limits:',num2str(1-mevd0minus(nwindow+1,1)^100),',',num2str(1-mevd0plus(nwindow+1,1)^100)];
disp(X)

