An asymptotically optimal transform of Pearson's correlation statistic
Iosif Pinelis
math.ST stat.TH 62E20, 62F03, 62F12





Optimal 
transform of Pearson's  
statistic

















Iosif Pinelis



Department of Mathematical Sciences

Michigan Technological University

Houghton, Michigan 49931, USA

E-mail: [ipinelis@mtu.edu]e1



It is shown that for any correlation-parametrized model of dependence and any given significance level ∈(0,1), 
there is an asymptotically optimal transform of Pearson's correlation statistic R, for which the generally leading error term for the normal approximation vanishes 
for all values ρ∈(-1,1) of the correlation coefficient. 

This general result is then applied to the bivariate normal (BVN) model of dependence and to what is referred to in this paper as the SquareV model.  
In the BVN model, Pearson's R turns out to be asymptotically optimal for a rather unusual significance level ≈0.240, whereas Fisher's transform R_F of R is asymptotically optimal for the limit significance level =0. 
In the SquareV model, Pearson's R is asymptotically optimal for a still rather high significance level ≈0.159, whereas Fisher's transform R_F of R is not asymptotically optimal for any ∈[0,1].
Moreover, it is shown that in both the BVN model and the SquareV model, 
the transform optimal for a given value of  is in fact asymptotically better than R and R_F in wide ranges of values of the significance level, including  itself. 

Extensive computer simulations for the BVN and SquareV models of dependence are presented, which suggest that, for sample sizes n≥100 and significance levels ∈{0.01,0.05}, the mentioned asymptotically optimal transform of R generally outperforms both Pearson's R and Fisher's transform R_F of R, the latter appearing generally much inferior to both R and the asymptotically optimal transform of R in the SquareV model. 







  



[class=
MSC2010]

62E20
62F03
62F12
















hypothesis testing
Pearson's 
correlation statistic
Fisher's z
transform
asymptotically optimal transform
models of dependence
copulas
bivariate normal distribution




chapter



subsubsection








plain 















equationsection








§ 
INTRODUCTION
 


A statistic closely related to Pearson's R is commonly known as the Fisher z transform, defined by the formula 

    R_F:=tanh^-1(R)=1/2ln1+R/1-R.

An advantage of using R_F (as opposed to R) in making statistical inferences about the true correlation coefficient ρ is usually ascribed to its variance-
stabilizing property in normal populations; see e.g. Fisher <cit.>, Gayen <cit.>, and Hotelling <cit.>; that is, n R_F→1 for all ρ∈(-1,1) as n→∞ (as opposed to n R →(1-ρ^2)^2) whenever 
the underlying distribution is bivariate normal. 
Everywhere here, n denotes the sample size. 


In his discussion of Hotelling's paper <cit.>, Kendall


provides heuristics suggesting that such variance stabilization of the distribution of a statistic may often result in it being closer to normality. 
Namely, if an approximate constancy of the variance of a statistic were the same as an approximate constancy of its distribution itself, and if the distribution is close to normality at least for one value of the parameter (say, ρ, as in the present case), then it would be close to normality for all values of ρ. 

For normal populations and large enough sample sizes, the Fisher z transform indeed brings the distribution of the correlation statistic closer to normality, 

and it is especially effective for values of ρ far from 0. 
However, it is well known 
(see e.g. <cit.>) 
that the closeness of the distribution of a statistic to normality is usually mainly determined, not by the variance, but by the third moments of the appropriately standardized statistic. 






In this paper, we shall see that for a general and most common class of models of dependence, including the bivariate normal (BVN) model, and
for each given significance level ∈(0,1) there is a certain transform Ψ_(R) of Pearson's statistic R that assures the vanishing of the generally leading term of the asymptotics of the probability that an approximately standardized version of the statistic Ψ_(R) exceeds the standard normal critical value 

    z_:=Φ^-1(1-);

here, as usual, Φ is the standard normal cumulative distribution function (cdf) and Φ^-1 is its inverse; unless otherwise specified, all the asymptotics here are for large sample sizes n.  
Thus, the 
transform Ψ_(R) of R is asymptotically optimal: its distribution is asymptotically the closest to normality exactly at the critical value. 

Once this optimality result is obtained for the general class of models of dependence, the rest of the paper is devoted to detailed analysis of the optimal transform Ψ_(R) of R in the BVN model and another specific model of dependence, referred to in this paper as the SquareV model. 

We shall 
see that, in the BVN model, the mentioned family (Ψ_(R))_∈(0,1) of transforms of Pearson's statistic R includes R itself: namely, 

    R=Ψ__P(R),

where 

    _P:=1-Φ(1/√(2))≈0.240.

Thus, Pearson's statistic R is asymptotically optimal for a significance level  of about 24%, but such a significance level is rather unusual in statistical practice. 

As for Fisher's transform R_F of R, we shall see that, again in the BVN model,  

    R_F=Ψ_0(R):=lim_↓0Ψ_(R),

which means that Fisher's transform is asymptotically optimal for the significance level 

    _F:=0.


Now one might explain the fact that for the usually rather small significance levels, such as 0.05 or 0.01, Fisher's statistic R_F is asymptotically closer to normality than Pearson's statistic R by noting that the significance level _F=0 (for which R_F is asymptotically optimal) is closer to 0.05 and especially to 0.01 than the significance level _R≈0.240 (for which R is asymptotically optimal). 

As for the SquareV model, there Pearson's R is asymptotically optimal for a still rather high significance level ≈0.159, whereas Fisher's transform R_F of R is not asymptotically optimal for any ∈[0,1].

