*! 1.0.0 NJC 27 Sept 2007
program bkrosenblatt, rclass
	version 9.2
	syntax varlist(min=2 max=2 numeric) [if] [in]
	marksample touse
	qui count if `touse'
	if r(N) == 0 error 2000

	tokenize `varlist'
	mata : bkr("`1'", "`2'", "`touse'")

	di _n as txt "n       = " as res r(n)            ///
	   _n as txt "n B_n   = " as res %5.4f r(n_B_n)  ///
	   _n as txt "z-score = " as res %5.4f r(z)      ///
           _n as txt "P-value = " as res %5.4f r(P_n_B_n)

	if "`crit'" != "" {
		tokenize `crit'
		di as txt "N.B. critical values at P = 0.05: " as res "`1'" ///
                _n as txt "                        P = 0.01: " as res "`2'"
	}

	return scalar P_n_B_n = r(P_n_B_n)
	return scalar z       = r(z)
	return scalar n_B_n   = r(n_B_n)
	return scalar n       = r(n)
end

mata :

void bkr(string scalar xname, string scalar yname,
        |string scalar tousename)
{
	real matrix x
	scalar n, i, n_B_n, b
	real colvector N1, N2, N3, N4

	if (args() == 3) x = st_data(., (xname, yname), tousename)
	else x = st_data(., (xname, yname))
	x = select(x, rowmissing(x) :== 0)
	n = rows(x)
	N1 = N2 = N3 = N4 = J(n, 1, .)

	for (i = 1 ; i <= n; i++) {
		N1[i] = sum(x[,1] :<= x[i,1] :& x[,2] :<= x[i,2])
		N2[i] = sum(x[,1] :>  x[i,1] :& x[,2] :<= x[i,2])
		N3[i] = sum(x[,1] :<= x[i,1] :& x[,2] :>  x[i,2])
		N4[i] = sum(x[,1] :>  x[i,1] :& x[,2] :>  x[i,2])
	}

	n_B_n = sum((N1 :* N4 :- N2 :* N3):^2) / n^4
	st_numscalar("r(n)", n)
	st_numscalar("r(n_B_n)", n_B_n)

	// note that as the transformation power h_n is negative,
	// the transformation reverses high and low values; hence
	// we negate the z-score from bkr_z()
	z = -bkr_z(n_B_n, n)
	st_numscalar("r(z)", z)
	st_numscalar("r(P_n_B_n)", normal(-z))

	if (n < 15) {
		if (n == 5)       st_local("crit", "0.0976 0.1408")
		else if (n == 6)  st_local("crit", "0.0872 0.1435")
		else if (n == 7)  st_local("crit", "0.0875 0.1312")
		else if (n == 8)  st_local("crit", "0.0842 0.1274")
		else if (n == 9)  st_local("crit", "0.0820 0.1250")
		else if (n == 10) st_local("crit", "0.0802 0.1218")
		else if (n == 11) st_local("crit", "0.0783 0.1195")
		else if (n == 12) st_local("crit", "0.0771 0.1172")
		else if (n == 13) st_local("crit", "0.0762 0.1163")
		else if (n == 14) st_local("crit", "0.0750 0.1137")
	}
}

real scalar bkr_z(real scalar n_B_n, real scalar n)
{
	real scalar mu, sigma, h_n, z

	if (n < 15) return(.)
	else {
		mu = n < 25 ? 4.663 - 1 / (0.2137 + 0.00448 * n) :
			3.823 - 1 / (0.193 + 0.01662 * n^0.8481)
		sigma = 0.614 - 1 / (1.187 + 0.0328 * n)
		h_n = -0.36 + 2.866 * n^(-0.775) - 0.683 * exp(-0.244 * n)
		z = (n_B_n^h_n - mu) / sigma
		return(z)
	}
}

end