Couldn't find it anywhere, so for future users who stumble on this thread, 
here is some code.
Note: If your data has missing values, delete those observations from the 
data first before running this code.

Place these functions in the global environment (ie run the code below).

Then, to obtain the test statistic: H_bisym_stat(x,y)  (x and y are the 
two vectors you are comparing).

To obtain a p-value: H_bisym_stat(x,y). The p-value is calculated via a 
Monte Carlo algorithm.

You can use the permutations argument to H_bisym_stat to increase the 
number of Monte Carlo samples.


#Hollander's test of bivariate symmetry

        #Calculating test statistic for Hollander's test of bivariate 
symmetry

H_bisym_stat <- function(x,y){

                        minXY <- pmin(x,y);
                        sortX <- x[order(minXY)];
                        sortY <- y[order(minXY)];
                        minXY <- pmin(sortX, sortY);
                        R <- as.numeric(sortX <= sortY);
                        maxXY <- pmax(sortX,sortY);
 
                        D1 <- matrix(0, length(x), length(x));
                        D2 <- matrix(0, length(x), length(x));

                        for (i in 1:length(x)){

                                for (j in 1:length(x)){

                                        D1[i,j] = as.numeric(minXY[j] < 
maxXY[i] & maxXY[i] <= maxXY[j]);
                                        D2[i,j] = as.numeric(minXY[i] <= 
minXY[j]);

                                        }

                        }

                        D = D1*D2;

                        S = 2*R - 1;

                        T = t( t(S) %*% D);

                        H_squared = (1/length(x)^2)*( t(T) %*% T );

                        H_squared};

        #Create a permutation sample for paired data;

permsamp_paired <- function(x,y){

                        new <- rbinom(length(x),1,0.5);
                        x_new <- x*new + y*(1 - new);
                        y_new <- y*new + x*(1 - new);

                        list(x_new = x_new, y_new = y_new)};


        #Calculate pvalue for Hollander statistic;

H_bisym_pv <- function(x, y, permutations = 1000){

                        H_squared = H_bisym_stat(x=x, y=y);

                        numb_greater = 0;

                        for (i in 1:permutations) {

                                newvars <- permsamp_paired(x=x, y=y);
                                x_new <- newvars$x_new;
                                y_new <- newvars$y_new;

                                newstat <- H_bisym_stat(x=x_new, y=y_new);

                                numb_greater = numb_greater + 
as.numeric(newstat >= H_squared);
 
                                }

                        pv <- numb_greater/permutations;

                        pv };


Joe Boyer
Statistical Sciences 
GlaxoSmithKline
	[[alternative HTML version deleted]]