It should be noted that in both the BVN model and the SquareV model, 
the transform optimal for a given value of  is in fact asymptotically better than R and R_F in wide ranges of values of the significance level, including  itself. 

We have also conducted extensive computer simulations for the BVN and SquareV models of dependence, which suggest that, for sample sizes n≥100 and significance levels ∈{0.01,0.05}, the mentioned asymptotically optimal transform of R generally outperforms both Pearson's R and Fisher's transform R_F of R, the latter appearing generally much inferior to both R and the asymptotically optimal transform of R in the SquareV model.

The rest of the paper is organized as follows. 

In Section <ref> we present an asymptotic expansion for statistics that are general smooth nonlinear functions of the sample mean of iid random vectors in ^p. This expansion, which may be viewed as a far-reaching refinement of the delta method, is a special case of results by Bhattacharya and Ghosh <cit.>. For Berry–Esseen-type bounds for general nonlinear statistics, see e.g. <cit.>.  

In Section <ref>, the mentioned asymptotic expansion is specialized for the cases of Pearson's correlation statistic R and its smooth enough transforms. A key observation there is that the main term of the asymptotic for such a transform of R differs from the corresponding main term for R itself only by a comparatively simple expression involving the first two derivatives of the transform function ψ. This allows one to obtain, for any correlation-parametrized model of dependence and for any given significance level ∈(0,1), a rather simple second-order ordinary differential equation (ODE) for the optimal transform function ψ that makes the main term of the asymptotic for the asymptotically optimal transform of R vanish for all values ρ∈(-1,1) of the correlation coefficient. This ODE can be explicitly solved for a number of models of dependence, including the important BVN model and models with a linear dependence of the correlation parameter. The mentioned SquareV model is a model with such a linear dependence. 

The BVN model is considered in detail in Section <ref>, where the corresponding results of the mentioned computer simulations are also presented and discussed. 

A similar treatment of the SquareV model is given in Section <ref>. 

Section <ref> is a summary of the results of this paper. 





§ ASYMPTOTIC EXPANSIONS FOR SMOOTH NONLINEAR STATISTICS
 

Let 

    V,V_1,V_2,…

be independent identically distributed (iid) zero-mean random vectors in ^p with 
V^3<∞, where p is a natural number and · is the Euclidean norm in ^p, which latter will be identified, as usual, with the space of all p×1 column matrices. 
Assume also that the Cramér-type condition lim sup_t→∞|exp(it^T V)|
<1 is satisfied, where i is the imaginary unit and ^T denotes the transposition, in this case of a column matrix t∈^p; for this Cramér-type condition to hold, it is enough  
that, for some natural k, the k-fold convolution of the distribution of V have a nonzero absolutely continuous component. 

Let  stand for the covariance matrix of V: 

    := VV^T.


Let f^p→ be a function which is twice continuously differentiable in a neighborhood of 0∈^p and such that f(0)=0. Let L and H denote, respectively, the gradient vector and the Hessian matrix of the function f at 0, so that 

    f'(0)(v)=L^T v  and   f”(0)(v,v)=v^T Hv

for all v∈^p. 
Since V is assumed to be zero-mean, one has L^T V=0. 
Introduce now 

    :=√((L^T V)^2),

which will be assumed to be nonzero, so that 

    :=L^T V/

is a well-defined r.v., with zero mean and unit variance. 
Consider the r.v.   

    T_n:=√(n)/ f(V̅),

where of course V̅:=1/n ∑_1^n V_i. 




Then, by Theorem 2 of the paper <cit.> by Bhattacharya and Ghosh, 

    sup_z∈|(T_n≤ z)-Ψ_3,n(z)|=o(1/√(n)),

where 





    Ψ_3,n(z)   :=Φ(z)+(z)/√(n), 
    (z)   :=-[(^3/6+a_3) (z^2-1)+a_1] (z) 
       =(Az^2+B) (z),


Φ and  denote, as usual, the distribution and density functions of N(0,1), and a_1, a_2, A, and B are constants depending only on L, H, , , and ^3 (but not on z or n): 

    a_1   :=1/2  H, 
     
    		a_3   :=1/4^3 (L^T L-^2)(H)+1/2^3 L^T H L, 
    
    	A   :=-(^3/6+a_3), 
    
    	B   :=-A-a_1,

with  denoting the trace of a matrix. 


If the condition V^3<∞ is strengthened to V^4<∞, then o(1/√(n)) in (<ref>) can be replaced by O(1/n). 




§ ASYMPTOTICS FOR THE PEARSON STATISTIC AND ITS TRANSFORMS
 
Let (Y,Z),(Y_1,Z_1),…,(Y_n,Z_n) be independent identically distributed random points in ^2 with a correlation coefficient ρ∈(-1,1) and (Y^6+Z^6)<∞. Pearson's sample correlation coefficient based on the observations (Y_1,Z_1),…,(Y_n,Z_n) is defined by the formula 

    R:=R_n:=Y̅Z̅-Y̅ Z̅/√(Y̅^̅2̅-Y̅^2)√(Z̅^̅2̅-Z̅^2),

where 
Y̅:=1/n∑_1^n Y_i, Z̅:=1/n∑_1^n Z_i, Y̅^̅2̅:=1/n∑_1^n Y_i^2, Z̅^̅2̅:=1/n∑_1^n Z_i^2, and Y̅Z̅:=1/n∑_1^n Y_iZ_i; 








let 
R take an arbitrarily assigned value in the interval [-1,1] if the denominator of the ratio in (<ref>) is 0. 





Let us assume that Y and Z are each standardized, that is, zero-mean and unit-variance. This assumption does not diminish generality, because R is invariant with respect to affine transformations Y_i↦ a+b Y_i and Z_i↦ c+d Z_i of the Y_i's and Z_i's, for any real constants a,b,c,d such that b>0 and d>0. 

