Validation of positive quadrant dependence

Insurance: Mathematics and Economics 56 (2014) 38–47

Contents lists available at ScienceDirect

Insurance: Mathematics and Economics journal homepage: www.elsevier.com/locate/ime

Validation of positive quadrant dependence Teresa Ledwina a , Grzegorz Wyłupek b,∗ a

Institute of Mathematics, Polish Academy of Sciences, ul. Kopernika 18, 51-617 Wrocław, Poland

b

Institute of Mathematics, University of Wrocław, pl. Grunwaldzki 2/4, 50-384 Wrocław, Poland

highlights • • • •

We provide a reparametrization of positive quadrant dependence. We propose a visualization of dependence structure of data and copula models. We introduce two tests for positive quadrant dependence. We find interesting properties of finite sample distribution of empirical copula.

article

info

Article history: Received November 2013 Received in revised form January 2014 Accepted 25 February 2014 JEL classification: C1 Keywords: Concordance ordering Copula Correlation Order preserving Rank test

abstract Quadrant dependence is a useful dependence notion of two random variables, widely applied in reliability, insurance and actuarial sciences. The interest in this dependence structure ranges from modeling it, throughout measuring its strength and investigations on how increasing the dependence effects of several reliability and economic indexes, to hypothesis testing on the dependence. In this paper, we focus on testing for positive quadrant dependence. We propose two new tests for verifying positive quadrant dependence. We prove novel results on finite sample behavior of power function of one of the proposed tests as well as evaluate and compare the two new solutions with the best existing ones, via a simulation study. These comparisons demonstrate that the new solutions are slightly weaker in detecting positive quadrant dependence modeled by classical bivariate models and outperform the best existing solutions when some mixtures, regression and heavy-tailed models have to be detected. Finally, the methods introduced in the paper are applied to real life insurance data, to assess the dependence and test them for positive quadrant dependence. © 2014 Elsevier B.V. All rights reserved.

1. Introduction and motivations There is growing evidence that in many areas of sciences and their applications such as insurance, mathematical finance and reliability the analysis of dependence structure cannot be neglected. For extensive discussion see Albers (1999), Denuit and Scaillet (2004), Dhaene et al. (2009) and references therein. Simultaneously, there is strong emphasis on going beyond classical approach based on linear correlations and normally distributed observations. See above, Cont (2001), Embrechts et al. (2002), and Kallenberg (2009) for further information. Among many different concepts of dependence of two random variables positive quadrant dependence (PQD in short) plays an

important role. It is not a very restrictive notion, has simple and intuitive probabilistic interpretation and induces natural ‘more positive’ ordering which is preserved by many important economic indexes and parameters; see Scaillet (2005) and Dhaene et al. (2009) for several comments and examples. We shall consider testing for PQD of two real random variables X and Y with joint distribution function H (x, y) and continuous marginal distribution functions F and G. We like to verify the null hypothesis

H0 : H (x, y) ≥ F (x)G(y) for all x, y ∈ R against the alternative

A0 : H (x, y) < F (x)G(y) for some x, y ∈ R. H0 means that the joint probability that X and Y are both small is

∗

Corresponding author. Tel.: +48 71 375 74 41. E-mail addresses: [email protected] (T. Ledwina), [email protected] (G. Wyłupek). http://dx.doi.org/10.1016/j.insmatheco.2014.02.008 0167-6687/© 2014 Elsevier B.V. All rights reserved.

larger than the corresponding one when X and Y are independent. The approach is nonparametric and we do not assume any prescribed parametric model for H. Since F and G are continuous

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

then there exists uniquely determined distribution function C on the unit square [0, 1]×[0, 1], called the copula, such that H (x, y) = C (F (x), G(y)). The function C (u, v) equals H (F −1 (u), G−1 (v)), u, v ∈ [0, 1], where F −1 is the generalized inverse of F defined as F −1 (u) = inf{t ∈ R : F (t ) ≥ u}, u ∈ [0, 1]. G−1 is defined similarly. Set U = F (X ), V = G(Y ). Then (U , V ) obeys the bivariate distribution function C . Kimeldorf and Sampson (1978) named it uniform representation of H. Additionally, set I (u, v) for the copula in case of independence of U and V . Obviously, I (u, v) = uv . In consequence, H0 and A0 can equivalently be stated as

H0 : C (u, v) ≥ I (u, v) for all u, v ∈ [0, 1], Introducing nonparametric approach and copulas allows us to restrict attention to studying PQD for distributions with the same marginals. In such a setup PQD is labeled as more concordant in Tchen (1980) and in many related papers, and more correlated in Dhaene et al. (2009). Testing H0 against A0 has recently been considered by Denuit and Scaillet (2004), Scaillet (2005), Gijbels et al. (2010), and Gijbels and Sznajder (2013a). Roughly speaking, the approach to construct test statistics in these papers is built upon a measure of discrepancy between a copula estimator and the independence copula I. Our starting point to test’s construction is different and is motivated by the approach elaborated in Kallenberg and Ledwina (1999), Janic-Wróblewska et al. (2004), and Ledwina and Wyłupek (2012a). Namely, we investigate the structure of C (u, v) and its relation to I by looking at the successive correlation coefficients of fr (U ) and fs (V ) for a rich class of functions fj ’s in L2 ([0, 1], du). To be specific, recall that Kallenberg and Ledwina (1999) considered independence testing. To enable a detection of higher order correlations, not a linear one, only, they found it useful to look at the theoretical and empirical correlation coefficients of functions br (U ) and bs (V ), where {bj }j≥1 are the orthonormal Legendre polynomials on [0, 1]. Moreover, in this system, for large class of underlying C ’s, independence is equivalent to CorrC [br (U ), bs (V )] = 0, r , s ≥ 1. This allowed them to replace testing of independence via checking whether the consecutive correlations between br (U ) and bs (V ), r , s = 1, 2, . . . , are 0. Besides, the first correlation coefficient (r = s = 1) is Spearman’s measure of linear dependence of U and V . In this way a much more flexible than Spearman’s rho solution was constructed. In Janic-Wróblewska et al. (2004) testing independence against PQD was considered. Then the system {bj }j≥1 is no longer useful, as the alternative is restricted one and, by the known result on PQD, it can be characterized by covariances in a system of non-decreasing functions. Such a system was proposed in that paper, but it does not provide a simple interpretation of the resulting correlations. Therefore, in this contribution we shall rely on a new system {lj }j≥1 , with the functions lj ’s patterned after Ledwina and Wyłupek (2012a,b). Some details are as follows. Recall that H0 is equivalent to the relation CovC [f1 (U ), f2 (V )] ≥ 0 for all non-decreasing f1 , f2 : [0, 1] → R, for which the covariance exists. We restrict attention to some selected net of such functions. Namely, to the non-decreasing functions lj ’s given by lj (z ) = −

