/* procedure to calculate Kendall's tau 
input: 
x 	n-vector or real numbers
y	n-vector of real numbers
return
t	Kendall's tau, or -9 in case of an error
*/

proc k_tau(x,y);
 local n_c, n_d, i, n, z;
 n=rows(x);
 if (rows(y)/=n);
  errorlog "k_tau.g";
  errorlog "x and y should be of equal size";
  retp(-9);
 endif;
 z=sortc(x~y,1);
 i=1;
 n_c=0;
 n_d=0;
 do while (i<=n);
  /* concordant pairs */
  n_c=n_c+sumc(z[i:n,2].>z[i,2]);
  /* discordant pairs */
  n_d=n_d+sumc(z[i:n,2].<z[i,2]);
  i=i+1;
 endo;
 retp( 2*(n_c-n_d)/(n*(n-1)) );
endp;