Observe that 

    R-ρ=f(V̅),
 
where 

    V:=(Y,Z,Y^2-1,Z^2-1,YZ-ρ)

and 

    f(v):=f_ρ(v):=ρ+v_5-v_1v_2/√(1+v_3-v_1^2)√(1+v_4-v_4^2)-ρ

if v=(v_1,…,v_5)∈^5 is such that 1+v_3-v_1^2>0 and 1+v_4-v_4^2>0; 
otherwise, let f(v):=
0. 
In this case, L=(0,0,-ρ/2,-ρ/2,1), whence 

    =√((YZ-ρ2 (Y^2+Z^2))^2)  and  =YZ-ρ2 (Y^2+Z^2)/.



As noted in <cit.>, the condition =0 is equivalent to the following exceptional situation: there exists some ∈ such that the random point (Y,Z) lies almost surely on the union of the two straight lines through the origin with slopes  and 1/ (for =0, these two lines should be understood as the two coordinate axes in the plane ^2). 


It will be assumed in what follows that the random point (Y,Z) is such that  is never 0. 
 
Then it is easy to check that all the conditions on f and V stated in Section <ref> are satisfied, with p=5; 



in particular, f(0)=0. 


Let now _R(z) denote (z) defined by (<ref>) with f as in (<ref>). 

Further, letting 

    μ_ij:= Y^iZ^j,


one has 

    =1/2 √(ρ^2 (μ _04+2 μ _22+μ _40)-4 ρ(μ _13+μ
       _31)+4 μ _22),


    ^3^3=-ρ^3/8(μ _06+3 μ _24+3 μ _42+μ_60)
    	+6 ρ^2 (μ _15+2 μ _33+μ _51) 
     
       -12 ρ(μ _24+μ _42)+8 μ _33,

and 

    96^3/(z)_R(z) = _R(z)
     
    :=16 [(z^2-1) (6 μ_12μ_21-μ_33)+3 ^2 z^2 (μ_13+μ_31)]  
     
    -12 ρ[(z^2-1) (4 μ_03μ_21+4
       μ_12μ_30+8 μ_12^2-2 μ_13μ_31+μ_13^2+8 μ_21^2-2 μ_24+μ_31^2-2 						μ_42) 
        +^2 ((2z^2+1) (μ_04+μ_40)+(4 z^2-2) μ_22
       )] 
     
    +12 ρ^2 (z^2-1) [2 μ_03 (3 μ_12+μ_30)+μ_04 (μ_13-μ_31)+10 μ_12μ_21-μ_13μ_40-μ_15
        +6 μ_21μ_30+μ_31μ_40-2 μ_33-μ_51] 
     	
    -ρ^3 (z^2-1) [24 μ _03μ _21+12 μ _03^2-6 
       μ_04μ_40+3 μ_04^2-2 μ_06
        +24 μ_12μ_30+12 μ_12^2+12 μ_21^2-6 μ_24+12 μ_30^2+3 μ_40^2-6 μ_42-2 			μ_60].


More generally, let now 

    ψ(-1,1)→

be a twice continuously differentiable function whose derivative ψ' does not vanish at any point of the interval (-1,1). Let then 

    g(v):=g_ρ(v):=ψ(f(v)+ρ)-ψ(ρ)/ψ'(ρ)

for v∈^5, with f as defined in (<ref>). 
Then, in view of (<ref>), 

    ψ(R)-ψ(ρ)/ψ'(ρ)=g(V̅),

Note that 
g(0)=f(0)=0 and g'(0)=f'(0), so that  and  for the function g are the same as in (<ref>) (given there for the function f). 
Thus, all the conditions stated in Section <ref> are satisfied with g in place of f. 

Let now _ψ(R)(z) denote (z) defined by (<ref>) with g as in (<ref>) in place of f. 
Then (<ref>) will hold with 

    τ_ψ,n:=ψ(R)-ψ(ρ)/ψ'(ρ)/√(n)

in place of T_n and _ψ(R)(z) in place of (z).  
One may note here that τ_ψ,n may be considered an asymptotically standardized version of ψ(R). 


A key observation is that 

    _ψ(R)(z)=_R(z)-ψ”(ρ)/2ψ'(ρ)  z^2 (z).


To begin using this observation, let us refer to any family (P_ρ)_ρ∈(-1,1) of distributions of the random pair (Y,Z) in ^2 parametrized by the correlation coefficient ρ of (Y,Z) as a correlation-parametrized model (CP) of dependence. 
The CP condition seems quite natural for parametric models of dependence. Indeed, let a real parameter þ represent the strength of the dependence between Y and Z. Then one should usually expect the correlation coefficient ρ to be a strictly increasing continuous function g of þ: ρ=g(þ). Replacing then þ by g^-1(ρ), one obtains a re-parametrization with ρ as the new parameter. 

In this regard, one may recall the formula 

    (Y,Z)=∬_^2a_Y,Z(y,z) dy dz

for the covariance of Y and Z, where a_Y,Z is the association function of r.v.'s Y and Z, given by the formula 

    a_Y,Z(y,z):=(Y>y,Z>z)-(Y>y)(Z>z)

for all (y,z)∈^2. 
Since (Y>y)=(-Y<-y), formula (<ref>) can be rewritten as 

    (Y,Z)=∬_^2[F_Y,Z(y,z)-F_Y(y) F_Z(z)] dy dz,

where F_Y,Z is the joint cdf of the random pair (Y,Z), and F_Y and F_Z are the corresponding marginal cdf's. 