1 − pj pj

1(0 ≤ z < pj ) +



pj 1 − pj

1(pj ≤ z ≤ 1),

(1)

where 1(A) is the   indicator of the set A while pj ∈ (0, 1). Note that lj (z )dz = 0, l2j (z )dz = 1, j = 1, 2, . . .. With such a choice of the underlying functions, for each r , s ≥ 1 set

ρr ,s = CovC [lr (U ), ls (V )].

Then

ρr ,s = CorrC [lr (U ), ls (V )] = EC [lr (U ) · ls (V )] = √

C (pr , ps ) − pr ps

pr ps (1 − pr )(1 − ps )

.

(2)

The last equation follows by immediate calculation of the integral EC [lr (U ) · ls (V )]. Details are given in the Appendix. If the points (pr , ps ), r , s ≥ 1, lie densely in [0, 1] × [0, 1] then, by uniform continuity of C (cf. Nelsen, 2006, p. 11), H0 and A0 can equivalently be expressed via the correlation coefficients as

H0 : ρr ,s ≥ 0 for all r , s ≥ 1,

A0 : C (u, v) < I (u, v) for some u, v ∈ [0, 1].



39

(3)

A0 : ρr ,s < 0 for some r , s ≥ 1. Obviously, the chain of relations (2) allows for equivalent phrasing the hypotheses in terms of the correlation coefficients of lr (U ) and ls (V ), or the Fourier coefficients of C in the system {lr · ls }r ,s≥1 , or in terms of the standardized differences C (pr , ps ) − pr ps , as well. Since for any copula it holds that C (u, v)− uv = C (u, v)[1 − u −v+ C (u, v)]−[u − C (u, v)][v− C (u, v)], therefore ρr ,s ≥ 0 if and only if the cross product ratio of probabilities of the respective quadrants, given by C (pr , ps )[1 − pr − ps + C (pr , ps )]/[pr − C (pr , ps )][ps − C (pr , ps )] is greater than or equal to 1. The reformulation of the testing problem opens some ways to solve it or to phrase existing solutions in these terms. To construct a test for (3) one can plug some empirical copula into (2) which shall provide the related empirical correlation coefficients. Then one of the traditional approaches is to consider a sum of squares or some quadratic form of the empirical coefficients. The solutions as that in Denuit and Scaillet (2004) could be included into this group. Recently developed approach by Gijbels et al. (2010) and Gijbels and Sznajder (2013a) can also be seen to follow such an idea in some refined form. The simplest variant is, roughly, the sum of squares of the empirical correlations calculated at jump points of some empirical copula. The solution is an example of a very useful variant of the integral Anderson–Darling statistic. See Section 4 for details. Another idea in this stream could be a natural extension of the solution QT elaborated in Ledwina and Wyłupek (2012a, 2013). Due to the reparametrization (3), still another of possible approaches to verify H0 is via multiple comparisons and we shall develop this idea in this contribution. Namely, we focus on two variants which can be viewed as infimum versions of the Anderson–Darling statistic. More precisely, in view of the form of A0 and (2), both new solutions shall be minima of the standardized differences between the empirical and independence copulas. The solutions shall differ by the choice of pj ’s. The first one can be seen to be a natural analogue of the Md statistic in Ledwina and Wyłupek (2012a) and we retain the same notation for it. In the second solution we consider random pj ’s being, as in the case of the above discussed integral Anderson–Darling statistic, the respective ranks corresponding to the jump points of the empirical copula. This statistic shall be denoted by AD∗ . Details are collected in Section 2. The paper is organized as follows. In Section 2, we describe the two new test statistics Md and AD∗ . Section 3 presents finite and large sample results on Md . In particular, important and interesting result on controlling, under fixed sample size n, the error of the first kind for the test based on Md is proved. This result is based on a new finding on finite sample distribution of empirical copula; cf. Theorem 1. It is also shown that the test based on Md is consistent against any sufficiently regular alternative from A0 , provided that the smoothing parameter d depends on n and tends to infinity at some not very fast rate. The related questions on AD∗ have not been solved yet and shall be a subject of further work. Section 4 reports the results of the extensive simulation study. We compare three tests already considered in Gijbels et al. (2010) with the two new

40

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

solutions, under several known and popular copula models. Our simulations confirm that the integral Anderson–Darling statistic works very well for dependence structures defined via many classical copulas. However, in cases of certain mixture models, some regression models, t-copulas and some other heavy-tailed bivariate models, the infimum version AD∗ of the Anderson–Darling statistic and Md are more sensitive. In Section 5, the set of Danish fire insurance data is examined for the presence of PQD. We give some concluding remarks in Section 6. Descriptions of models used in our simulation experiments are gathered in the Appendix. Moreover, derivation of the formula (2) is provided there. 2. New solutions Throughout (X1 , Y1 ), . . . , (Xn , Yn ) are independent identically distributed random vectors drawn from two-dimensional distribution function H. Furthermore, Ri is the rank of Xi , i = 1, . . . , n, in the sample X1 , . . . , Xn whereas Si denotes the rank of Yi , i = 1, . . . , n, in the sample Y1 , . . . , Yn . 2.1. Estimates of ρr ,s ’s via classical empirical copula Standard definition of the empirical copula Cn is as follows Cn (u, v) =

n 1   Ri

n i=1

1

n

 Si ≤ u, ≤ v , n

u, v ∈ [0, 1].

L(r , s) =

√

(4)

(5)

nρˆ r ,s

1 pr ps (1 − pr )(1 − ps )

Rn (r , s),

1≤r ,s≤d

(7)

is too small. L(r , s) is defined in (5) while d is a given fixed number of dyadic points in which we control the difference in (2). We allow d to depend on the sample size n. Note that we consider the critical region of the form {Md < c } for some c and calculate finite sample quantiles of Md as in Ledwina and Wyłupek (2012a). See p. 735 there for more details. The statistic Md is based on the components which have nice and useful interpretation; see Propositions 2 and 3, Corollary 1 and Figs. 1–4. In the simulations and examples we fixed d in Md to 31, i.e. the statistic inspects 961 empirical correlations. It is also possible to exploit Md to construct some, possibly more specialized, data driven statistic for H0 . For constructions of such kind, for case of verification of stochastic order, see Ledwina and Wyłupek (2012a, 2013)). They were denoted by QT , QT 1 , . . . therein. We are not elaborating this idea here. Instead, in Section 2.3 we propose an alternative omnibus solution AD∗ , which checks differences in (2) only in n jump points of slightly refined empirical copula.

