function [tau,z] = kendall_tau(x1,x2)

%KENDALL_TAU computes Kendall's tau statistic for two distributions.

%INPUTS:  x1 = first n-vector of data
%         x2 = second n-vector of data

%OUTPUT: tau = Kendall's tau
%          z = z-statistic [dist N(0,1) under null hyp of independence]

%Reference: Kendall, M.G. (1962) RANK CORRELATION METHODS, 3rd Ed., Hafner Pub, New York.

%NOTE: This program assumes that neither x1 or x2 are constant vectors( i.e.,
%      of the form a*ones(n,1).

n = length(x1);

if n ~= length(x2)
   
   error('vectors must be of same length');
   
end


%*******************************

%First compute ranks

%*******************************

L1 = [[1:n]',x1];

L1 = sortrows(L1,2);

R1 = zeros(n,1);

L2 = [[1:n]',x2];

L2 = sortrows(L2,2);

R2 = zeros(n,1);

%Now convert x1 to proper ranks

rows = L1(:,1);

R = [1:n]';

L = L1(:,2);

i = 1;

D1 = n*(n-1)/2;

while i < n
   
   ctr = i;
   
   j = 1;
   
   while (i + j <= n)&(L(i+j) == L(i))
      
      ctr = ctr + (i + j);
      
      j = j + 1;
      
   end
   
   if j > 1
      
      R(i:i+j-1) = ctr/j;  %Revise ranks
      
      D1 = D1 - j*(j-1)/2; %Revise D1
      
      i = i + j;
      
   else
      
      i = i + 1;
      
   end
   
end

R1(rows) = R;


%Next convert x2 to proper ranks

rows = L2(:,1);

R = [1:n]';

L = L2(:,2);

i = 1;

D2 = n*(n-1)/2;

while i < n
   
   ctr = i;
   
   j = 1;
   
   while (i + j <= n)&(L(i+j) == L(i))
      
      ctr = ctr + (i + j);
      
      j = j + 1;
      
   end
   
   if j > 1
      
      R(i:i+j-1) = ctr/j;  %Revise ranks
      
      D2 = D2 - j*(j-1)/2; %Revise D2
      
      i = i + j;
      
   else
      
      i = i + 1;
      
   end
   
end

R2(rows) = R;


%*******************************

%Now compute Kendall's Tau

%*******************************

S = 0;

for i = 1:n-1
   
   x1 = R1(i);
   
   y1 = R2(i);
   
   for j = i+1:n
      
      x2 = R1(j);
      
      y2 = R2(j);
      
      S = S + sign(x1 - x2)*sign(y1 - y2);
      
   end
   
end

if D1*D2 == 0
   
   D1 = D1;
   
end


tau = S/sqrt(D1*D2);

%Compute corresponding z value

if S > 0 %Continuity correction
   
   S = S - 1;
   
elseif S < 0
   
   S = S + 1;
   
end

var = (1/18)*n*(n-1)*(2*n + 5);  %Kendall, p.50

sig = sqrt(var);

z = S/sig;