*! version 4.1.0  26sep97 TJS   STB-41 sbe19  STB-44 sbe19.1
*  modification of ktau to allow N==2, un-continuity-corrected
*  z and p values, and to pass more parameters
program define ktau2
	version 4.0
        local varlist "req ex min(2) max(2)"
        local if "opt"
        local in "opt"
        parse "`*'"
	parse "`varlist'", parse(" ")
        local x "`1'"
        local y "`2'"
	tempname k N NN pval score se tau_a tau_b
	tempname xt xt2 xt3 yt yt2 yt3
	tempvar doit nobs order work
	mark `doit' `in' `if'
	markout `doit' `x' `y'
	quietly count if `doit'
	scalar `N' = _result(1)
if `N' < 2 { error 2001 } 
	local Nmac = `N'
	quietly {
		gen long `order' = _n  /* restore ordering at end */
		replace `doit' = -`doit'
		sort `doit'  /* put obs for computation first */
		gen double `work' = 0  /* using type double is fastest */
                scalar `k' = 2
                while (`k' <= `N') {
			local kk = `k' - 1
			#delimit ;
			replace `work' = `work'
				+ sign((`x' - `x'[`k'])*(`y' - `y'[`k']))
				in 1/`kk' ;  /* using "in" is fastest */
			#delimit cr
                       	scalar `k' = `k' + 1
                }
		replace `work' = sum(`work') in 1/`Nmac'
		scalar `score' = `work'[`N']
	/* Calculate ties on `x' */
		egen long `nobs' = count(`x') in 1/`Nmac', by(`x')
tempvar nobsxm
egen `nobsxm' = max(`nobs')
	/* Calculate correction term for ties on `x' */
		replace `work' = sum((`nobs' - 1)*(2*`nobs' + 5)) in 1/`Nmac'
		scalar `xt' = `work'[`N']
	/* Calculate correction term for pairs of ties on `x' */
		replace `work' = sum(`nobs' - 1) in 1/`Nmac'
		scalar `xt2' = `work'[`N']
	/* Calculate correction term for triplets of ties on `x' */
		replace `work' = sum((`nobs' - 1)*(`nobs' - 2)) in 1/`Nmac'
		scalar `xt3' = `work'[`N']
	/* Calculate ties on `y' */
		drop `nobs' 
		egen long `nobs' = count(`y') in 1/`Nmac', by(`y')
tempvar nobsym
egen `nobsym' = max(`nobs')
	/* Calculate correction term for ties on `y' */
		replace `work' = sum((`nobs' - 1)*(2*`nobs' + 5)) in 1/`Nmac'
		scalar `yt' = `work'[`N']
	/* Calculate correction term for pairs of ties on `y' */
		replace `work' = sum(`nobs' - 1) in 1/`Nmac'
		scalar `yt2' = `work'[`N']
	/* Calculate correction term for triplets of ties on `y' */
		replace `work' = sum((`nobs' - 1)*(`nobs' - 2)) in 1/`Nmac'
		scalar `yt3' = `work'[`N']
	/* Compute Kendall's tau-a, tau-b, s.e. of score, and pval */
		scalar `NN'    = `N'*(`N' - 1)
		scalar `tau_a' = 2*`score'/`NN'
		scalar `tau_b' = 2*`score'/sqrt((`NN' - `xt2')*(`NN' - `yt2'))
		#delimit ;
scalar `se' = `NN'*(2*`N' + 5);
tempname tmax;
scalar `tmax' = max(`nobsxm', `nobsym');
if `tmax' > 1 { scalar `se' = `se'
  - (`xt' - `yt')
  + `xt3'*`yt3'/(9*`NN'*(`N' - 2))
  + `xt2'*`yt2'/(2*`NN') } ;
scalar `se' = sqrt((1/18)*`se');
		#delimit cr
		local zcc = (abs(`score') - 1) / `se'
local z = `score '/ `se'
tempname pvalcc
		if `score' == 0 {
             scalar `pval' = 1
             scalar `pvalcc' = 1 
          }
		else scalar `pvalcc' = 2*(1 - normprob((abs(`score') - 1)/`se'))
		else scalar `pval' = 2*(1 - normprob(abs(`score')/`se'))
	/* Restore original ordering of data set */
		sort `order'
	}
/* Print results */
	#delimit ;
       	di _n
		in gr "  Number of obs = " in ye  %7.0f `N' _n
        	in gr "Kendall's tau-a = " in ye %12.4f `tau_a' _n
        	in gr "Kendall's tau-b = " in ye %12.4f `tau_b' _n
        	in gr "Kendall's score = " in ye  %7.0f `score' _n 
        	in gr "    SE of score = " in ye %11.3f `se' _c ;
	if `xt2' > 0 | `yt2' > 0 { di in gr "   (corrected for ties)" _c } ;
       	di _n(2)
		in gr "Test of Ho: `x' and `y' independent" _n
in gr "             z  = " in ye %12.2f `z' _n
       		in gr "       Pr > |z| = " in ye %12.4f = `pval' _n(2)
in gr "             z  = " in ye %12.2f sign(`score')*`zcc' _n
       		in gr "       Pr > |z| = " in ye %12.4f = `pvalcc'
		in gr "  (continuity corrected)" ;
	#delimit cr
local c = 0
if `xt2' > 0 | `yt2' > 0 { local c = 1 }
	global S_1 = `N'
	global S_2 = `tau_a'
	global S_3 = `tau_b'
	global S_4 = `score'
	global S_5 = `se'
	global S_6 = `pval'
global S_7 = `z'
global S_8 = `pvalcc'
global S_9 = `zcc'
global S_10 = `c'
end