Suppose now that we have a dependence model (_þ), where þ is a strength of the association/dependence parameter, so that the association function a_þ;Y,Z of the pair (Y,Z) with respect to the probability measure _þ is increasing in þ on the average in the sense that 
the integral in (<ref>) with a_þ;Y,Z in place of a_Y,Z is increasing in þ. Suppose also that the Y- and Z-marginals of the distribution of the pair (Y,Z) with respect to _þ do not depend on þ. Then the correlation coefficient ρ of (Y,Z) will be an increasing function of þ. 
 
A more specific, but still rather general way to construct a CP model is as follows. By Sklar's theorem (see e.g. <cit.>), 


    F_Y,Z(y,z)=C(F_Y(y),F_Z(z))

for some copula C and all (y,z)∈^2; recall that a copula can be defined as the joint cdf of a random pair with values in the unit square [0,1]^2 whose marginals are uniform on the interval [0,1]. 
Let now (C_þ) be any family of copulas increasing in þ; a large number of such families can be found in <cit.>. Fix the marginal cdf's F_Y and F_Z, and for each value of the parameter þ let F_þ;Y,Z(y,z):=C_þ(F_Y(y),F_Z(z)), again for all (y,z)∈^2. Then the correlation coefficient ρ corresponding to the joint cdf F_þ;Y,Z will be an increasing function of þ. 

In view of (<ref>), in any correlation-parametrized model of dependence and for any given real z0,  


the expressions in (<ref>) and (<ref>) for  
 and _R(z) 
will depend on ρ only. 
Then, by the key observation (<ref>), the condition _ψ(R)(z)=0 can be rewritten as the second-order ordinary differential equation (ODE)

    ψ”(ρ)/ψ'(ρ)=h_z(ρ)

for the function ψ, where 

    h_z(ρ):=_R(z)/48^4 z^2.


Solving now ODE (<ref>) with the natural initial conditions

    ψ(0)=0  and  ψ'(0)=1,

we have 

    ψ'(ρ)=exp∫_0^ρ dr  h_z(r)

and 

    ψ(ρ)=ψ_z(ρ):=∫_0^ρ dr exp∫_0^r ds  h_z(s)

for ρ∈(-1,1); in (<ref>) and (<ref>), 

we use the common convention 
∫_0^s:=-∫_s^0 for s<0. 
Thus, we obtain 


In any correlation-parametrized model of dependence and for any given real z0, 
the generally leading error term for the normal approximation for ψ_z(R) vanishes:  

    _ψ_z(R)(z)=0

for all ρ∈(-1,1). 


Letting now 

    Ψ_:=ψ_z_,

we can rewrite (<ref>) as 

    _Ψ_(R)(z_)=0

for all ∈(0,1), with z_=Φ^-1(1-), as defined in (<ref>). 

One may note here that for a rather large class of models of dependence the functions h_z will be rational, and hence, according to (<ref>), ψ'_z will be an elementary, closed-form function. This class of models with rational functions h_z includes 
the bivariate normal model and models with linear dependence of the joint cdf F_þ;Y,Z on þ. In particular, the class of models with linear dependence of F_þ;Y,Z on þ contains Farlie's model <cit.>. 

We shall consider the bivariate normal model and a particular simple model with linear dependence of F_þ;Y,Z on þ in the following sections, to compare the performance of Pearson's R itself, its Fisher transform R_F, and the asymptotically optimal transform Ψ_(R) of R in non-asymptotic settings, for specific sample sizes. 






§ BIVARIATE NORMAL MODEL (BVN)




 §.§ Asymptotically optimal transform  in the BVN model


Here it is assumed that the random point (Y,Z) has the bivariate normal (BVN) distribution with zero means, unit variances, and an arbitrary correlation coefficient ρ∈(-1,1). Then the expressions for _R and , and thus for h_z(ρ), in formulas (<ref>), (<ref>), and (<ref>) can be greatly simplified. 

Indeed, in this case the pair (Y,Z) equals (Y,ρ  Y+√(1-ρ^2) Y_1) in distribution, whence, by (<ref>),   

    μ_ij=∑ _k=0^j jkρ^k(1-ρ^2)^(j-k)/2 m(i+k) m(j-k)

for all i,j=0,1,… and 

    m_j:= Y^j,

so that (m_0,…,m_6)=(1,0,1, 0, 3, 0, 15). 
As the result, ODE (<ref>) becomes  

    ψ”(ρ)/ψ'(ρ)=
    	p_z -2ρ/1-ρ ^2,

where 

    p_z:=1/2z^2-1.

ODE (<ref>) is easily solved, yielding 

    ψ'(ρ)=(1 - ρ^2)^p_z

and 

    ψ(ρ)=ψ_z(ρ):=∫_0^ρ(1 - r^2)^p_z dr
    	=ρ  _2F_1(12,-p_z;32;ρ ^2),

where _2F_1 is the ordinary hypergeometric function, given by the formula 

    _2F_1(a,b;c;x)=∑_k=0^∞(a)_k(b)_k/(c)_k x^k/k!

for x with |x|<1, where (q)_k:=∏_i=0^k-1(q+i) is the Pochhammer symbol. 
The last equality in (<ref>) can be obtained by expanding the integrand (1 - r^2)^p_z into the Maclaurin series in powers of r and then integrating the series term-wise. 

Thus, recalling (<ref>), we see that in the bivariate normal case the transform Ψ_(R)=ψ_z_(R) of Pearson's R with ψ_z as in (<ref>) is asymptotically optimal for any given significance level ∈(0,1). 

In particular, choosing

    z=1/√(2)≈0.707,

