% The script was written by Dr. Qiang Zhang from
% Nanjing Institute of Geography and Limnology,
% Chinese Academy of Sciences on March 26, 2006.
% Financially supported by Alexander von Humboldt Stiftung, Germany
% if you use this script for research, please cite the references below:
%% Reference:
% Gerstengarbe, F.W. and Werner, P.C., 1999. Estimation of the beginning and end of
%% recurrent events within a climate regime. Climate Research 11: 97-107.
% Qiang Zhang, Chunling Liu, Chong-yu Xu, Youpeng Xu, Tong Jiang. Observed trends of annual maximum
% water level and streamflow during past 130 years in the Yangtze River basin, China. Journal of Hydrology, 2006, 324, 255-265.
function [z]=mk(x)
% x为一时间序列(列向量)
% define the time
time=1873+(0:length(x)-1); % please alter the starting year of this series.
n=length(x); %the size of the series
[rows,cols]=size(x);
%---------JUDGEMENT OF TIME SERIES------------
if rows==1
error('x should be a column vector')
end
if cols~=1
error('x should be a column vector')
end
% -----------------COMPUTATION in normal sequence--------------------
for k=2:n
count1=0;
for j=1:k-1
if x(k)>x(j)
count1=count1+1;
end
%junzhi=k*(k-1)/4;
% fangcha=k*(k-1)*(2*k+5)/72;
%z(1,1)=0;
% z(k,1)=(count1-junzhi)/sqrt(fangcha);
end
junzhi1(k,1)=k*(k-1)/4;
fangcha1(k,1)=k*(k-1)*(2*k+5)/72;
zonghe1(1,1)=0;
zonghe1(k,1)=count1;
end
%-------------------COMPUTATION IN ADVERSE SEQUENCY-----------------
x1=flipud(x);
for m=2:n
count2=0;
for i=1:m-1
%bijiao=x1(m)-x1(i);
if x1(m)>x1(i)
count2=count2+1;
end
%junzhi=m*(m-1)/4;
%fangcha=m*(m-1)*(2*m+5)/72;
%z(1,2)=0;
%z(m,2)=(count-junzhi)/sqrt(fangcha);
end
junzhi2(m,1)=m*(m-1)/4;
fangcha2(m,1)=m*(m-1)*(2*m+5)/72;
zonghe2(1,1)=0;
zonghe2(m,1)=count2;
end
%------------------------autocorrelation analysis------------------------
kkk = fix(n/3);
for ii = 0:kkk;
rr(ii+1) = sum(((x(1:n-ii,1)-mean(x))/std(x)).*((x(1+ii:n,1)-mean(x))/std(x)))/(n-ii);
end
%---------------confidence level for autocorrelation coefficient--------
for jj = 0:kkk
pk(jj+1) = length(x(1:n-jj));
end
lower = (-1-1.96*(pk-2).^0.5)./(pk-1);
upper = (-1+1.96*(pk-2).^0.5)./(pk-1);
%----------计算调整系数----------------
adjust=1+2*(rr(2)^(n+1)-n*rr(2)^2+(n-1)*rr(2))/(n*(rr(2)-1)^2);
%-----------调整系数计算完毕-----------
%对原MK值根据自相关系数计算情况进行调值-------------
if rr(2)>upper(2)
%-------------------- Z2 and Z1 value computation ----------------
junzhi1(1,1)=0;
fangcha1(1,1)=0;
z(1,1)=0;
zonghe1=cumsum(zonghe1);
z(2:n,1)=(zonghe1(2:end)-junzhi1(2:end))./sqrt(fangcha1(2:end)*adjust);
%% Z2 value %%
z(1,2)=0;
junzhi2(1,1)=0;
fangcha2(1,1)=0;
zonghe2=cumsum(zonghe2);
z(2:n,2)=(zonghe2(2:end)-junzhi2(2:end))./sqrt(fangcha2(2:end)*adjust);
else
junzhi1(1,1)=0;
fangcha1(1,1)=0;
z(1,1)=0;
zonghe1=cumsum(zonghe1);
z(2:n,1)=(zonghe1(2:end)-junzhi1(2:end))./sqrt(fangcha1(2:end));
%% Z2 value %%
z(1,2)=0;
junzhi2(1,1)=0;
fangcha2(1,1)=0;
zonghe2=cumsum(zonghe2);
z(2:n,2)=(zonghe2(2:end)-junzhi2(2:end))./sqrt(fangcha2(2:end));
end
z(:,2)=-z(:,2);
xx=flipud(z(:,2));
z(:,2)=xx;
%----------------PLOT------------------
%plot(time,z(:,2),'o-k','MarkerSize',3)
%hold on
%plot(time,z(:,1),'o-k','MarkerFaceColor','k',...
% 'MarkerEdgeColor','None',...
% 'MarkerSize',3)
%plot(time,ones(size(time))*1.96, ':k',time,ones(size(time))*-1.96, ':k')
%plot(time,zeros(size(time)),'-k')
%hold off
%xlabel('Time (year)','FontSize',12)
%ylabel('Z value','FontSize',12)