Now, given n, we take in (1) pj = j/(n + 1). Moreover, as in Gijbels et al. (2010), instead of Cn (u, v), in this construction we consider C˜ n (u, v) =

for a natural, from the testing viewpoint, estimator of ρr ,s and its standardized version, respectively. The correction for continuity in (5) plays the same role as the transformation Ri /n → Ri /(n + 1) applied e.g. in Gijbels et al. (2010). Namely, it improves finite sample behavior of ρˆ r ,s . Such a correction was applied throughout the book of Behnen and Neuhaus (1989) and in some our earlier papers. It is easy to see that ρˆ 1,1 is a variant of the Blomqvist rank statistic; see Blomqvist (1950) and Yanagimoto and Okamoto (1969). In general, ρˆ r ,s is a linear combination of numbers of pairs of ranks (Ri , Si ) falling into the respective four sets, naturally defined by lr · ls , cf. the Appendix for the related details. In consequence, ρˆ r ,s can be seen to be a linear function of the number of pairs of ranks (Ri , Si ) falling into the quadrant of the form [1, npr ]×[1, nps ]. More precisely,

ρˆ r ,s = √

Md = min L(r , s)

2.3. The second solution

Set n 1   Ri − 0.5   Si − 0.5  ρˆ r ,s = lr ls and n i =1 n n

Our first test rejects H0 if

(6)

Rn (r , s) = Cn (p∗r , p∗s ) − pr ps − ps (p∗r − pr ) − pr (p∗s − ps ) while p∗j = p∗j (n) = (⌊npj + 0.5⌋ − 1)/n if npj + 0.5 is a natural number, otherwise p∗j = p∗j (n) = (⌊npj + 0.5⌋)/n. The symbol ⌊•⌋ denotes the integer part of •. 2.2. The first solution In this section we follow one of the approaches developed in Ledwina and Wyłupek (2012a) in case of testing for stochastic dominance. We consider lj ’s in (1) with pj = aj , where a1 , a2 , . . . are the successive points of the form (2i − 1)/2k+1 , k = 0, 1, . . . , i = 1, 2, . . . , 2k . Such lj ’s can be seen to be the (standardized) projections of elements of the Haar basis into the cone of all nondecreasing functions; cf. Ledwina and Wyłupek (2012b) for more comments.

n  1   Ri Si 1 ≤ u, ≤v . n i =1 n+1 n+1

(8)

Here the idea is to measure minimal standardized difference between C˜ n and I only in jump points of C˜ n . Therefore, we define

√   Rj Sj  , n C˜ n j=1,...,n n+1 n+1   Rj Sj Rj Sj  Rj  Sj  − / 1− 1− . (9) n+1n+1 n+1n+1 n+1 n+1

AD∗ = min

As in the case of Md , we reject H0 for small values of AD∗ , i.e. if {AD∗ < c } for some c. The statistic AD∗ exploits more efficiently available information on empirical copula than Md does. Moreover, such a construction can be particularly useful for some potential extensions for dimensions higher that two. Note also that introducing absolute value in the definition of AD∗ and changing min into max yields a new test of independence. On the other hand, AD∗ is a more complicated object than Md and its investigation deserves separate careful considerations. 3. Some theoretical results For fixed r , s, ρˆ r ,s of the form (5) is the so-called linear rank statistic and a rejection region of the form {ρˆ r ,s < c } provides a good test for the one-sided testing problem ρr ,s ≥ 0 versus ρr ,s < 0. Obviously, there is vast literature on linear rank tests for independence. For an illustration see Behnen (1971), Ruymgaart (1978), and references therein. However, we are interested in some specific properties of such tests and therefore some further considerations are needed. In Section 2.1 we noticed that ρˆ r ,s is a linear function of Cn (p∗r , p∗s ). Asymptotic behavior of Cn has extensively been studied; see Fermanian et al. (2004) for recent account and the related literature. In contrast, some useful finite sample properties of the process {Cn (u, v), u, v ∈ [0, 1]} and some related test statistics have not been established yet, to the best of our knowledge. In particular, the question is: whether and when H0 implies that

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

finite sample distributions of the process {Cn (u, v), u, v ∈ [0, 1]} are stochastically ordered in a useful sense and for which test statistics Φ (Cn ) and the related critical regions Kn , say, under C ∈ H0 and fixed n it holds that PC (Φ (Cn ) ∈ Kn ) ≤ PI (Φ (Cn ) ∈ Kn ).

introduce H∗ and H∗∗ being the cumulative distribution functions corresponding to C∗ and C∗∗ , respectively. Hence, PC∗ (Q11 ≥ q11 ) =

 =

n



w1

R2n

Proposition 1. Suppose that a copula C∗ has larger quadrant dependence than a copula C∗∗ , i.e. (11)

Then, under any fixed u1 , v1 ∈ [0, 1], any c11 ∈ R, and any n it holds that (12)

Proof. Set w1 = ⌊nu1 ⌋, z1 = ⌊nv1 ⌋, q11 = nc11 . Let #[A] stand for the number of elements of the finite set A. Define

(13)

A(J1 , w1 ) = { x : Ri1 = 1, . . . , Riw1 = w1 }

(15)





n

w1

R2n

1Ao (x) 1B o (y )

n 

dH∗∗ (xi , yi )

i=1



The property (12) means that, under the underlying C∗ satisfying (11), the random variable Cn (u1 , v1 ) is stochastically larger than under C∗∗ . This immediately implies that to control the error of the first kind of the critical region {ρˆ r ,s < c }, for any fixed sample size n, it is enough to do it under I. Since, however, we like to combine the statistics ρˆ r ,s , r , s ≥ 1, or, equivalently, the values of the empirical copula in some points, therefore we need a stronger result than (12). For this purpose, along with (12), we prove the following statement. Theorem 1. Assume (11). Then for any d1, d2 ≥ 1, any 0 ≤ u1 < · · · < ud1 ≤ 1, 0 ≤ v1 < · · · < vd2 ≤ 1, and arbitrary cij , i = 1, . . . , d1, j = 1, . . . , d2, it holds that



d1  d2  

Cn (ui , vj ) ≥ cij





≥ PC∗∗

 d1  d2  

Cn (ui , vj ) ≥ cij



 .

(16)

i=1 j=1