we have p_z=0. Hence, in view of the integral expression in (<ref>), ψ_z(ρ)=ρ, so that we have   
(<ref>), confirming that the family (Ψ_(R))_∈(0,1) of transforms of Pearson's statistic R includes R itself. 

On the other hand, letting z→∞, we have p_z→-1, so that, using again the integral expression in (<ref>) (and, say, the dominated convergence theorem), we see 

    ψ_z(ρ)z→∞⟶ψ_∞(ρ):=1/2ln1+ρ/1-ρ

for all ρ∈(-1,1), thus confirming (<ref>). 

In view of formula (<ref>), the calculation of values of the functions ψ_z or, equivalently, of the functions Ψ_ mainly reduces to the calculation of values of the hypergeometric function _2F_1. In general, this hypergeometric function is not elementary. However, there are a number of highly efficient ways to compute values of _2F_1. It takes only about 1.7×10^-5 sec 

on an average to compute a value of ψ_2(ρ) (on a standard computer), which may be compared with the corresponding execution time of about 0.45×10^-5 sec for Fisher's ψ_∞(ρ)=1/2ln1+ρ/1-ρ. 
Therefore and because usually in statistical practice the value of the transform Ψ_(R)=ψ_z_(R) of the statistic R needs to be computed only once, the use of the hypergeometric function _2F_1 should not cause any complications. 

Also, according to (<ref>), the derivarive ψ' of the function ψ=ψ_z is a simple elementary expression, which makes it easy to obtain various analytical properties of ψ_z. For instance, using the special l'Hospital-type rule for monotonicity (see e.g. <cit.>), we can immediately see that the ratio ψ_z_1(ρ)/ψ_z_2(ρ) is decreasing in ρ^2 for any real z_1 and z_2 such that 0<|z_1|<|z_2|. In particular, it follows that the ratio of  ψ_z(ρ) to Fisher's ψ_∞(ρ)=1/2ln1+ρ/1-ρ is decreasing in ρ^2 for any real z0. One may also note that the values  

    ψ_z(±1)=±√(π) Γ (p_z+1)/2 Γ(p_z+3/2)

at the endpoints of the interval [-1,1] 
are finite for all real z0, in contrast with Fisher's limit values ψ_∞(±(1-))=±∞; here p_z is as defined in (<ref>). 
It is also clear that ψ_z(ρ) is odd in ρ, for each z0. 
 




 §.§ The transform  in the BVN model is asymptotically better than  and  in wide ranges of values of the significance level, including  itself
 


According to Theorem <ref>, for any given correlation-parametrized model of dependence and any given significance level ∈(0,1), the transform Ψ_(R) of R is asymptotically optimal for all ρ∈(-1,1) as the sample size n goes to ∞. In particular, for any given significance level ∈(0,1), the transform Ψ_(R) of R is asymptotically better than both R itself and its Fisher transform R_F. In fact, Ψ_(R) is asymptotically better than R and R_F for rather wide ranges of values (say ) of the significance level; of course, these ranges include the value  itself. 
Indeed, one can see that in the BVN model 

    _ψ(R)(z)=ρ/2 (2 z^2-1)-z^2/2(1-ρ ^2)ψ”(ρ )/ψ'(ρ).

In particular, 

    _R(z)=ρ/2 (2 z^2-1),
 

    _R_F(z)=_ψ_∞(R)(z)=-ρ/2,

where ψ_∞ is as defined in (<ref>), and 

    _Ψ_(R)(z)=_ψ_z_(R)(z)=ρ/2(z^2/z_^2-1).
 
So, 
|_Ψ_(R)(z_)|<|_R_F(z_)| for ρ0 if 0<z_<z_√(2). That is,
the transform Ψ_(R), which is asymptotically optimal for the given significance level ∈(0,1), will still be asymptotically better than Fisher's transform R_F for any significance level ∈(0,1) such that 0<z_<z_√(2). For instance, if =0.05, then Ψ_(R) will be asymptotically better than R_F, not just for the significance level =0.05, but for any significance level ∈(0.01000,0.5) – because 0<z_<z_0.05√(2) for all ∈(0.01000,0.5). Similarly, if =0.01, then Ψ_(R) will be asymptotically better than R_F for any significance level ∈(0.00050,0.5). 

As for the comparison of the asymptotically optimal transform Ψ_(R) with R itself, we can similarly see that, for instance, if =0.05, then Ψ_(R) will be asymptotically better than R, not just for the significance level =0.05, but for any significance level ∈(0,0.17912); if =0.01, then Ψ_(R) will be asymptotically better than R for any significance level ∈(0,0.16933). 



 §.§ Comparison of asymptotically optimal transform  in the BVN model with  and  for finite sample sizes
 


According again to Theorem <ref>, for any given correlation-parametrized model of dependence and any given significance level ∈(0,1), the transform Ψ_(R) of R is asymptotically optimal for all ρ∈(-1,1) as the sample size n goes to ∞. However, one may ask: how does this asymptotically optimal transform of R compare in performance with R itself and the Fisher transform of R for finite sample sizes n? 

To address this question for the BVN model, we have done the following. For each value ∈{0.01,0.05} of the significance level, each value ρ∈{0,0.1,0.5,0.9} of the correlation coefficient, and each value n∈{10,10^2,10^3,10^4} of the sample size, N:=10^6 samples of size n from the BVN distribution with zero means, unit variances, and correlation coefficient ρ∈(-1,1) were simulated using the command  of the computing system 
Mathematica, on each of the 12 working in parallel processing cores of a computer workstation. Thus, altogether 12×10^6×2×4×(10+10^2+10^3+10^4)≈10^12 realizations of the BVN random pair (Y,Z) were simulated, which took about 3.6 hours. This suggests that Mathematica's simulation of the BVN distribution is highly effective. (In view of the symmetry (Y,-Z)=-(Y,Z), negative values of ρ have not been considered here.)

