function(x=c(8 , 26 , -7 , -1 , 2 , 9 , 0 , -4 , 13 , 3 , 3 , 4),
  y=c(14 , 18 , 2 , 4 , -5 , 14 , -3 , -1 , 1 , 6 , 3 , 3),Nsim=10000){
N=length(x)
dyx=y-x ; adyx=abs(dyx)
signyx=rep(0,length(dyx)) ; signyx[dyx>0]=1;signyx[dyx<0]=-1
rd=rank(adyx)*signyx
Vstar=sum(rd[rd>0])
rdr=abs(rd[signyx!=0]) ; Nr=length(rdr)
Vvec=NULL
for(i in 1:Nsim){
signyxr=rbinom(Nr,1,.5)*2-1
Vvec[i]=sum(rdr[signyxr>0])
}
pval.sim=mean(Vvec>=Vstar)
pval.sim
}