function SCEIExample
% Compute the  Standardized Compound Event Indicator (SCEI) of dry and
% hot condition based on precipitation and temperature
% Please check the paper for detail: 
% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% Hao Z. et al. (2019). A monitoring and prediction system for compound dry and
% hot events. Environ. Res. Lett. , in press
% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 clc;clear all;close all
PD=rand(468,1); % Monthly precipitation (39 year)
TD=rand(468,1); % Monthly temperature(39 year)
% Time scale
ts=3;
% Compute accumulation
APD=zeros(size(PD)).*nan;
ATD=zeros(size(TD)).*nan;
for i=1:(length(PD)-ts+1)
    APD(ts+i-1)=mean(PD(i:(i+ts-1)));
    ATD(ts+i-1)=mean(TD(i:(i+ts-1)));
end

% SCEI computation of each month

SCEI=zeros(size(PD)).*nan;
for k=1:12
    clear dp dt
    if k<ts
        dp=APD(k+12:12:end);
        dt=ATD(k+12:12:end);
    else
        dp=APD(k:12:end);
        dt=ATD(k:12:end);
    end
    
    jp=empdis2([dp dt]);
    
    if k<ts
        SCEI(k+12:12:end)=norminv(empdis(jp));
    else
        SCEI(k:12:end)=norminv(empdis(jp));
    end
end

n=length(SCEI);
figure(1);set(gcf,'position',[100,-100,900,500])
plot(1:n,SCEI,'b','LineWidth',2);hold on;

hold on; plot(0:(n+1),zeros(n+2,1),'k-.')
xlabel('Year','FontSize',12);ylabel('Index Value','FontSize',12);
ylim([-2.5 2.5])
xk=[1980:3:2018];for i=1:length(xk) xl{i}=num2str(xk(i)); end
for i=1:length(xk) xt(i)=(xk(i)-1980).*12+1; end
set(gca,'xtick',xt,'xticklabel',xl)
title('Standardized Compound Event Indicator (SCEI)','FontSize',12)



function y=empdis(d)
%  Comput the univariate empirical probability
n=length(d);
bp=zeros(n,1);
for i=1:n
    bp(i,1)=sum(d(:,1)<=d(i,1));
end

% y=(bp)./(n+1); Weibull
% Gringorten plotting position
y=(bp-0.44)./(n+0.12);

function y=empdis2(d)
% Comput the bivariate empirical probability of dry and hot 
% P(X¡Üx,Y>y)
n=length(d);
bp=zeros(n,1);
for i=1:n
    td=zeros(n,3);
    td(find(d(:,1)<=d(i,1)),1)=1;
    td(find(d(:,2)>d(i,2)),2)=1;
    td(:,3)=td(:,1).*td(:,2);
    bp(i,1)=sum(td(:,3));
end
% y=(bp)./(n+1); Weibull 
% Gringorten plotting position
y=(bp-0.44)./(n+0.12);