For each quadruple 
    (,ρ,n,k)∈{0.01,0.05}×{0,0.1,0.5,0.9}×{10,10^2,10^3,10^4}×{1,…,12}
 
(where k indexes the 12 processing cores) and for each of the three transforms of R (R itself, its Fisher transform R_F, and the asymptotically optimal transform Ψ_(R) of R), the relative frequencies – say _ψ(R)(,ρ,n,k) – of the values of 

τ_ψ,n exceeding z_ were computed, where ψ(ρ)≡ψ_1/√(2)(ρ)≡ρ for R, ψ=ψ_∞ for R_F, ψ=ψ_z_ for Ψ_(R); and τ_ψ,n is as defined in (<ref>). So, each of these relative frequencies _ψ(R)(,ρ,n,k) is an estimate of the significance level , with the corresponding estimates 

    _ψ(R)(,ρ,n,k):=_ψ(R)(,ρ,n,k)/-1

of the relative errors 
of the approximation of  by _ρ(τ_ψ,n>z_), for each ψ∈{ψ_1/√(2),ψ_∞,
ψ_z_}. 
Then, for each triple (,ρ,n)∈{0.01,0.05}×{0,0.1,0.5,0.9}×{10,10^2,10^3,10^4} and each ψ∈{ψ_1/√(2),ψ_∞,ψ_z_}, the mean 


    _ψ(R)(,ρ,n):=1/12∑_k=1^12_ψ(R)(,ρ,n,k)

and standard deviation 

    s_ψ(R)(,ρ,n):=√(1/12-1∑_k=1^12(_ψ(R)(,ρ,n,k)-_ψ(R)(,ρ,n))^2, 
    )

of the 12 values 

    _ψ(R)(,ρ,n,1),…,_ψ(R)(,ρ,n,12)

were computed. Thus, again for each triple (,ρ,n) and each ψ∈{ψ_1/√(2),ψ_∞,ψ_z_}, _ψ(R)(,ρ,n) is the best available in this setting estimate of the relative error of the approximation of  by 
_ρ(τ_ψ,n>z_), and a reasonable estimate of the standard error of this estimate of the relative error of the approximation is 

    s̃_ψ(R)(,ρ,n):=s_ψ(R)(,ρ,n)/√(12).




The values of _ψ(R)(,ρ,n) and s̃_ψ(R)(,ρ,n) are presented in Table <ref>, along with the corresponding values of s̃_ψ(R)(,ρ,n), all of those values round to the precision of 3 decimal digits. 

For instance, the entry -0.236±0.00253 in the first row of Table <ref>, under the heading _R(,ρ,n)±s̃_R(,ρ,n), means that for the triple (,ρ,n)=(0.01 , 0 , 10) and ψ(R)=ψ_1/√(2)(R)=R, we have _ψ(R)(,ρ,n)=_R(,ρ,n)≈-0.236 and 
s̃_ψ(R)=s̃_R(,ρ,n)≈0.00253. For the same triple (,ρ,n)=(0.01 , 0 , 10), we have _R_F(,ρ,n)=_ψ_∞(R)(,ρ,n)≈1.63 and 
s̃_R_F(,ρ,n)=s̃_ψ_∞(R)(,ρ,n)≈0.00463, and also 
s̃_Ψ_(R)(,ρ,n)=s̃_ψ_z_(R)(,ρ,n)≈0.00439. 

For an easier grasp, graphs of the values of the means _ψ(R)(,ρ,n) are given in Fig. <ref>. The 8 pictures in Fig. <ref> correspond to the pairs (,ρ)∈{0.01,0.05}×{0,0.1,0.5,0.9} and are labeled accordingly. Each of these 8 pictures contains, for the given pair (,ρ) and for the three transforms – R, R_F, and R_*:=Ψ_(R) of R, the successively connected graphs of _ψ(R)(,ρ,n) as functions of n∈{10,10^2,10^3,10^4}. 

By Remark <ref>, it should be expected that the means _ψ(R)(,ρ,n) decrease in n roughly as 1/√(n). Therefore, to see the differences between the three graphs in each of the 8 pictures more clearly for the larger values of n, the corresponding graphs of the values of  _ψ(R)(,ρ,n) √(n) are given in Fig. <ref>. 


Looking at Fig. <ref> and back at Table <ref>, we see that for all considered triples (,ρ,n) the asymptotically optimal transform Ψ_(R) provides for a better normal approximation than Fisher's transform R_F does. The advantage of Ψ_(R) over R_F is smaller for the smaller values of ρ and , which should be understandable. 
Indeed, since all the three versions of the function ψ (namely, ψ_1/√(2), ψ_∞, and ψ_z_) satisfy the same initial conditions (<ref>) at ρ=0, all the three corresponding transforms of R (R itself, R_F, and Ψ_(R)) are close to one another if ρ is close to 0. Also, the smaller the value of  is, the closer is the critical value z_ to Fisher's value z_0+=∞. 