Proof. Analogously as above, set wi = ⌊nui ⌋, zj = ⌊nvj ⌋, qij = ncij and Qij = #[ (Rt , St ) : Rt ≤ wi , St ≤ zj ]. Then (16) is equivalent to

and



B (J1 , z1 , q11 ) = { y : #[ Sij : ij ∈ J1 , Sij ≤ z1 ] ≥ q11 }. Note that, given w1 , for different possible choices of 1 ≤ i1 < · · · < iw1 ≤ n the related sets A(J1 , w1 )’s are disjoint. Moreover,

A(J1 , w1 ) ∩ B (J1 , z1 , q11 ).

dH∗ (xi , yi ).

i=1 j=1

We shall justify (13) using the approach developed in Yanagimoto and Okamoto (1969). Namely, given w1 , put J1 = J1 (i1 , . . . , iw1 ) for a set of indexes i1 , . . . , iw1 satisfying 1 ≤ i1 < · · · < iw1 ≤ n. Next set x = (x1 , . . . , xn ), y = (y1 , . . . , yn ) and define

{Q11 ≥ q11 } =

n  i=1

This concludes the proof.

With these notations (12) reads as



PH∗ (Ao ∩ B o )

1Ao (x) 1B o (y )

PC∗ (Q11 ≥ q11 ) ≥

P C∗

Q11 = #[ (Rt , St ) : Rt ≤ w1 , St ≤ z1 ].

PC∗ (Q11 ≥ q11 ) ≥ PC∗∗ (Q11 ≥ q11 ).

w1

= PC∗∗ (Q11 ≥ q11 ).

To formulate our first result, recall that ρˆ 1,1 is a linear function of the number of pairs of ranks (Ri , Si ) falling into the quadrant [1, (n + 1)/2) × [1, (n + 1)/2). This observation, along with the argument of the proof of Theorem 8.1 in Yanagimoto and Okamoto (1969), leads us to the following interesting conclusion.

PC∗ (Cn (u1 , v1 ) ≥ c11 ) ≥ PC∗∗ (Cn (u1 , v1 ) ≥ c11 ).



n

Now observe that both indicators in (15) are concordant with respect to the ith coordinate for each i ≤ n i.e. considered as functions of the ith coordinate (with all other coordinates held fixed), they are monotone in the same direction. More precisely, when i ≤ w1 then both indicators are non-increasing while for i > w1 both are non-decreasing. Moreover, EH∗ 1Ao (X ) = EH∗∗ 1Ao (X ) and similar relation holds for 1B o (Y ) as the distribution of ranks is independent on parent distribution of observations. Hence, Proposition 7.4(i) of Yanagimoto and Okamoto (1969) yields EH∗ 1Ao (X )1B o (Y ) ≥ EH∗∗ 1Ao (X )1B o (Y ) and, in consequence,

3.1. Finite sample results on Cn

C∗ (u, v) ≥ C∗∗ (u, v) for all u, v.



(10)

Note that the last question has already been addressed and solved for some interesting classes of Φ ’s and for stronger orders than that implied by positive quadrant dependence. For some evidence see Lehmann (1966), Yanagimoto and Okamoto (1969), and Schriever (1987), to mention basic references. Schriever (1985) speculated about such a question for PQD and his conclusion was very pessimistic: ‘‘. . . this ordering is unsatisfactory because the order preserving property does not carry over to samples in any useful sense’’.; cf. p. 62 ibidem. Below, we show that it was too far reaching inference. The property (10) is very useful for controlling the error of the first kind of Kn , as, under I, Cn is distribution free. It implies that, for every fixed n, the critical value corresponding to Kn and found by a simple simulation experiment under I is valid for H0 .

41

(14)

all J1 ’s

Since (X1 , Y1 ), . . . , (Xn , Yn ) are independent and identically distributed then each summand in (14) has the same probability equal to that for the choice i1 = 1, . . . , iw1 = w1 . Therefore, set Ao = A(J1 (1, . . . , w1 ), w1 ), B o = B (J1 (1, . . . , w1 ), z1 , q11 ). Besides,

P C∗

d1  d2  

 Qij ≥ qij

i=1 j=1



 ≥ PC∗∗

d1  d2  

Qij ≥ qij

 

.

(17)

i =1 j =1

Given w1 < · · · < wd1 , denote by J1 = J1 (i1 , . . . , iw1 ), J2 = J2 (iw1 +1 , . . . , iw2 ), . . . the consecutive sets of blocks of indexes

1 ≤ i1 < · · · < iw1 < iw1 +1 < · · · < iw2 < · · · < iwd1 ≤ n and set J = (J1 , . . . , Jd1 ), w = (w1 , . . . , wd1 ), z = (z1 , . . . , zd2 ), and q = (q11 , . . . , qd1d2 ). With these notations define

A(J , w ) = { x : Ri1 = 1, . . . , Riw1 = w1 , Riw

1 +1

= w1 + 1, . . . , Riwd1 = wd1 }

42

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

Proof. Set

and

B ( J, z , q ) =

d2  d1 



 y : # ik ∈

j=1 i=1

i 

 Jt : Sik ≤ zj



 w(r , s) = 1/ pr ps (1 − pr )(1 − ps ).

≥ qij .

By the definition of pj and Cn one obtains

t =1

Note that, given w, for the successive possible J ’s the resulting sets A(J , w )’s are disjoint. Besides, F and G are continuous and

all J ’s

Similarly as above, consider one particular J , say, J o , defined by the indexes 1, 2, . . . , wd1 . Then

 P C∗

d1  d2  {Qij ≥ qij }



=



n

wd1

R2n

1A(Jo ,w ) (x) 1B (Jo ,z ,q) (y )

n 

dH∗ (xi , yi ).

i =1

The rest of the argument is analogous as in the case of (15).



The relation (16) can be read as: the vector (Cn (u1 , v1 ), . . . , Cn (ud1 , vd2 )), under the underlying distribution C∗∗ , is smaller in the upper orthant order than the same vector in the situation when the underlying distribution is C∗ . For the notion of orthant order and its properties see Shaked and Shanthikumar (2007), Section 6.G.1. 3.2. Finite sample results on L(r , s) and Md Corollary 1. By (6) and Theorem 1, for any C in H0 , any n, any c ∈ R, any fixed r , s, and any d ≥ 1, we get PC (L(r , s) < c ) ≤ PI (L(r , s) < c ) and PC (Md < c ) ≤ PI (Md < c ). So, the tests based on L(r , s) and Md satisfy (10). Another example is a variant of the Kolmogorov–Smirnov statistic defined below and patterned after Scaillet (2005). A statistic of such kind shall be considered in our simulation study. To be specific, under notations of Theorem 1 set

√ KS =

max

1≤i≤d1,1≤j≤d2

2 n

  1 1 ≤ Cn (ar , as ) − Cn ar + , as + n

n{ui vj − Cn (ui , vj )}.

Corollary 2. For any C in H0 , any n, and any c ∈ R it holds that PC (KS ≥ c ) ≤ PI (KS ≥ c ). This shows that critical values of KS can be obtained by a simple simulation experiment, only. 3.3. Asymptotic results on L(r , s) and Md In this section, in both L(r , s) and Md , we restrict attention to pj = aj , where aj ’s result from the dyadic partition of [0, 1]. We √ start with a simple result on ρˆ r ,s = L(r , s)/ n, showing that this estimate of ρr ,s is strongly consistent.

with probability one.

n

1 ≤ r , s ≤ d,

(19)

while ρˆ r ,s = w(r , s){Cn (ar , as ) − ar as } + o(1). Since Cn (u, v) can be expressed via bivariate empirical distribution function of the sample (X1 , Y1 ), . . . , (Xn , Yn ), therefore the proof is concluded by an application of the Glivenko–Cantelli theorem in two dimensions. For details of such an argument see Groeneboom et al. (1976), p. 121. Next proposition, though asymptotic by nature, shall be useful in interpreting finite sample results. Proposition 3. Under C = I, every fixed r , s, and n → ∞ it holds that D

L(r , s) → N (0, 1),

(20)

D

where → denotes the convergence in distribution. Proof. By (6), the statement (20) follows, for example, from Theorem 3 of Fermanian et al. (2004); cf. p. 851 and p. 854, and the form of the asymptotic covariance function, given on p. 851 ibidem, which implies that, under √ independence, ar as (1 − ar )(1 − as ) is the asymptotic variance of nRn (r , s).  The statement (20) along with Corollary 1 make L(r , s) diagnostic tool in detecting departures from H0 , as then a significance of negative values of L(r , s) can be easily estimated. Next result gives rough information on a rate in which Md tends to −∞ under H0 and n → ∞, in the case when d = d(n) and d(n) → ∞. More precisely, for some sequence of natural k(n)’s now we consider splitting [0, 1] into intervals with usage of 2k(n)+1 − 1 points aj , j = 1, . . . , 2k(n)+1 − 1, of the form aj = (2i − 1)/2k+1 , k = 0, . . . , k(n), i = 1, . . . , 2k . We set d(n) = 2k(n)+1 − 1 for the number of the resulting points ai ’s. Proposition 4. Under C = I, d = d(n) → ∞, as n → ∞, implies min

L(r , s) = Md(n) = OP (d(n)).

(21)

Proof. Note that ai satisfies 1/[d(n) + 1] ≤ ai ≤ d(n)/[d(n) + 1] and the same bound is valid for 1 − ai . Hence, for each 1 ≤ r , s ≤ d(n), it holds that

w(r , s) ≤ 2[d(n) + 1].

(22)

By (6) and (19), we infer that for |Md(n) | = | max1≤r ,s≤d(n) {−L(r , s)}| it holds that

|L(1, 1)| ≤ |Md(n) | ≤

√ n[ max

1≤r ,s≤d(n)

w(r , s)]

| Cn (ar , as ) − ar as + Rn∗ (r , s) | √ √ ≤ 2[d(n) + 1][ n sup |Cn (u, v) − uv| + 8/ n ]. ×

Proposition 2. If H is continuous distribution function on R × R then, as n → ∞, it holds that

n

Hence, the expression Rn (r , s), appearing in (6), has the structure Cn (ar , as ) − ar as + Rn∗ (r , s), where

1≤r ,s≤d(n)

ρˆ r ,s → ρr ,s

n

n

|Rn∗ (r , s)| ≤ 8/n,

i=1 j=1



−

≤ Cn (ar , as ) − Cn (p∗r , p∗s )   2 4 2 ≤ . ≤ Cn (ar , as ) − Cn ar − , as −

d1  d2   {Qij ≥ qij } = A(J , w ) ∩ B (J , z , q). i =1 j =1

(18)

∗

max

1≤r ,s≤d(n)

0≤u,v≤1

By the above, and Theorem 3 of Fermanian et al. (2004) the conclusion follows. 

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

The rate in (21) is presumably highly sub-optimal; cf. Theorem 1 in Einmahl (1996) for an adequate result in the case of weighted empirical process. However, we are not aware of a ready counterpart of this theorem in the case of weighted copula process, while such extension is beyond the scope of this contribution. On the other hand, Proposition 4 is sufficient to get qualitative results on the behavior of Md(n) , such as consistency of the related test, saying that more dense partitioning results in greater precision of inference. Proposition 5. Suppose that the alternative A0 holds true and the copula function C (u, v) has continuous partial derivatives. If d(n) is a non-decreasing sequence, d(n) → ∞ as n → ∞, and d(n) = o(n1/2 ) then the test rejecting H0 for small values of Md(n) is consistent under C . Proof. Since C ∈ A0 and C is uniformly continuous then there exist r0 and s0 such that C (ar0 , as0 ) < ar0 as0 . Let n be so large that both r0 and s0 are smaller than d(n). Given α , let cn,α denote the critical √ value of the test on the level α . Then, by Proposition 4, cn,α = o( n). Moreover, by (6) it holds that PC Md(n) < cn,α





= PC



min

1≤r ,s≤d(n)

L(r , s) < cn,α



  ≥ PC L(r0 , s0 ) < cn,α PC

√

n w(r0 , s0 ){Cn (ar0 , as0 ) − ar0 as0 } < cn,α + o(1)



= PC

n w(r0 , s0 ){Cn (ar0 , as0 ) − C (ar0 , as0 )}

√

 + ∆n (r0 , s0 ) < cn,α + o(1) , where

∆n (r0 , s0 ) =

√

n w(r0 , s0 ){C (ar0 , as0 ) − ar0 as0 }.

√

√

By ∆n (r0 , s0 ) < 0, ∆n (r0 , s0 ) = O( n), cn,α = o( n) and the weak convergence of the empirical copula process, see Theorem 3 of Fermanian et al. (2004), the proof is complete.  4. Simulation study To see how the new solutions Md and AD∗ work in some standard situations and how they compare to some competitors we conducted the extensive simulation study. The screening is partially patterned after the ones reported in Scaillet (2005) and Gijbels et al. (2010). We took into account most of the alternatives from these papers. Moreover, to cover some further important and interesting cases, not included there, we added two additional extreme value models, some regression models and bivariate distributions with heavy tails of marginals. For comparison, we took the two following statistics introduced in Gijbels et al. (2010). The Cramér–von Mises statistic C v M, calculated as C vM =

n   R  Rj Sj Sj 2 j − C˜ n , , n+1n+1 n+1 n+1 + j=1

where { · }+ = max{ · , 0 }, and the integral Anderson–Darling statistic AD, given by

2 S , n+j 1 +   . Rj Sj Rj Sj 1 − 1 − n +1 n +1 n +1 n +1 

AD =

n  j =1

Rj

Sj

n +1 n +1

− C˜ n



Rj

n+1

Gijbels et al. (2010) have proved that, under some assumptions on the underlying C , the tests rejecting H0 for large values of C v M and AD are asymptotically unbiased under H0 and consistent under A0 .

43

Moreover, similarly as in Gijbels et al. (2010), we included the following Kolmogorov–Smirnov test statistic

√ KS =

n

max

i,j=1,...,19

 i

j

20 20

− C˜ n

 i

,

j 

20 20

which returns the largest standardized difference between I and the estimated copula (8) on the grid of points at [0, 1]2 . The grid is made of the values (u, v) evenly spaced inside {0.05, 0.10, . . . , 0.95} × {0.05, 0.10, . . . , 0.95}. The same grid was applied in the construction of the Kolmogorov–Smirnov test in Scaillet (2005). KS rejects H0 for large values. We compared C v M , AD, and KS, under n = 200 and the significance level α = 0.05, with AD∗ and Md with d = 31. As in Gijbels et al. (2010), the critical values were calculated under independence. We used 100 000 repetitions. We have got: KS : 0.6364, C v M : 0.0316, AD : 0.8634, AD∗ : −3.0100, Md : −3.8844. The empirical powers were computed on the basis of 10 000 Monte Carlo runs. All computations were done using R. In the course of calculations of the critical values the seed was set as 1 while the seed 2 was used when powers were computed. Fig. 1 presents the outcomes of the experiments labeled E1 −E19 . In each case we show a heat map of the estimated and standardized correlations L(r , s)’s, cf. (5), for the dyadic grid corresponding to d = 15. Each square of size 0.0625 × 0.0625 represents the respective value of L(r , s) in its upper-right corner. By Propositions 2 and 3, and Corollary 1, the map, in an easily readable manner, shows where is the main source of departures from the null hypothesis. Each separate display is accompanied with the related empirical powers of the five tests under consideration. Parameters of the alternatives were fixed to guarantee the highest power in each case Ej , j = 1, . . . , 19, to lie in [0.6, 0.8]. The first 9 examples are taken from Gijbels et al. (2010). The first four are related to classical copulas of: bivariate normal distribution, Farlie–Gumbel–Morgenstern family, Frank family, and Clayton distribution. Next row contains known mixtures. The third row (E9 − E11 ) illustrates departures from independence of three extreme value distributions: Gumbel, Marshall and Olkin, and Mai and Scherer. To avoid positive dependence their uniform representation (U , V ) is replaced by (U , 1 − V ). The consecutive row (E12 − E15 ) is patterned after Thas and Ottoy (2004) and shows situations related to quadratic, cubic, fourth-order and sinusoidal regression cases, respectively. The last row contains results for skewed t-distribution, symmetric t-distribution and two heavytailed cases. The three last mentioned cases were taken from Kallenberg and Ledwina (1999), where t-distribution is labeled as Pearson type VII. More details on all the cases under consideration are contained in the Appendix. The simulation experiments clearly show some tendencies. AD is the best among the three already elaborated solutions: KS, C v M, and AD. AD works very well when underlying C lies entirely below I (cf. E1 − E4 , E9 − E11 , E13 ). When C and I intersects and the area where C > I is getting larger (white space in the heat maps) then AD is getting weaker. In such situations AD∗ outperforms AD. Md proves to be also useful in the last mentioned situations. On the average, AD∗ proves to be more stable than both AD and Md with d = 31. To see whether the above described tendencies are retained when the sample size is much larger, we calculated empirical powers under n = 1500, α = 0.05, d = 31 in four representative cases: Clayton family, Mardia distribution, quadratic regression ′ and skewed t-distribution. The cases are labeled by E4′ , E6′ , E12 , and ′ E19 , respectively. The parameters of the above families are adjusted to have again maximal empirical powers in [0.6, 0.8]. The results collected in Fig. 2 show that the above described tendencies in the case n = 200 are also present under n = 1500 and are even stronger.

44

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

Fig. 1. Heat maps of the standardized empirical correlation coefficients L(r , s); r , s = 1, . . . , 15 for examples E1 –E19 and related empirical powers. n = 200, α = 0.05.

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

45

′ ′ Fig. 2. Heat maps of the standardized empirical correlation coefficients L(r , s); r , s = 1, . . . , 15 for examples E4′ , E6′ , E12 , E19 and related empirical powers. n = 1500,

α = 0.05.

Table 1 Danish fire insurance data: values of the test statistics for (B, C), (B, P), and (C, P) pairs along with the corresponding p-values, based on 10 000 MC runs. Test

KS

Pair

Obtained value

(B, C) (B, P) (C, P)

1.7266 0.5988 −0.0399

Pair

p-value

(B, C) (B, P) (C, P)

0.0000 0.0804 1.0000

C vM

AD

AD∗

Md

P

S

K

B

Obtained value 0.2932 0.0107 0.0000

12.0429 0.4345 0.0000

−10.6238 −3.6387

−10.8887 −3.7266

0.6094

0.1661

0.0000 0.3154 1.0000

0.0000 0.2518 1.0000

0.0000 0.0253 1.0000

0.0000 0.0395 1.0000

0.5462 0.7816 0.6159

0.1415 0.2842 0.6525

0.0854 0.1944 0.4687

4.7477 4.8261 12.2068

1.0000 1.0000 1.0000

1.0000 1.0000 1.0000

1.0000 1.0000 1.0000

p-value

5. Real data example We shall analyze Danish fire insurance data set available at http://www.ma.hw.ac.uk/~mcneil/data.html. This data were tested for some forms of positive dependence in Denuit and Scaillet (2004) and Gijbels and Sznajder (2013a,b). This is a multivariate data set. For its thorough description see Gijbels and Sznajder (2013b). The set collects three types of claims referring to losses to buildings (B), their contents (C) and the profit (P) they generated. There are 1502, 529, and 604 strictly positive observations for the respective pairs of (B, C), (B, P), and (C, P).

1.0000 1.0000 1.0000

In Fig. 3 we graphically present the data and related standardized empirical correlations L(r , s)’s. Again, for readability, we restrict attention to d = 15 and r , s = 1, . . . , 15. The displays of the correlations show that for (B, C) and (C,P) the situation is obvious. This is also reflected by p-values of all the tests already considered in Section 4; see Table 1. As before, we took d = 31 in Md . For (B, P) pair there are a few L(r , s)’s well below −3, what invalidates H0 to a large extent. This is reflected by p-values of both AD∗ and Md . Since there are not many strong negative correlations, AD and C v M have problems with detecting them. For illustration, we have also included into Table 1 the related results for classical Pearson’s,

Fig. 3. Scatter plots and empirical correlation coefficients for the Danish fire insurance data. First row: scatter plots of (Ri /(n + 1), Si /(n + 1)), i = 1, . . . , n. Second row: heat maps of the standardized empirical correlation coefficients L(r , s); r , s = 1, . . . , 15.

46

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

Fig. 4. Scatter plots and empirical correlation coefficients for the Danish fire insurance data; positive claims in n = 517 cases. First row: scatter plots of (Ri /(n + 1), Si /(n + 1)), i = 1, . . . , n. Second row: heat maps of the standardized empirical correlation coefficients L(r , s); r , s = 1, . . . , 15.

Spearman’s, Kendall’s, and empirical Blomqvist’s measures of correlation. These statistics are denoted there by P , S , K , and B , respectively, where B = L(1, 1). It is seen that their outcomes are not adequate for (B, C) and (B, P) pairs. In Gijbels and Sznajder (2013a) the analysis of these data was restricted to positive claims in all three variables. In consequence, the three samples were reduced to 517 observations. Such a reduction does not change essentially the inference. For illustration, we provide in Fig. 4 the scatter plots and the standardized correlations. The related p-values are similar as before and we skip their presentation. 6. Conclusions In this paper, we investigate the problem of assessment and validation of positive quadrant dependence. Our approach relies on reparametrization of the testing problem in terms of infinite matrix of some specific correlation coefficients. The estimated correlation coefficients serve to visualize a dependence structure of data at hand and to construct the two new tests for PQD. The test statistics are flexible and work well in detecting many dependence models important in present-day applications. Acknowledgments We thank to the anonymous reviewer and J. Mielniczuk for useful remarks and D. Sznajder for sharing a code. Our work was supported by grant N N201 608440 from the National Science Centre, Poland. Calculations have been carried out in Wrocław Centre for Networking and Supercomputing (http://www.wcss.wroc.pl) under Grant No. 199. Appendix Derivation of (2) √ Since the relation ρr ,s = {C (pr , ps ) − pr ps }/ pr ps (1 − pr )(1 − ps ) plays a key role in this paper, we provide here an elementary justification of it. A similar argument applies when deriving (6).

We have lj (z ) = αj IAj (z ) + βj IBj (z ), where Aj = [0, pj ), Bj =

  [pj , 1], αj = − (1 − pj )/pj , βj = pj /(1 − pj ). Set K1 = Br × Bs , K2 = Ar × Bs , K3 = Ar × As , K4 = Br × As , ω1 = βr βs , ω2 = αr βs , ω3 = αr αs , ω4 = βr αs . With these notations  1 1  4  dC . ρ r ,s = lr (u)ls (v)dC (u, v) = ωi 0

0

i=1

Ki

Since C is uniformly continuous bivariate distribution function satisfying C (u, 0) = C (0, v) = 0, C (u, 1) = u, C (1, v) = v, u, v ∈ [0, 1], the above yields ρr ,s = ω1 + [ω1 − ω2 + ω3 − ω4 ]C (pr , ps ) + pr [−ω1 + ω2 ] + ps [−ω1 + ω4 ]. Noting that, by (18), ω1 − ω2 + ω3 − ω4 = w(r , s), −ω1 + ω2 = −ps w(r , s), −ω1 + ω4 = −pr w(r , s), the conclusion follows.  Information on generation of E1 –E19 We provide a concise description of the alternatives, in a form useful for their generation. E1 : Gauss(θ ), θ ∈ [−1, 1], stands for bivariate normal distribution with 0 means, unit variances and correlation θ . Standard algorithm is taken from Devroye (1986), pp. 566–567. E2 : FGM (θ ), θ ∈ [−1, 1], is Farlie–Gumbel–Morgenstern distribution, see Exercise 3.23, p. 87, Nelsen (2006). E3 : Frank(θ ), θ ∈ R, denotes Frank family generated according to Genest (1987), p. 554. E4 : Clayton(θ ), θ ∈ [−1, ∞), see Exercise 4.16, p. 134, Nelsen (2006). E5 , E7 : mix(θ , γ ), θ ∈ [−1, 1] and |γ | ≤ 1/θ 2 , is mixture copula introduced in Gijbels et al. (2010), p. 567, formula (8). mix(·, ·) is a reparametrized variant of the known Fréchet distribution; cf. Nelsen (2006), pp. 14–15, Exercise 2.4. E6 : Mardia(θ ), θ ∈ [−1, 1], stands for Mardia copula; cf. Nelsen (2006), p. 15, formula (2.2.9) or Gijbels et al. (2010), p. 567. E8 : mixFrank(θ ), θ ∈ R, denotes mixture of Frank copulae: (0.5)Frank(θ ) + (0.5)Frank(−θ ). In the next three cases original uniform representation of (U , V ) is replaced by (U , 1 − V ). E9 : Gumbel(θ ), θ ∈ [1, ∞). We used the following packages implemented in CRAN : copula, mnormt, mvtnorm, scatterplot3d, sn.

T. Ledwina, G. Wyłupek / Insurance: Mathematics and Economics 56 (2014) 38–47

E10 : MO(λ1 , λ2 , λ12 ), λ1 > 0, λ2 > 0, λ12 > 0, denotes uniform representation of Marshall–Olkin distribution, see Nelsen (2006), p. 52. E11 : MS (a, b), a ∈ (0, 1) and b ∈ (0, 1), cf. Example 1.9, p. 25, Mai and Scherer (2012). Next four examples are patterned after Thas and Ottoy (2004). Each generator returns (X , Y ), where X is random explanatory variable, Y = fi (X )+ε, i = 1, . . . , 4, while X and ε are independent standard normal. E12 : qrm(β), β ∈ R, is quadratic regression model Y = β X 2 + ε . E13 : crm(β), β ∈ R+ , is cubic regression model Y = −β X 3 + ε . E14 : horm(β), β ∈ R, is higher order regression model Y = β X 4 + ε. E15 : srm(β), β ∈ R, stands for sinusoidal regression model Y = β sin(2π X ) + ε . E16 : stable(a1 , b1 , a2 , b2 , a3 , b3 , ∆), ai ∈ (0, 2], bi ∈ [−1, 1], i = 1, 2, 3, ∆ ∈ R, denotes bivariate stable distribution with parameters (a1 , b1 , a2 , b2 , a3 , b3 , ∆), defined on p. 293 of Kallenberg and Ledwina (1999). E17 : tν –symmetric, ν > 0, denotes bivariate t-Student distribution with ν > 0 degrees of freedom, see Demarta and McNeil (2005), p. 112. E18 : sub–Gauss(ρ, a), ρ ∈ [−1, 1] and a ∈ (0, 2), stands for sub-Gaussian distribution described on p. 290 of Kallenberg and Ledwina (1999). E19 : tν –skew ed(ρ, γ1 , γ2 ), ν > 0, ρ ∈ [−1, 1], γ1 ∈ R, and γ2 ∈ R, denotes skewed t-copula; cf. the formula (20), p. 119, Demarta and McNeil (2005) with g (w) = w and γ = (γ1 , γ2 ). References Albers, W., 1999. Stop-loss premiums under dependence. Insurance Math. Econom. 24, 173–185. Behnen, K., 1971. Asymptotic optimality and ARE of certain rank-order tests under contiguity. Ann. Math. Stat. 42, 325–329. Behnen, K., Neuhaus, G., 1989. Rank Tests with Estimated Scores and their Application. B.G. Teubner, Stuttgart. Blomqvist, N., 1950. On a measure of dependence between two random variables. Ann. Math. Stat. 21, 593–600. Cont, R., 2001. Empirical properties of asset returns: stylized facts and statistical issues. Quant. Finance 1, 223–236. Demarta, S., McNeil, A.J., 2005. The t copula and related copulas. Internat. Statist. Rev. 73, 111–129. Denuit, M., Scaillet, O., 2004. Nonparametric tests for positive quadrant dependence. J. Financ. Econom. 2, 422–450. Devroye, L., 1986. Non-Uniform Random Variate Generation. Springer-Verlag, New York.

47

Dhaene, J., Denuit, M., Vanduffel, S., 2009. Correlation order, merging and diversification. Insurance Math. Econom. 45, 325–332. Einmahl, J.H.J., 1996. Extension to higher dimensions of the Jaeschke–Eicker result on the standardized empirical process. Comm. Statist. Theory Methods 25, 813–822. Embrechts, P., McNeil, A., Straumann, D., 2002. Correlation and dependency in risk management: properties and pitfalls. In: Dempster, M., Moffatt, H. (Eds.), Risk Management: Value at Risk and Beyond. Cambridge University Press, Cambridge, pp. 176–223. Fermanian, J.-D., Radulović, D., Wegkamp, M., 2004. Weak convergence of empirical copula process. Bernoulli 10, 847–860. Genest, C., 1987. Frank’s family of bivariate distributions. Biometrika 74, 549–555. Gijbels, I., Omelka, M., Sznajder, D., 2010. Positive quadrant dependence tests for copulas. Canad. J. Statist. 38, 555–581. Gijbels, I., Sznajder, D., 2013a. Positive quadrant dependence testing and constrained copula estimation. Canad. J. Statist. 41, 36–64. Gijbels, I., Sznajder, D., 2013b. Testing tail monotonicity by constrained copula estimation. Insurance Math. Econom. 52, 338–351. Groeneboom, P., Lepage, Y., Ruymgaart, F.H., 1976. Rank tests for independence with best strong exact Bahadur slope. Z. Wahrscheinlichkeitstheor. Verwandte Geb. 36, 119–127. Janic-Wróblewska, A., Kallenberg, W.C.M., Ledwina, T., 2004. Detecting positive quadrant dependence and positive function dependence. Insurance Math. Econom. 34, 467–487. Kallenberg, W.C.M., 2009. Estimating copula densities, using model selection techniques. Insurance Math. Econom. 45, 209–223. Kallenberg, W.C.M., Ledwina, T., 1999. Data-driven rank tests for independence. J. Amer. Statist. Assoc. 94, 285–301. Kimeldorf, G., Sampson, A.R., 1978. Monotone dependence. Ann. Statist. 6, 895–903. Ledwina, T., Wyłupek, G., 2012a. Nonparametric tests for stochastic ordering. TEST 21, 730–756. Ledwina, T., Wyłupek, G., 2012b. Two-sample test against one-sided alternatives. Scand. J. Stat. 39, 358–381. Ledwina, T., Wyłupek, G., 2013. Tests for first-order stochastic dominance. Preprint 746, Institute of Mathematics, Polish Academy of Sciences. http://www.impan.pl/Preprints/p746.pdf. Lehmann, E.L., 1966. Some concepts of dependence. Ann. Math. Stat. 37, 1137–1153. Mai, J-F., Scherer, M., 2012. Simulating Copulas. Stochastic Models, Sampling Algorithms, and Applications. Imperial College Press, Singapore. Nelsen, R.B., 2006. An Introduction to Copulas. Springer, New York. Ruymgaart, F.H., 1978. Asymptotic Theory of Rank Tests for Independence. Mathematical Centre Tracts 43, Amsterdam. http://oai.cwi.nl/oai/asset/19566/ 19566A.pdf. Scaillet, O., 2005. A Kolmogorov–Smirnov type test for positive quadrant dependence. Canad. J. Statist. 33, 415–427. Schriever, B.F., 1985. Order Dependence. Mathematical Centre Tracts 20, Amsterdam. http://oai.cwi.nl/oai/asset/12728/12728A.pdf. Schriever, B.F., 1987. An ordering for positive dependence. Ann. Stat. 15, 1208–1214. Shaked, M., Shanthikumar, J.G., 2007. Stochastic Orders. Springer, New York. Tchen, A.H., 1980. Inequalities for distributions with given marginals. Ann. Probab. 8, 814–827. Thas, O., Ottoy, J.P., 2004. A nonparametric test for independence based on sample space partitions. Comm. Stat. Simulation Comput. 33, 711–728. Yanagimoto, T., Okamoto, M., 1969. Partial orderings of permutations and monotonicity of a rank correlation statistic. Ann. Inst. Statist. Math. 21, 489–506.