*! 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