For the sample size n=10, the asymptotically optimal transform Ψ_(R) appears to have little or no advantage over R itself, which is also understandable, since the sample size n=10 is small, whereas the transform Ψ_(R) is guaranteed to be optimal only for large enough sample sizes; actually, for n=10 all the three transforms perform poorly or, mostly, very poorly – except for the rather strange case of (,ρ,n)=(0.05,0.1,10), when R performs surprisingly well. Also, the advantage of Ψ_(R) over R appears to be smaller for smaller values of ρ, a possible reason for which is that, as explained in the previous paragraph, all the three transforms of R are close to one another if ρ is close to 0; in fact, R appears to be somewhat better than Ψ_(R) for ρ=0 (except when =0.01, ρ=0, and n=10^4) and for =0.05, ρ=0.1, and n=10^2. However, for ρ∈{0.5,0.9} and n∈{10^2,10^3,10^4}, the transform Ψ_(R) appears to be much better than R itself. 

For the same values of ρ and n as in the last sentence, Fisher's transform R_F also appears to be much better than R, but even in these cases, the advantage of R_F over R is significantly less than that of Ψ_(R) over R. 



§ SQUAREV MODEL


Here we shall consider the dependence model, which is the family (P_ρ)_-1<ρ<1 of distributions of the random pair (Y,Z) on the vertices of the square [-1,1]×[-1,1] given by the following formulas: 

    _ρ((Y,Z)=(1,1))=_ρ((Y,Z)=(-1,-1))=1+ρ4, 
       _ρ((Y,Z)=(1,-1))=_ρ((Y,Z)=(-1,1))=1-ρ4.

In other words, the distribution P_Y,Z=P_ρ;Y,Z of (Y,Z) under _ρ is the mixture 

    {2   (1-ρ)P__1,_2+ρ P__1,_1       if ρ≥0, 
       (1+ρ)P__1,_2-ρ P__1,-_1       if ρ<0,  
    .

where _1,_2 are independent Rademacher r.v.'s, with (_j=±1)=1/2 for j=1,2. 
Then (Y,Z)=ρ under _ρ, so that the use of the symbol ρ to denote the parameter is consistent. 

The just described model of dependence will be referred to as the SquareV model, where “V” stands for “vertices”. 

In view of the previously mentioned symmetry (Y,-Z)=-(Y,Z), negative values of ρ will not be further considered in this section. 



 §.§ Asymptotically optimal transform  in the SquareV model


Using (<ref>) and (<ref>), we obtain the following expressions for the joint moments of (Y,Z) as defined in (<ref>):   

    μ_ij   =1+ρ/4(1+(-1)^i+j) +1-ρ/4((-1)^i+(-1)^j)  
       =(1-ρ) 1+(-1)^i/2 1+(-1)^j/2 
       +ρ 1+(-1)^i+j/2

for all i,j=0,1,… and all ρ∈[0,1). 


As the result, ODE (<ref>) becomes  

    ψ”(ρ)/ψ'(ρ)=
    	q_z -2ρ/1-ρ ^2,

where 

    q_z:=1/3z^2-1/3.

We see that formulas (<ref>)–(<ref>) are rather similar to 
(<ref>)–(<ref>). 
Hence, quite similarly to (<ref>) and (<ref>), 
here we have 

    ψ'(ρ)=(1 - ρ^2)^q_z

and 

    ψ(ρ)=ψ_4;z(ρ):=∫_0^ρ(1 - r^2)^q_z dr
    	=ρ  _2F_1(12,-q_z;32;ρ ^2).

Here the subscript 4 in ψ_4;z refers to the four points of the distribution of the random point (Y,Z) in the SquareV model, currently under consideration; thus, one can distinguish between the function ψ_4;z in (<ref>) and the function ψ_z in (<ref>). 

Accordingly, recalling again (<ref>), we see that in the SquareV model the transform 

    Ψ_4;(R)=ψ_4;z_(R)

of Pearson's R 
is asymptotically optimal for any given significance level ∈(0,1). 

In particular, choosing

    z=1,

we have q_z=0. Hence, in view of the integral expression in (<ref>), ψ_4;z(ρ)=ρ, so that the family (Ψ_4;(R))_∈(0,1) of transforms of Pearson's statistic R includes R itself. 
More specifically, 

    R=Ψ_4;(R)  for  =1-Φ(1)≈0.159.


However, in order for Fisher's transform R_F of R to belong to the family (Ψ_4;(R))_∈(0,1) of asymptotically optimal transforms of R in the SquareV model, one would have to have q_z=-1 for some real z, which is impossible, because, in view of (<ref>), q_z is always greater than -1/3. 
So, in contrast with the BVN model (where, according to (<ref>), R_F is asymptotically optimal in the limit case with =0 and z_=∞), in the SquareV model Fisher's transform R_F is not asymptotically optimal for any significance level ∈[0,1], even if the endpoints =0 and =1 are included as limit cases.  



 §.§ The transform  in the SquareV model is asymptotically better than  and  in wide ranges of values of the significance level, including  itself
 


This subsection is similar to Subsection <ref>.
One can see that in the SquareV model 

    _ψ(R)(z)=ρ/3√(1-ρ ^2) (z^2-1)
    -z^2/2 √(1-ρ ^2) ψ”(ρ)/ψ'(ρ).


In particular, 

    _R(z)=ρ/3√(1-ρ ^2) (z^2-1),
 

    _R_F(z)=_ψ_∞(R)(z)=-ρ/3√(1-ρ ^2) (2z^2+1),

where ψ_∞ is as defined in (<ref>), and 

    _Ψ_4;(R)(z)=_ψ_4;z_(R)(z)
    =ρ/3√(1-ρ ^2) (z^2/z_^2-1).
 

So, 
|_Ψ_4;(R)(z_)|<|_R_F(z_)| for ρ0 whenever z_≥1/√(2) or, equivalently, ∈(0,1-Φ(1/√(2))), with 1-Φ(1/√(2))=0.2397…. 
Therefore, 
the transform Ψ_4;, which is asymptotically optimal for the given significance level ∈(0,1), will still be asymptotically better than Fisher's transform R_F for any significance level ∈(0,0.5) provided that ∈(0,0.2397). 

As for the comparison of the asymptotically optimal transform Ψ_(R) with R itself, we can similarly see that, for instance, if =0.05, then Ψ_(R) will be asymptotically better than R, not just for the significance level =0.05, but for any significance level ∈(0,0.11344); if =0.01, then Ψ_(R) will be asymptotically better than R for any significance level ∈(0,0.096927). 

We see that the advantage of the asymptotically optimal transform of R over R itself is substantially less in the SquareV model than in the BVN model. 
Vice versa, the advantage of the asymptotically optimal transform of R over the Fisher transform R_F of R is much greater in the SquareV model than in the BVN model. 




 §.§ Comparison of asymptotically optimal transform  in the SquareV model with  and  for finite sample sizes
 


Computer simulations 
mainly similar to those for finite sample sizes in the BVN model were carried out for the SquareV
model as well, the main difference with the BVN case being the use of the Mathematica command  instead of 
. However, instead of about 3.6 hours for the BVN model, the simulations of the same number, 
12×10^6×2×4×(10+10^2+10^3+10^4)≈10^12, of realizations of the random pair (Y,Z) in the SquareV model took about 17 hours. It appears that the command , used to simulate distributions on finite sets, is not as highly optimized as  . 

After these simulations for the SquareV model had been completed, it was realized that the SquareV model can obtained by a transformation and a subsequent re-parametrization of the BVN model, as follows: if a random pair (U,V) has the BVN distribution with zero means, unit variances, and an arbitrary correlation coefficient þ∈(-1,1), then the random pair 

    (Y,Z):=( U, V)

will have the SquareV distribution described by formula (<ref>), provided that the parameters þ and ρ are related by the formula 

    (-1,1)∋cos[π2 (1-ρ)]=þ⟷ρ=1-2π arccosþ∈(-1,1).

This can be shown, for instance, by a standard change of variables in the relevant integrals to reduce the BVN density with correlation coefficient þ to the standard BVN density, and then integrating in polar coordinates.  

So, the simulations for the SquareV model were mainly reduced to those for the BVN model, with the reduction of the execution time from about 17 hours to about 4 hours. Thereby, altogether 2×12×10^6×2×4×(10+10^2+10^3+10^4)≈2×10^12 realizations of the random pair (Y,Z) in the SquareV model were obtained. 
Thus, one was able to reduce the corresponding estimates s̃_·(,ρ,n) of the standard error of the  estimates _·(,ρ,n) of the relative error of the approximation of  using R, R_F, and Ψ_4;(R) by a factor of about √(2)≈1.4.  

Except for the differences described above, the format of Table <ref> for the SquareV model is similar to that of Table <ref> for the BVN model. Correspondingly, Figs. <ref> and <ref> for the SquareV model are similar to  Figs. <ref> and <ref> for the BVN model. 

 


Looking at Fig. <ref>, we see that generally Fisher's transform R_F performs much worse than R and Ψ_4;(R), at least when ρ∈{0.1,0.5,0.9} and n∈{10^2,10^3,10^4}. Therefore, to make the comparisons between R and Ψ_4;(R) clearer, in Fig. <ref> the graphs corresponding to R_F were removed. 

Somewhat similarly to the case of the BVN model, here we can see from Fig. <ref> that, 
if n=10 or if ρ∈{0,0.1}, the asymptotically optimal transform Ψ_4;(R) appears to have comparatively little or no advantage over R itsel;  
however, for ρ∈{0.5,0.9} and n∈{10^2,10^3,10^4}, the transform Ψ_4;(R) appears to be mostly better or much better than R itself. 



§ CONCLUSION
 
The main result of this paper is Theorem <ref>, which shows for any correlation-parametrized model of dependence and for any given significance level ∈(0,1), there is an asymptotically optimal transform of Pearson's correlation statistic R, for which the generally leading error term for the normal approximation vanishes 
for all values ρ∈(-1,1) of the correlation coefficient. 

It is also shown that in the BVN model Pearson's R turns out to be asymptotically optimal for a rather unusual significance level ≈0.240, whereas Fisher's transform R_F of R is asymptotically optimal for the limit significance level =0. 
In the other specific model of dependence considered in this paper – the SquareV model, Pearson's R is asymptotically optimal for a still rather high significance level ≈0.159, whereas Fisher's transform R_F of R is not asymptotically optimal for any ∈[0,1]. 

Moreover, we saw that in both the BVN model and the SquareV model, 
the transform Ψ_(R), asymptotically optimal for a given value of , is in fact asymptotically better than R and R_F in wide ranges of values of the significance level, including  itself. 

Recall that Fisher's transform R_F of R was designed for the BVN case, with the purpose of making the asymptotic variance constant with respect to the correlation coefficient ρ. That R_F usually turns out to be asymptotically closer to normality than R in the BVN  model might now be explained by the observation that the significance level =0 (for which R_F is asymptotically optimal in the BVN case) is closer to such usual in statistical practice values of the significance level as 0.05 than to the significance level _R≈0.240 (for which R is asymptotically optimal in the BVN case). 

Extensive computer simulations for the BVN and SquareV models of dependence presented in this paper suggest that, for sample sizes n≥100 and significance levels ∈{0.01,0.05}, the mentioned asymptotically optimal transform of R generally outperforms both Pearson's R and Fisher's transform R_F of R, the latter appearing generally much inferior to both R and the asymptotically optimal transform of R in the SquareV model.  




































imsart-number














