Mahzad Kalantari
email: mahzad.kalantari@ensg.eu

Franck Jung

PRECISE, AUTOMATIC AND FAST METHOD FOR VANISHING POINT DETECTION

Keywords: vanishing point detection, circles detection, Thales theorem, uncertainty propagation

A new approach for vanishing point detection is described. This method is based on the theorem of Thales. The main contribution of this paper is the automatic and simultaneous detection of all vanishing points of the image that consists in detecting circles in a complex cloud of points. This extraction is performed without any prior knowledge of the internal calibration. An analysis of the error propagation is done in order to give quantitative elements on the precision of the detected vanishing points.

INTRODUCTION

AUTOMATIC detection of vanishing point constitutes a problem widely studied in computer vision for key applications ranging from robotics to camera calibration, or automatic recovery of image orientation. The 3D kind of information issued from the vanishing point localisation inside the image can be used in almost any natural image or video description. The human vision, or photographic images as well, present a characteristic conical geometry, transforming real world parallel lines into bundles of lines, in the recorded image, that intersect on vanishing points. These points are typical features of images from manmade objects, especially buildings, in which almost all the visible or detectable lines correspond to strictly horizontal or vertical elements. The goal of this paper is the presentation of a new algorithm allowing the automatic detection of all vanishing points generated by the segments in a 2D image space. The method presented is simple and robust. Its precision is provided by computing the standard deviation on the vanishing point localization. The model used is based on basic geometric features, working directly in the image space. Section "state of art" reviews the classical approaches. All methods start with the segment detection in the image, an aspect which is not investigated in this paper. The segment classifier uses either a 3D geometric domain based on the Gaussian sphere (or equivalently the Hough transform) or a 2D description from more recent works. Section "Algorithm and implementation" describes the various steps of the algorithm based on Thales theorem, and the variance analysis. Section "Assessments and Results" presents experimental results assessing the method in terms of precision and reliability of detection, and a comparison with a method that works on the unit Gaussian sphere. STATE OF ART [START_REF] Barnard | Interpreting perspective images[END_REF] introduced the most popular algorithm for the detection of vanishing points based on the construction of the Gaussian sphere. Its principle is to embed all vanishing lines of a 2D image inside a 3D projective space to explain the under-relationships between the scene and the projection point of view of a 2D photo. Each vanishing line intersects with the Gaussian sphere. The advantage of this method is to bring back all vanishing points in a finite space. Magee and Aggarwal [START_REF] Magee | Determining vanishing points from perspective images[END_REF] completed the approach by accumulating the projections of the intersection of the image segments on the Gaussian sphere leading to very heavy computations (in the 3D space) and very precise results. Lutton et al. (Lutton 1998) proposed a new adaptation of the Hough transform for the detection of vanishing points and sometimes later Tuytelaars [START_REF] Tuytelaars | The cascaded Hough transform as an aid in aerial image interpretation[END_REF] introduced an interactive method also based on the same finite space of accumulation. Shufelt [START_REF] Shufelt | Performance evaluation and analysis of vanishing point detection[END_REF] tries to find vanishing point on the oblique aerial images using Gaussian sphere.

The second class of methods works directly on the image space. Quan and Mohr (Quan,1989) and Den Heuvel (Den Heuvel, 1998) introduced a detection method based on geometric constraints. In 2003, Almansa [START_REF] Almansa | Vanishing points detection without any a priori information[END_REF] developed a new method of vanishing point detection requiring no a priori information, but using complex probabilistic models. Brauer and Voss (Brauer, 2000) detected vanishing points in highly noisy images. F. Schaffalitzky and A. Zisserman [START_REF] Schaffalitzky | Planar grouping for automatic vanishing lines and points[END_REF] extracted vanishing lines and points by planar grouping.

One can also mention Rother (Rother, 2000) who proposed a very simple method based on the intersection of the segments two by two followed by a majority vote, simple to implement but very heavy in terms of calculation time.

Our work is also based on a simple 2D geometry of the image space. It aims at a fully automatic detection of vanishing points for city buildings and street scenes. The modelling for the 3D scene is not required. Nonetheless, the only a priori requirement is to have enough segments in the image that are automatically detected. These segments correspond to lines, in urban images, almost totally relative to horizontal or vertical features.

ALGORITHM AND IMPLEMENTATION

The different steps of the algorithm are presented in the Fig. 1.

Geometry of the method.

Let S 1 , S 2 , S 3 , … be parallel lines in the real space, their images after a conical projection are converging lines that join on the point (P). This point is conventionally known as the vanishing point. Now a new bundle of lines with O as center (the origin of the image) is built so that, one by one, the lines of this bundle join with a right angle the lines S 1 , S 2 , S 3 , … on the points H 1 , H 2 , H 3 , …. According to the theorem of Thales, H 1 , H 2 , H 3 , … and the origin O define a circle (Fig. 2). In order to determine the point P (the vanishing point) it is enough to determine the parameters of the circle, which is a relatively simple operation. Therefore there are as many circles on the image as vanishing points, and thus to look for the vanishing points that correspond to the intersection of the lines carrying the segments in the images is equivalent to detect these circles (Fig. 3). The exploitation of the theorem of Thales for the research of the vanishing points has already been proposed by Brauer and Voss (Brauer, 2000), but without treating the problem put by the extraction of the mixed different circles, and without any variance analysis. The present work has therefore been completed in this sense. As seen on Fig. 4, the origin O has been chosen in the top left corner of the image, the x-axis being horizontal and the y-axis on the left vertical, the reference formed is not therefore of direct orientation, but every pixel receives positive coordinates. The angle θ is counted from the x-axis. An automatic detection of the segments of the image is performed [START_REF] Deriche | Using Canny's criteria to derive an optimal edge detector recursively implemented[END_REF]. At this step, the covariance matrix of each segment is provided [START_REF] Deriche | From Noisy Edges Points to 3D Reconstruction of a Scene : A Robust Approach and Its Uncertainty Analysis[END_REF].

Among the main results of this step, there are: 1. The extraction of all segments of the image 2. The variance-covariance matrix for each segment. The polar equation of the line, support of each detected segment, is the following : sin cos x y

ρ θ θ = - + (1)
The point H i is defined by : .sin .cos

i i i i i H ρ θ ρ θ  -  =     (2)
It is important to be able to propagate the variance-covariance matrix of each segment on its corresponding point H i . Using the general law of the propagation of variance for non linear functions, and using the Taylor theorem:

, i T i i H J J θ ρ × × Σ = Σ (3)
Where J is the Jacobian matrix that contains the derivatives of H i according to ρ i and θ i . The variance-covariance matrix i H Σ on every point H i is then: 2 2 cos sin cos sin sin cos sin cos

i i i i i i i i i i i i T i i i ρ θ θ σ σ - - - - - -         = × ×              Σ , so that 2 2 , , , , , ,
H i x i y i x H i x i y i y H H i H H H σ σ σ σ σ σ     =       Σ (4)
One can classically illustrate the variance covariance matrixes of the points H i by their error ellipses. The Fig. 6 represents the error ellipses obtained from the segments found in the Fig. 5 and the corresponding points H i .

Circles extraction.

There are two distinct problems: 1 -Which model to choose for the best circle fitting a cloud of points? 2 -How to extract different circles in mixed clouds of points?

Choice of the model of the best circle fitting a cloud of points.

The different ways of modelling a circle have been reviewed; it is not useful to enter into details, as numerous references for this topic exist [START_REF] Gander | Least squares fitting of circles and ellipses[END_REF]. A circle being defined by three points, different obvious ways exist to make a circle pass by a set of more than three points, either by minimizing the algebraic distance between the circle and the different points, and either by minimizing the geometric distance. Of course, each method has its own advantages and drawbacks.

Minimization of the algebraic distance.

Advantage of this method: it is a linear system, with a reduced calculation complexity, and therefore a reduced calculation time. Drawback: this minimization doesn't have any geometrical sense, and besides the evaluation of the parameters of the circle is not very precise.

Minimization of the geometric distance.

The advantage of this minimization is that it is more precise than the previous one [START_REF] Gander | Least squares fitting of circles and ellipses[END_REF]. Besides, it possesses a geometrical sense. On the other hand, its inconvenience is that it implies the resolution of a non linear system, which therefore requires an approached solution. For it, we have used the previously described algebraic minimization. Besides, it is an iterative calculation that requires more computation power.

According to the advantages and drawbacks of each method, thereafter one or the other will be used.

Extraction of the circles in several mixed clouds of points.

The core of the algorithmic implementation of the method and its proper working relies entirely on the efficiency of this step. A method of extraction of the different circles has been inspired from the RanSac method [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. Some modifications had to be brought in order to adapt it to the present situation. The different steps of the algorithm are the following:

1 -Choice of two points H i at random from the whole set of available points H, since two points are sufficient to get a circle passing by the origin.

2 -Calculation of the circle passing by the sampled points and the origin, the mathematical model chosen here is the minimization of the algebraic distance.

3 -Search for all others points H i that are susceptible to contribute to this circle. For that it is necessary to have first fixed a capture threshold, that is to say a strip around this circle: if a point H i fall in this strip of width ε, it is kept for the definition of the circle, and otherwise it is discarded. The way of calculation of this threshold level, from the variance of the points H i , is clarified farther.

4 -Identification of the number of selected points, re-iteration of the steps 1, 2, 3, and finally conservation during these successive tests of the one that captured most points. This operation is reiterated many times in order to insure that the number of points captured is really maximal (here, 1000 times).

5 -At the end of this process, one keeps this set, and the corresponding captured points are then withdrawn from the initial set. Thus for the selected points, one looks for the best circle by a least squares adjustment, at that time with the help of the geometric minimization, with for approached value the parameters obtained by the calculation of the algebraic minimization. This calculation is led taking in account the uncertainties bound to every point H i .

6 -The extraction of the circles is stopped when it remains only, for example, 5 of the initially detected points H i .

So an efficient extraction of the circles relies entirely on the good definition of the value of the width of the capture strip ε. If this value is too small, a considerable number of circles without really different physical senses will be found, and on the other hand if this value is too large, one will make contribute to a given circle, and therefore to the corresponding vanishing point, points that actually should be associated to another circle.

It proves therefore to be fundamental to be able to define the capture threshold in an automatic way for a given image, and it would probably not be very adapted to use the same threshold for images with for example different resolutions. The use of a manual method to define it would cause a heavy loss of time and make lose a major interest of the method. In the optics of a complete automation of the process we did therefore try to extract parameters independent ǫ ǫ from the format of the image, the parameter ε being bound obviously to the precision of the determination of the points H i . In the Fig. 6 it can be seen that the error ellipses form a headband around the points H i , and therefore the median value of the main axes of the ellipses (to 39 % of confidence) can be chosen as a satisfactory threshold level for this capture using RanSaC. Therefore (i being the indication designating every point H i ), ε is the median of the semi-major axes of all error ellipses :

2 2 2 2 2 2 2 2 , , , , , , , , 1 2 ( ( ) 4 (
) ) Calculation of the best origin.

H H H H H H H H i x i y i x i y i x i y i x i y i median σ σ σ σ σ σ σ ε + + + - - = (5)
One of the key points in the algorithm is the choice of the origin. In a first time, because of simplicity, the origin can be chosen on the top left, i. e. on the origin of the image. The choice of the origin has also been discussed in the paper of (Brauer, 2000). However it was only for a given circle. The problem here it is that several circles exist and that it is necessary to choose an unique origin, optimal for all circles. Thus it is suggested to take the origin on the barycentre of the circles that have been detected when the origin is on the top left. It requires to proceed in two steps: In a first time, to detect the circles with the origin in the top left. Then in a second time, to calculate the barycentre of the circles detected and re-build the circles in this configuration. The major advantage is that the distances between the origin and the other circles are equal, so as to get the best geometric intersection between the circles. For example one may see that any fixed arbitrary origin, e. g. the middle of the image (Fig. 8), may be inappropriate in some specific situations, as it creates strong geometric difficulties to separate the circles. Once obtained the parameters of the circle, i. e. its center coordinates (X c , Y c ) since the circle passes by the origin, the coordinates of the P vanishing point are calculated very easily because P is on the OP diameter (Fig. 2); these coordinates are the following : 2

2 p c p c X X Y Y = = (6)
The segments detection in the images, with the help of operators such as Canny-Dériche one, is never reliable to 100%. This does not make any problem as long as one knows how to model this rate of confidence. Besides, the final goal is to be able to calculate the impact of the imprecision of the segments detected on the localization of the vanishing point. This localization will therefore be expressed accompanied of its uncertainty matrix. In this section we demonstrate how to model its impact on the vanishing point, using the uncertainty on the segment.

If we suppose that the error distribution on the segments follows a Gaussian law, the different steps of the calculation are the following [START_REF] Helmert | Die Ausgleichungsrechnung nach der Methode der kleinsten Quadrate. Mit Anwendungen auf die Geodäsie[END_REF], [START_REF] Cooper | Control Surveys in Civil Engineering[END_REF]:

1 -Propagation of the uncertainty of the segment on the coordinates of the corresponding point H i .

2 -Choice of a mathematical model and definition of the observations and unknowns. The chosen mathematical model is the minimization of the geometric distance.

(

)

0 ) , ( = - - = c H c H c X X X X X f (7) H X
is the vector of the points H i , in other words the vector of the observations. In the equation ( 7) the model is defined for exact values of the observations, obviously fictional. For the differentiation in the following, the values c X and H X will be defined as the exact theoretical values.

3 -To get the exact values of the observations and unknowns, the theorem of Taylor is applied:

) .( ) ( ) .( ) ( ) , ( ) , ( 0 0 0 0 0 0 H H H c c c H c H c X X X f X X X f X X f X X f - ∂ ∂ + - ∂ ∂ + = (8)
And if we note

0 0 0 0 0 0 ( ) , B ( ) , ( , ), ( ), (
)

c H c c H H c H f f A b f X X x X X v X X X X ∂ ∂ = = = = - = - ∂ ∂ (9)
the formulation may be simplified in :

b = A x + B v (10)
In other words, A is the Jacobian matrix related to the unknowns, i. e. the coordinates of the center of the circle, B is the Jacobian matrix related to the observations, i. e. the coordinates of the points H i . However, as we have the variance covariance matrix of the points H, we inject this information in the equation ( 9). Thus the weight matrix for the observations is created, named W. For N observations, the weight matrix has a dimension of N x 2 N. If one supposes that the points H i are statistically independent from each other, the matrix will have the following configuration :

      =           Σ Σ Σ (11) i H
Σ is the variance covariance matrix of the points H i provided by the segment detector, matrix that was previously computed (equation 4).

The resolution of this system is made by the minimisation of

T i i i v W v × ×

∑

where v i is the vector of observations. The resolution of the system uses the Lagrange multiplier in an iterative way. At each iteration the vectors x and v are calculated. The convergence is considered as being obtained when the resulting values don't evolve anymore. The center of the circle after compensation is given by : ( )

1 1 1 1 1 ( ) ( ) T T T T c X A B W B A A B W B b - ∧ - - - - = (12)
The variance covariance matrix of the coordinates of the center of the circle is then equal to:

( )

1 1 1 ( ) T T c X A B W B A - - - = Σ (13)
And the uncertainty on the vanishing point is given by:

2 2 , 2 , 2 0 2 0 0 2 0 2 c c c c c c x x y p x y y X σ σ σ σ         =                × × Σ (14)
The equation 14 allows therefore to calculate the uncertainty on the vanishing point P.

ASSESSMENTS AND RESULTS

Some examples obtained with the method presented here are shown in Fig. 9. A visual validation can be obtained by automatically ortho facade including the vanishing point corresponding to the vertical and a vanishing point corresponding to an horizontal direction, proving that the automatic identification of the three vanishing points was correct (cf Fig. 10). This additional step can be required to automatically recognize the building facades. A visual validation can be obtained by automatically ortho-rectifying every facade including the vanishing point corresponding to the vertical and a vanishing horizontal direction, proving that the automatic identification of the three vanishing points was correct (cf Fig. 10). This additional n be required to automatically recognize the building facades. The calculation part of the vanishing points without the error calculation may be performed in real time (less than 0•04 second with a processor Xeon Intel, 1•60 Ghz and 2 Go of Ram, image size 3008 x 2000). Therefore for real time applications, where in a first time it is not necessary to have the maximal precision on the vanishing point, the algorithm works properly. The error calculation part varies, with this not very powerful computer, between 0•5 to 1 minute, according to the number of segments of the image, but this part would be obviously skipped for real time, where the algorithm based on the circles works quite correctly in terms of efficiency and speed. The essential of the calculation time is due to the detection of the segments, as the present algorithm of assignment is comparatively extremely fast. For example for 5600 segments detected on an image of 3008 x 2000 pixels, the algorithm (without the detection of the segments) requires 0•6 seconds without any optimization in the program.

Experimental validation of the variance analysis of the extracted vanishing points.

In order to validate the calculations bound to the extraction of the circles and the analysis of variance of the vanishing points, this method was applied on sets of theoretical segments with known vanishing points coordinates and bundles of segments and lines as shown in Fig. 13. Then a Gaussian noise with a known standard deviation was applied at the segments extremities. The robustness of the algorithm was tested on the two images (cf. Fig. 13) with different segment lengths so as to estimate the effect of the length on the vanishing point accuracy. The set of operations that intervenes downstream to the process of segments extraction was then validated. Comparison with other methods of detection of vanishing points.

It is important to be able to confront the method especially with other methods of detection of vanishing points working on the Gauss sphere, as these methods are reputed as quite efficient for the extraction of the vanishing points close to the infinity. In the examples here under, it is shown that with the method of the circles, the determination of such vanishing points is very satisfactory, and fully comparable to the methods that work in the space of the Gauss sphere.

A synthetic image was chosen with an optical axis perpendicular to the façade, so as to get vanishing points at infinity, which is considered as a quite unfavourable situation for a method working in the plane. The detection using the circles is compared with an algorithm working on the Gauss sphere [START_REF] Kalantari | Robust and Automatic Vanishing Points Detection with their Uncertainties from a Single Uncalibrated Image, by Planes Extraction on the Unit Sphere[END_REF]. As it can be seen on the results, the cases where the vanishing point is at infinite create no problem to classify the different families of segments. Here is another comparison on a complex real case, as the image contains 4 vanishing points, one of them being nearly at infinity (Fig. 16).

Figure 16. Comparison of the method of the circles with a method based on the Gauss sphere on a complex real image

Assessment of the impact of the distance of the vanishing point from the image.

Assessments with the help of synthesized segments have been done, so as to demonstrate the efficiency of the method even when the vanishing point is extremely far. The size of the image chosen is 3008 x 2000 pixels, and the focal distance is 3000 pixels. For each image, about one hundred segments, around 200 pixels long, are synthesized. The assessments take in account two levels of Gaussian noise, one with a standard deviation, judged as realistic, of 0.2 pixel (for the points being the extremities of each segment), and the other with a 1 pixel standard deviation, in order to value the impact of situations with very high noise. The mean values of 30 consecutive results are presented so as to facilitate the reading of the results (Fig. 17 and Fig. 18). The reference system is the same as in the Fig. 4, with the origin of the system in the top left of the image. The x coordinate of the vanishing point is fixed to 1504 pixels. For each new iteration, 1000 pixels are added to the y coordinate of the artificial vanishing point, that varies thus from 0 to 1 000 000 pixels. One can consider that for an image of 3008 x 2000 pixels, 1 000 000 pixels is an infinite value. The angular error is defined as the angle difference between the theoretical direction of the vanishing point and the direction of the valued vanishing point. In this simulation the method presented in this paper has been confronted to the method of the Gauss sphere [START_REF] Kalantari | Robust and Automatic Vanishing Points Detection with their Uncertainties from a Single Uncalibrated Image, by Planes Extraction on the Unit Sphere[END_REF].

Figure 17. Results of the assessment, using the method of the circles Figure 18. Results of the assessment, using the method of the Gauss sphere As one can see in Fig. 17 and Fig. 18, the important result is that the more or less distant localization of the vanishing point doesn't have any impact on the precision of its direction so determined, whatever the used method. In other words the method of the circles, as well as the method using the Gauss sphere (one may note that it is clearly less sensitive to the noise of the segments), operate in a satisfactory way whatever the localization of the vanishing point. The major factor that influences the determination of the localization of the vanishing point is, as anticipated, the noise on the segments. In the absence of noise the angular error is always null, showing the stability of the calculation method. One notes without surprise that for a noise of 0.2 pixel, the error found for the direction of the vanishing point is about 5 times weaker than for a noise of 1 pixel.

For this reason it is important to have an evaluation of the localization error of the vanishing point, with the help of the error propagation presented in the section xx. This propagation allows us therefore to know accurately the impact of the segments noise on the localization of the vanishing point.

To conclude this part, it is important to note that the method of the circles, even working in the image space, doesn't meet the limitations found in methods based on the intersection of the segments on an image such as Rother's one.

CONCLUSION

For images acquired in urban environment, the presented algorithm allows to recover vanishing points entirely automatically, and to simultaneously determine their variance. One major characteristic of this algorithm is its time efficiency when compared to a Gauss sphere method. Its efficiency is comparable, even when the vanishing points are at infinity. Besides, its accuracy can still be measured. As a result, this gives a simple way to calculate the relative orientation of these images as well as a partial ortho-rectification for a real time process.

ACKNOWLEDGMENT

These works have first been led thanks to the support of the Urban Community of Le Mans, of the General Council of Sarthe and the Regional Council of the Pays de Loire, within the Ecole Supérieure des Geomètres et Topographes (ESGT) of the Conservatoire National des Arts et Métiers.

Fig. 1 .

 1 Fig. 1. The different steps of the algorithm presented in this paper.

Fig. 2 .

 2 Fig. 2. Geometric configuration based on the Thales theorem

Fig. 3 .

 3 Fig. 3. Example of a building, and the set of corresponding points H. One notes that the corresponding circles of the different perspective families go through the origin of the image

Fig. 4 .

 4 Fig. 4. Definition of the image reference system

Fig. 5 .

 5 Fig. 5. Example of segments detection in an urban scene

Fig. 6 .

 6 Fig. 6. : Error ellipses (at 95% confidence level) showing the error determination of the points H, depending of the efficiency of the tool of automatic extraction of the segments (from the image presented in Fig. 5).

Fig. 7 .

 7 Fig.7 . Points Hi kept for a given circle of width ε represented here by the zone between the two dashed curves

Fig. 8 .

 8 Fig. 8. : Examples of points Hi for different configurations of the origin. Using the barycentre (b), the geometric intersection of the circles is optimal. a) origin in the top left. b) origin in the barycentre. c) origin in the middle of the image (unfavourable configuration) (a) (b) (c)

Figure 9 .

 9 Figure 9. Different examples of automatic identification of families of segments corresponding to parallels in the real space, with the superimposition of the points vanishing point) extracted, the segments are presented with the same colour.

Fig. 10 .Figure 9 .

 109 Fig. 10. Based on the automatic extraction of various vanishing points, an automatic orthorectification of two facades is shown, as a simple and efficient visual test

Fig. 10 .

 10 Fig. 10. Based on the automatic extraction of various vanishing points, an automatic orthorectification of two facades is shown, as a simple and efficient visual test (using the image of Fig. 3).

Fig. 11 .

 11 Fig. 11.These 8 images are typical of the images where the automatic detection of vanishing points has been tested, and for which all expected points have been found

Figure 13 :Figure 14 :

 1314 Figure 13: Theoretical segments used for checking the variance analysis of the automatically extracted vanishing points. In (13 a) the length of the segments range from 100 to 150 units (arbitrary units, the image of Fig. 13 measures 3000 x 2000 units). In (13 b) the length ranges from 1500 to 2200 units.

Fig. 15 .

 15 Fig. 15. Comparison of the method of the circles (left) with a method based on the Gauss sphere (right) in a case where the vanishing points are at infinity

Reliability of the detection of the vanishing points.

The previous algorithms were applied to a set of different urban images, in order to identify the limits of the proposed method. To avoid lenient validations, a database has been created and is on free access on our web site http://mahzad.kalantari.free.fr/pdf.htm On Fig. 11, some images are presented where the detection functioned satisfactorily, whereas Fig. 12 shows images where no vanishing points were provided by the algorithm. Here the results are only assessed by visual inspection of the correct number of vanishing points that must be found. These results can be considered as satisfactory.


On the relation between extremal dependence and concomitants

Keywords: Asymptotic Theorems, Concomitants, Copula, (Tail) Dependence, Extremes, Gaussian, Logistic, Maxima, Order Statistics, Pareto, Slowly/Regularly Varying Functions, Tail Equivalence, Weak Convergence. 2010 AMS classification: 60G70, 62E20, 62G30

The study of concomitants has recently met a renewed interest due to its applications in selection procedures. For instance, concomitants are used in rankedset sampling, to achieve efficiency and reduce cost when compared to the simple random sampling. In parallel, the search for new methods to provide a rich description of extremal dependence among multiple time series has rapidly grown, due also to its numerous practical implications and the lack of suitable models to assess it. Here, our aim is to investigate extremal dependence when choosing the concomitants approach. In this study, we show how the extremal dependence of a vector (X, Y) impacts the asymptotic behavior of the maxima over subsets of concomitants. Furthermore, discussing the various conditions and results, we investigate how transformations of the marginal distributions of X and Y influence the degeneracy of the limit.

Let (X i , Y i ), i = 1, . . . , n, be a sequence of independent and identically distributed (i.i.d.) bivariate random vectors with common cumulative distribution (cdf) function F(x, y). We define the concomitants of order statistics as follows. Arrange the X variates in ascending order, thus obtaining the order statistics X (1)

X (2) • • • X (n)
for the X variable. Then, the Y variable corresponding to the i-th order statistic X (i) is called the concomitant of the i-th order statistic, and is denoted by Y [i] . For example, if the maximum of the X variates is X (n) ≡ X 3 , then its concomitant is Y [n] ≡ Y 3 . Concomitants of order statistics are of notable interest in practical applications: It is worth mentioning their utility in problems concerning the estimation of parameters for multivariate datasets affected by type II censoring, in the development of selection procedures (where it is more convenient to measure a characteristic linked to a measure of interest), in the setting of ranked-set sampling schemes. In particular, the latter deals with situations where the measurement of the variable of interest is time-consuming or expensive, but the ranking of a set of observations linked to the target variable is easy to be performed. For such cases, ranked-set sampling is a better alternative with respect to the simple random sampling, both in terms of efficiency and reduced cost, hence this renewed attention on the study of concomitants. We report the seminal works [START_REF] Dell | Ranked set sampling theory with order statistics background[END_REF], [START_REF] Stokes | Characterization of a ranked-set sample with application to estimating distribution functions[END_REF], and more recently [START_REF] Wang | Concomitants of multivariate order statistics with application to judgment poststratification[END_REF], [START_REF] Shao | The dependent wild bootstrap[END_REF], [START_REF] Balci | Modified maximum likelihood estimators using ranked set sampling[END_REF], [START_REF] Wang | Using ranked set sampling with cluster randomized designs for improved inference on treatment effects[END_REF], [START_REF] Wang | Unbalanced ranked set sampling in cluster randomized studies[END_REF], Zamanzade and Wang (2018), [START_REF] Ozturk | Two-stage cluster samples with ranked set sampling designs[END_REF] among others. Our motivation is to investigate extremal dependence when choosing the concomitant approach. Looking at the literature, this topic has been tackled in the 90's. One can mention the seminal paper [START_REF] Nagaraja | Distribution of the Maximum of Concomitants of Selected Order Statistics[END_REF] where, given a number k, they study the distribution of the maximum of the concomitants V 1 = max Y [n] , . . . , Y [n-k+1] . Finite and asymptotic results are given in both cases k a fixed integer (extreme case) and k = [np], p ∈ (0, 1) (quantile case). This leads to the natural question of how much V 1 is close to the maximum Y (n) , a fundamental problem in practical applications, tackled in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF]. In their work, authors propose to link this question to the study of the joint distribution of two maxima of concomitants. More precisely, given an i.i.d. sequence of bivariate random vectors with parent random vector (X, Y), define, for k 1,

V 1 = max Y [n] , . . . , Y [n-k+1] and V 2 = max Y [n-k] , . . . , Y [1] .
(1) [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF] introduced the random variable W k := V 1 /Y (n) and proved that its cdf satisfies

P(W k w) = P(V 1 wV 2 , V 2 > 0) ,
for which they need to compute the joint distribution of (V 1 , V 2 ).

Before briefly reporting their result, let us introduce some notations that will be used throughout the paper.

Let G(z; µ, σ, ξ), for µ, ξ ∈ R, σ > 0, denote a Generalised Extreme Value distribution (GEV). A distribution function F is said to belong to the Maximum Domain of Attraction (MDA) of G, written as F ∈ D(G), if there exist normalizing sequences of constants a n > 0 and b n ∈ R such that F n (a n z + b n ) → G(z), as n → ∞ for any z ∈ Γ G = {z : ξ(z -µ)/σ > 0}. The structure of a MDA is studied in [START_REF] Fisher | Limiting forms of the frequency distribution of the largest or smallest member of a sample[END_REF], [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une serie aleatoire[END_REF], among others, and more recently, with another view, in [START_REF] Leonetti | The maximum domain of attraction of multivariate extreme value distributions is small[END_REF]. Sufficient conditions for a distribution F to belong to the MDA of a certain GEV G are the so called Von Mises conditions: We refer to Resnick (1987, Propositions 1.15 -1.17) for a thorough description and provide a brief recall in Appendix C for the paper to be self-contained.

A measurable function h : R + → R + is regularly varying at infinity with index a (written h ∈ RV a ) if for x > 0

lim t→∞ h(tx) h(t) = x a .
The index a is called the exponent of variation. When a is 0, h is said slowly varying (SV).

The upper-tail dependence coefficient of a random vector (X, Y) with copula C and marginal distributions F X and F Y , respectively, is defined as [START_REF] Coles | Dependence measures for extreme value analyses[END_REF])

λ u := lim α→1 P Y > F -1 Y (α) | X > F -1 X (α) = lim α→1 1 -2α -C(α, α) 1 -α , (2)
F -1 denoting the generalized inverse function of a cdf F. Note that λ u ∈ [0, 1], where λ u = 0 means extremal independence between X and Y.

Let us also introduce the following notation, borrowed from [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF]. Given a vector (X, Y), we define

F 1 (y | x) := P(Y y | X > x) , (3)
F 2 (y | x) := P(Y y | X x) , (4)
F 3 (y | x) := P(Y y | X = x) . ( 5
)
The conditional distributions F 1 , F 2 and F 3 play a fundamental role in computing the distributions of the concomitants of the order statistics as well as the distributions of their transformations.

Theorem 2 in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF]: Suppose F X satisfies one of the Von Mises conditions and, for all x and y, assume

F X ∈ D(G X ), F Y ∈ D(G Y ), and
nP(X > a n x + b n , Y > c n y + d n ) → n→∞ 0, ∀x, y ∈ R, (6)
where the norming constants

(a n ) n , (b n ) n , (c n ) n and (d n ) n are such that F 3 (c n y + d n | a n x + b n ) → 1 as n → ∞. (7)
Further, suppose there exist constants ãn > 0 and bn ∈ R such that

F 1 ( ãn y + bn | a n x + b n ) → n→∞ H 1 (x, y). (8)
Then, the cdf of (V 1 , V 2 ) satisfies

F (V 1 ,V 2 ) ( ãn v 1 + bn , c n v 2 + d n ) → n→∞ H(v 1 )G Y (v 2 ), ∀v 1 , v 2 ∈ R, (9)
where

H(v 1 ) = (k!) -1 H k 1 (x, v 1 )(-log G X (x)) k g X (x) dx.
We begin our analysis by focusing on Condition (7). We first show in Lemma 1 the relation between ( 6) and ( 7), which proof is given in Appendix A.

Lemma 1. Let x, y ∈ R. Then, we have

nP(X > a n x + b n , Y > c n y + d n ) → n→∞ 0 ⇔ F 3 (c n y + d n | a n z + b n ) → n→∞ 1, ∀z x. (10)
Now, we can see that ( 7) can be directly understood in terms of the upper-tail dependence coefficient, namely:

Proposition 1. Consider a bivariate random vector (X, Y) with copula C. Then,

F 3 (c n y + d n | a n x + b n ) → n→∞ 1, ∀x y ⇔ λ u = 0
where the normalizing constants are chosen depending on the MDA to which X and Y belong.

Note that Proposition 1 is then simply another way to express Proposition 5.27 in Resnick (1987, p. 296), which states that

λ u = 0 ⇔ nP(X > a n x + b n , Y > c n y + d n ) → n→∞ 0.
Using Lemma 1 and Proposition 1, we can now rewrite Theorem 2 in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF] in terms of λ u . This reformulation, given in Corollary 1, is helpful since it clarifies the connection between the asymptotic independence of (V 1 , V 2 ) and that of (X, Y).

Corollary 1. Suppose F X satisfies one of the Von Mises conditions and, for all x, y,

F X ∈ D(G X ) and F Y ∈ D(G Y ).
Moreover, assume that the upper tail dependence coefficient between X and Y is 0: λ u = 0. If there exist constants ãn > 0, bn s.t.

F 1 ( ãn y + bn | a n x + b n ) → n→∞ H(x, y), then the cdf F (V 1 ,V 2 ) of (V 1 , V 2 ) satisfies (9).
The assumption of asymptotic independence between X and Y is generally too restrictive in practical applications, especially in this highly interconnected world. For instance, over the past 15 years, we have had to live through a major global financial crisis, then a pandemic, combined with an increase in cyber attacks. The presence of systemic risk, with a strong dependence among extremes, is one characteristic of all those events (e.g. [START_REF] Dacorogna | Living in a stochastic world and managing complex risks[END_REF]). Hence, investigating this (asymptotic) dependence aspect is of primary interest. This is the main objective of the next section.

We aim at studying the limiting distribution of the concomitants (V 1 , V 2 ) when X and Y may exhibit some asymptotic dependence. To do so, we take advantage of the joint-tail model representation given in [START_REF] Ledford | Concomitant tail behaviour for extremes[END_REF], then we state the main result (in Subsection 2.2) and discuss the choice of the normalizing constants to avoid degeneracy in the limit.

. The LT joint-tail model

First, let us recall some definitions related to the slowly varying notion in the bivariate case. It will be needed to define the joint-tail model.

Definition 1. A function L that satisfies lim n→∞ L(nx, ny) L(n, n) = r(x, y), with r(ax, ay) = r(x, y), (11)
for all a > 0 and (x, y) ∈ R 2 + , is said to be bivariate slowly varying (BSV).

It can be shown that there exists a univariate function r * such that

r(x, y) = r * (w), with w = x x + y . ( 12
)
Definition 2. A BSV function L for which the function r * defined in (12) satisfies r * (w) r * (1 -w) is SV at w = 0 and w = 1, is said to be quasi-symmetric.

Let us now present the joint-tail model, named the LT joint-tail model, introduced in [START_REF] Ledford | Concomitant tail behaviour for extremes[END_REF] to study the asymptotic distribution of the concomitant of the max, Y [n] , i.e. the limit distribution of the marginal distribution of V 1 (case k = 1), and to compute P Y

(n) = Y [n] .
Let (X, Y) be a bivariate random vector with unit Fréchet marginals (i.e. F X (x) = e -1/x , x > 0) and joint distribution function F(x, y). Suppose that there exist functions c(t) and ψ(x, y) for which the joint survival function

F(x, y) = 1 -e -1/x -e -1/y + F(x, y) satisfies lim t→∞ t F(tx, ty) c(t) = ψ(x, y) for all (x, y) ∈ R 2 + . ( 13
)
Additionally suppose that there exists γ ∈ R such that F(0, t)/ F(t, 0) is RV γ at infinity and

F(t, 0)/ F(0, t) is RV -γ at infinity. ( 14
)
Then the following holds.

Definition 3 (The LT joint-tail model). Assume (X, Y) have unit Fréchet marginals. Let η ∈ (0, 1] be the coefficient of tail dependence used to determine the decay rate of F -1 (t, t), as t → ∞. Under Assumptions ( 13) and ( 14), the joint survival distribution of (X, Y) is given by

P(X > x, Y > y) = L(x, y)x -α y -β , ∀x, y, (15)
where α, β > 0, α + β = η -1 and the function L(x, y) is a quasi-symmetric BSV function (see Definition 2).

The coefficient η describes the type of limiting dependence between X and Y, while the function L its relative strength given a particular value of η; see Ledford and Tawn (1998) and [START_REF] Heffernan | A directory of coefficients of tail dependence[END_REF] for further comments. As pointed out in de [START_REF] De Haan | Extreme residual dependence for random vectors and processes[END_REF], the η = 1/2 corresponds to the case where X and Y are independent, while the case η > 1/2 indicates a positive association of the extremes of (X, Y). Values of η lower than 1/2 represent cases where the extremes of X and Y are negatively associated. The bounding cases of perfect negative and positive dependence correspond respectively to η → 0 and η = 1 with L(•) = 1. [START_REF] Coles | Dependence measures for extreme value analyses[END_REF] introduce an elementary measure of dependence: Since the survival copula C is given by C(u, v) = 1 -u -v + C(u, v), they consider the following indices to describe the upper-tail dependence:

χ := lim u→1 2 log(1 -u) log C(u, u) -1 = 2η -1 (16)
λ u := 2 -lim u→1 log C(u, u) log u =      c if χ = 1, L(t) → c > 0, as t → ∞, 0 if χ = 1, L(t) → 0, as t → ∞, 0 if χ < 1. (17)
Since χ = 1 if and only if η = 1, we have a link between λ u and η as well. In other words,

• X and Y are asymptotically independent when α + β = 1 and L(n, n) → n→∞ 0, or when α + β 1 (with no condition on the limit as n → ∞ of L(n, n));

• X and Y are asymptotically dependent when α

+ β = 1 and L(n, n) → n→∞ c > 0. Note that, in this case, c = λ u .
This turns out to be useful when considering extensions of Theorem 2 in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF], relaxing Condition (7).

For the sake of simplicity, we will use the notation C ⊥ ⊥ to indicate the asymptotic dependence case, that is α + β = 1 and L(n, n) → n→∞ λ u > 0. Thus, the indicator 1 C ⊥ ⊥ = 1 points out the presence of asymptotic dependence.

.

Main result

Let us state our main result on the joint asymptotic behaviour of maxima of concomitants, allowing for asymptotic dependence of X and Y:

Theorem 1.

Let (X, Y) be a bivariate random vector with unit Fréchet marginals and upper-tail dependence coefficient λ u . Assume that (X, Y) follows the LT joint-tail model given in Definition 3, with survival cdf (15), and that the BSV function L is such that L(nx, ãn y + bn ) → Then, the joint distribution of the concomitants maxima (V 1 , V 2 ) defined in (1) satisfies:

F (V 1 ,V 2 ) ( ãn v 1 + bn , nv 2 ) ---→ n→∞ +∞ 0 H k 1 (v 1 | x) H 2 (v 2 | x) x -k-2 k! e -1/x dx, ∀v 1 , v 2 ∈ R, (19)
where

H 1 (y | x) := 1 -c(x, y) x 1-α y -β and H 2 (y | x)
is defined by the product of the limits given below in (20) and ( 21), respectively, with r(x, y) = λ -1 u c(x, y).

To prove Theorem 1, we need to evaluate the asymptotic behavior of the three conditional distributions F 1 , F 2 and F 3 defined in (3), ( 4) and ( 5), respectively, each adequately transformed to obtain a non-degenerated limit distribution. This is what is presented in the following lemma.

Lemma 2. Assume the LT joint-tail model (15) holds and call r the limit function corresponding to L. Then, as n → ∞,

F n 2 (ny | nx) → exp - 1 y 1 -λ u r(x, y) y x α 1 C ⊥ ⊥ , ( 20
)
F 3 (ny | nx) → 1 + λ u x -α+2 y α-1 ∂ ∂x r(x, y) - α x r(x, y) 1 C ⊥ ⊥ . ( 21
)
Moreover, under Condition (18), we have

F 1 ( ãn y + bn | nx) → n→∞ c(x, y). ( 22
)
Note that taking 1 C ⊥ ⊥ = 0 in ( 20) and ( 21) provides lim n→∞ F n 2 (ny | nx) = e -1/y and lim n→∞ F 3 (ny | nx) = 1, which characterise copulas with asymptotic independence, as considered in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF]. The proofs of Theorem 1 and Lemma 2 are developed in Appendix A. It is worth noticing that Condition (18) is needed to ensure that H 1 (v 1 | x) is a proper cdf for v 1 ∈ (0, +∞). Moreover, our Theorem 1 gives precise formulations on the choice of the normalizing constants ãn and bn . This is in fact a delicate point, reason why we dedicate a comprehensive discussion on it: Explanations relative to Condition (18) are given in the next subsection, while Section 4 addresses the fundamental example of the Gaussian bivariate distribution, for which Condition (18) does not hold.

. Discussion on the choice of ãn and bn

In the following, we carefully deduce which are the only possible asymptotic behaviours of the constants ãn and bn . We highlight the importance of the present subsection, as finding suitable constants which lead to non degenerate limits represents generally an issue in asymptotic theories.

We start by noticing that the limit of F 1 ( ãn y + bn | nx) strongly depends on that of L(nx, ãn y + bn ). Indeed, we can write

F 1 ( ãn y + bn | nx) = 1 - L(nx, ãn y + bn )(nx) -α ( ãn y + bn ) -β 1 -e -1/(nx) ∼ n→∞ 1 - L(nx, ãn y + bn ) n α-1 x α-1 ( ãn y + bn ) β (23) using the approximation 1 -e -1/(nx) ∼ n→∞ (nx) -1 .
In the following, we discuss which sequences ãn and bn should be chosen to have

F 1 ( ãn y + bn | nx) → n→∞ 1, i.e. V 1 not degenerate.
Given a couple (x, y), we study the limit of ( 23) depending on the behavior of L, namely:

(i) Assume L(nx, ãn y + bn ) → 0 as n → ∞. To have F 1 ( ãn y + bn | nx) → 1, we need L(nx, ãn y+ bn ) n α-1 x α-1 ( ãn y+ bn ) β → 0.
But, since we are in the case L(nx, ãn y + bn ) → n→∞ 0, this would mean, either

• n α-1 x α-1 ( ãn y + bn ) β →
n→∞ 0 at the same rate of L, which is impossible because of the assumptions of the model (L ∼ SV), or

• n α-1 x α-1 ( ãn y + bn ) β → n→∞ 0 faster than L, which is impossible because we would have

L(nx, ãn y+ bn ) n α-1 x α-1 ( ãn y+ bn ) β → n→∞ ∞, that denies F 1 is a cdf. Hence case (i) is not possible. (ii) Assume L(nx, ãn y + bn ) → ∞ as n → ∞. So, as n → ∞, the term n α-1 x α-1 ( ãn y + bn ) β tends to infinity faster than L(nx, ãn y + bn ) does, which yields that F 1 ( ãn y + bn | nx) → n→∞ 1.
Hence, this case (ii) has also to be discarded.

(iii) The previous reasoning implies that we can only have, as n → ∞,

L(nx, ãn y + bn ) → n→∞ c(x, y) ∈ (0, ∞).
Which conditions on ãn and bn can be given? Here, we have to distinguish whether c(x, y) depends on y or not.

(A) Suppose c(x, y) depends on y. Then, L(nx, ãn y + bn )

n α-1 x α-1 ( ãn y + bn ) β ∼ n→∞ c(x, y) n α-1 x α-1 ( ãn y + bn )
β , so we have to consider the relative behaviour of ãn and bn .

(a) If ãn = o( bn ), then ãn y + bn ∼ bn , as n → ∞. This implies c(x,y)

n α-1 x α-1 ( ãn y+ bn ) β ∼ n→∞ c(x,y) n α-1 x α-1 bβ n
, so, in order to maintain the dependence on y, we need n α-1 bβ n → n→∞ 0, ∞. This condition obviously depends on the value of α: Call γ = (1 -α)/β, then

• If 0 < α < 1, n α-1 → n→∞ 0, the only possibility is bn = O (n γ )
with γ ∈ (0, 1), since α + β 1.

• If α = 1, then bn would be independent of n; it is then impossible to propose ( ãn , bn ) in such a case.

• If α > 1, then bn = O (n γ ) with γ < 0. However, since we are interested in the dependence in the extremes, we need to consider

nP(X > nx, Y > n γ y) n 1-γ • n γ P(X > n γ x, Y > n γ y) → n→∞ 0,
which implies γ ∈ (0, 1). So, the case α > 1 cannot be considered.

(b) If bn = o ( ãn ) or bn = O ( ãn ), then ãn y + bn ∼ n→∞ ãn y and c(x, y) n α-1 x α-1 ( ãn y + bn ) β ∼ n→∞ 1 n α-1 ãβ n c(x, y) x 1-α .
Hence, we need to consider ãn = O (n γ ) and α ∈ (0, 1).

(B) Suppose c(x, y) does not depend on y. In this case, L(nx, ãn y + bn )

n α-1 x α-1 ( ãn y + bn ) β ∼ n→∞ n 1-α max ãβ n y β , bβ n cx 1-α ,
and since we do not want to loose the dependence on y, we need to assume max ãβ n y β , bβ n = ãβ n y β . Moreover, we need to ask

n 1-α ã-β n → 0 (otherwise F 1 ( ãn y + bn | nx) → n→∞ 1), which implies ãn = O (n γ ) and α ∈ (0, 1).
Therefore, to obtain a non degenerate limit for F 1 ( ãn y + bn | nx), we need to assume α ∈ (0, 1), ãn = O (n γ ) and bn = o ( ãn ) or bn = O ( ãn ). Now that the main result is provided with a justification of the constants involved, let us illustrate the theoretical results on standard examples of the risk literature.

In this section, we consider two examples: the Pareto -Lomax with Survival Clayton copula and the Bivariate Symmetric Logistic Extremal copula. We apply our asymptotic approximations given in Theorem 1 and compare it with the empirical results obtained via simulation. The choice of such examples lies in the fact that the extremal dependence is ruled by parameters, so easy to control. Note also that for the following examples we can achieve the same results by direct computation (with a considerable amount of time and tedious calculations), thus it is possible to make comparisons by means of the mathematical expressions obtained. We simulate 5000 replicates of 100000 bivariate samples. The dimension k of the subsample used to compute V 1 is set to 10 (Appendix B contains simulations relative to other values of k, to have insights of its impact on the asymptotic result).

.

Pareto -Lomax marginals with Survival Clayton copula

Let (X P , Y P ) be a random vector with Survival Clayton(θ) copula and Pareto -Lomax(ν, 1) marginals (the subscript P indicates the Pareto marginals): First, we transform te Pareto marginals into the unit Fréchet ones, to match the assumptions of the model given in Theorem 1. We can write where, for the Pareto distribution, F -1

FX P (x) = (1 + x) -ν , ∀x > 0, FY P (y) = (1 + y) -ν , ∀y > 0, ν > 0, FP (x, y) = P(X P > x, Y P > y) = (1 + x) νθ + (1 + y) νθ -1 -1/θ , θ > 0. (a) Survival Clayton copula (θ = 2) (b) Logistic copula (γ = 1/2)
F(x, y) = F P F -1 X P e -1/x , F -1 Y P e -1/y , (a) Bivariate distribution of (V 1 , V 2 ).
X P (q) = F -1 Y P (q) = (1 -q) -1/ν -1. We obtain F(x, y) = 1 -e -1/x -θ + 1 -e -1/y -θ -1 -1/θ . ( 24
)
Now, observe that the joint survival (24) can be rewritten, for x > x 0 > 0 and y > y 0 > 0, as

F(x, y) ∼ x θ + y θ -1 -1/θ = x θ + y θ -1 -1/θ (xy) 1/2 (xy) -1/2 = L(x, y)x -1/2 y -1/2 , where L(x, y) = x θ + y θ -1 -1/θ (xy) 1/2 (xy) and clearly α = β = 1/2. Note that L(n, n) = (2n θ -1) -1/θ n → 2 -1/θ as n → ∞.
This is coherent with the value of the upper-tail dependence coefficient, since

λ U = lim q→1 P Y > F -1 Y (q) | F -1 X (q) = lim q→1 (1 -q) -1 (2(1 -q) -θ -1) -1/θ = 2 -1/θ .
According to Theorem 1, since 1 -α = β, we have ãn = a n = n. Hence,

L(nx, ny) = n θ x θ + n θ y θ -1 -1/θ n(xy) 1/2 → n→∞ x θ + y θ -1/θ (xy) 1/2 = c(x, y).
Note that, for each a > 0, c(ax, ay) = (ax) θ + (ay) θ -1/θ (a 2 xy) 1/2 = c(x, y), which means that c(x, y) and also r(x, y) := 2 1/θ c(x, y) satisfy ( 11). Moreover it is straightforward to show that c(x, y) satisfies Condition (18) in Theorem 1.

Since it holds

∂ ∂x r(x, y) = 2 1/θ ∂ ∂x c(x, y) = 2 1/θ x -1 c(x, y) 1 2 -1 + y x θ -1
, applying Lemma 2 and Theorem 1 gives that

H 1 (y | x) = 1 -1 + y x θ -1/θ , ( 25a
)
H 2 (y | x) = 1 -1 + y x θ -1/θ-1 e -1 y + 1+( y x ) θ -1/θ . ( 25b
)
The plot in Figure 2a shows the bivariate surface computed by applying Theorem 1 (the integration is performed numerically) and the values of the bivariate empirical cdf computed in the simulated V 1 and V 2 . These points lie on the surface, thus displaying the very good quality of the asymptotic approximation of Theorem 1. It can also be assessed from Figures 2b and2c, showing the boxplot of the absolute value of the errors between the black cloud and the surface in Figure 2a, the asymptotic approximation of the marginal distributions of V 1 and the errors in V 2 , respectively. For what concerns the error (lower plot of Figure 2c), we observe that it is small for low values of v 2 (the maximum error is close to 1%) and converges very quickly to 0. In particular, it is 0% for all the values of v 2 we are considering.

. Unit Fréchet marginals with Logistic Extremal copula

Let (X, Y) be a random vector with unit Fréchet marginals and a bivariate symmetric logistic extremal dependence, i.e.

F(x, y) = e -V(x,y) , where V(x, y) = x -1/γ + y -1/γ γ , x, y > 0, γ ∈ (0, 1).

By considering the bivariate joint tail approximation (15) proposed in [START_REF] Ledford | Modelling dependence within joint tail regions[END_REF] characterised by a general extremal dependence function V, we can write in our case

L(x, y) = (xy) -1/2 (x + y) -(xy) 1/2 V(x, y). ( 26
)
Note that

(a) Bivariate distribution of (V 1 , V 2 ).
L(n, n) = 2 -V(1, 1),
since V is homogeneous of order -1 (i.e. V(nx, ny) = n -1 V(x, y)). Moreover, 1 -α = β, so that ãn = a n = n. Hence, we can write

L(nx, ny) = n -1 (xy) -1/2 n(x + y) -n(xy) 1/2 n -1 V(x, y) = (xy) -1/2 (x + y) -(xy) 1/2 V(x, y) = c(x, y).
Note that, for each a > 0, c(ax, ay) = c(x, y), so that the candidate L in ( 26) is a proper BSV function. Simple algebra helps to prove that c(x, y) satisfies Condition (18) in Theorem 1. Since it holds

∂ ∂x c(x, y) = (xy) -1/2 2 1 - y x -yV(x, y) 1 -2(xV(x, y)) -1/γ ,
by applying Lemma 2 and Theorem 1, we obtain

H 1 (y | x) = x V(x, y) -y -1 , (27a)
H 2 (y | x) = 1 + y x -1 γ γ-1 e -1 x 1+( y x ) -1 γ γ -1 . ( 27b
)
After having obtained the functions H 1 and H 2 , we can now apply a numerical integration to compute the joint cdf (19).

Whilst our focus so far was on obtaining the joint asymptotic cdf of (V 1 , V 2 ) when (X, Y) exhibits asymptotic or extremal dependence, we question if Condition (18) of Theorem 1 is strongly related to the extremal dependence assumption. In other words, what does occur if the BSV function L goes to 0 or to infinity? This is why we return to the case when (X, Y) is Gaussian, as it is well known that (X, Y) is then asymptotically independent. This case has already been treated in the literature, due both to practical interest of the distribution and to its known manageability. Although direct computations in this case become simple, problems arise in applying models involving approximations with SV functions. In fact, even if such approximations are provided, they are then not used to compute for instance the conditional distributions we encountered above (as far as we could observe in the literature). The objective of the present section is to question this apparent paradox and provide tools to overcome this issue, thus bridging the gap between direct computations and those made through asymptotic models, with some attention put on the asymptotic behaviour of the BSV L.

.

Existing literature

Let us start to review the existing literature on the Gaussian case.

(a) [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF]: By direct computations, the authors are able to provide constants ãn , bn such that the rescaled distribution of (V 1 , V 2 ) converges to a product of a standard normal cdf to the power k and a Gumbel cdf. When applying those constants in our case, after the transformation of the Gaussian marginals into unit Fréchet, Condition (18) in Theorem 1 is not satisfied. In fact, we obtain that F 1 degenerates for n → ∞.

Hence, our question: Could we find other transformations of the marginals that let us circumvent the case L → 0, and possibly find a new method or solution to handle this case? The relevance of this problem lies in the fact that relations between asymptotic dependence and the behavior of L are generally unknown: We know by formula ( 16) that the behaviour of L(n, n) and the values η takes are linked to the upper-tail dependence coefficient λ u , but the case L → 0 and η < 1 is not tackled. Moreover, we also need to consider L(n, ãn ) instead of L(n, n).

(b) [START_REF] Ledford | Concomitant tail behaviour for extremes[END_REF]: The focus here is only on the asymptotic distribution of the concomitant of the maximum X (n) , that is the asymptotic distribution of V 1 when k = 1 (in our notation), meaning to study F 1 . The authors explain that there can be cases where unit Fréchet marginals are not appropriate and suggest a standard Gumbel transformation for the variable Y, arriving to a new BSV function L 1 (in fact the first order approximation of L) that avoids the degeneracy of the limit of H 1 (with our notations). However, to compare this approach with the case where both the marginals are unit Fréchet, they need to assume y → ∞.

Let (X N , Y N ) be a bivariate normal random vector with standard normal marginals and with correlation 0 < ρ < 1. Theorem 3 in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF] states that, in the extreme case (k fixed, n → ∞), it holds that

P V 1 -ãn bn v 1 , V 2 -a n b n v 2 → n→∞ Φ k (v 1 )e -e -v 2 , ( 28
)
where

a n = (2 log n) -1/2 b n = 2 log n - 1 2 log (4π log n) 2 log n (29a) ãn = 1 -ρ 2 bn = ρb n . ( 29b
)
Formula ( 28) represents an asymptotic independence property of V 1 and V 2 , strongly connected to the fact that

F 1 ( ãn y + bn | a n x + b n ) = P Y N ãn y + bn | X N > a n x + b n → n→∞ Φ(y).
(30) Recall that our model in Theorem 1 is not applicable here. In fact, Condition (18) is not satisfied as L(n, ãn ) → 0 for any choice of the normalizing constants. Thus, further study would be needed in such a case. A way to look at the problem may be through [START_REF] Wadsworth | A new representation for multivariate tail probabilities[END_REF]. The authors study the effect of letting the components of a bivariate random vector grow at different rates, with the result of stating a new class of regular variation conditions and providing a new characterization of the link between the multivariate tail decay and the considered marginal growth rates. More precisely, the authors transform the marginals of the original vector into standard exponentials, called (X E , Y E ), and assume that, for all

(ζ 1 , ζ 2 ) ∈ R 2 + \ {0}, P(X E > ζ 1 log n, Y E > ζ 2 log n) ∼ n→∞ L(n; ζ 1 , ζ 2 )n -κ(ζ 1 ,ζ 2 ) ( 31
)
where L is a univariate SV function in n as n → ∞, and the function κ(ζ 1 , ζ 2 ) maps the different marginal growth rates to the joint tail decay rate with the role of providing information about the level of dependence between variables at subasymptotic levels. Assumption ( 31) is then used to describe models for conditional probabilities of extreme exceedances, i.e.

P(X

E > ζ 1 log n | Y E > log n) ∼ n→∞ L(n; ζ 1 , 1)e (1-κ(ζ 1 ,ζ 2 )) log n . ( 32
)
Studying ( 32) is proven by the authors to be equivalent to the search for functions a * , b * and h such that

lim n→∞ P X E -b * (log n) a * (log n) > x Y E > log n = e -h(x) . ( 33
)
The use of the exponential marginals is helpful to identify the normalizing constants for the bivariate Gaussian example, thus partially filling the gap left open in [START_REF] Ledford | Concomitant tail behaviour for extremes[END_REF], where, as previously explained, a way to tackle the problem is given in the case of Fréchet-Gumbel marginals, albeit that solution can be used only for some high values of Y. However, this discussion

does not yet answer our main question: Why, when using unit Fréchet marginals, aren't we able to find normalizing constants that lead to a non-degenerate limit for L(nx, ãn y + bn )? Our understanding of the problem is made more complete by the insightful paper by [START_REF] Heffernan | Limit laws for random vectors with an extreme component[END_REF]. There, the authors explore the implications of assuming the existence of a scaling function c * 1 (•), a centering function c * 2 (•), and a non-null Radon measure µ on Borel subsets of [-∞, ∞] × (-∞, ∞] such that, for each fixed y ∈ E ξ := {y ∈ R : 1 + ξy > 0} (where ξ is the shape parameter of the extreme value distribution relative to F Y ), it holds a) µ([-∞, x] × (y, +∞)) is a finite and non-degenerate distribution function in x b) the limit of the cdf of the normalized X given that Y is extreme exists and

P X -c * 2 • b ← (t) c * 1 • b ← (t) x Y > t → t→y + µ([-∞, x] × (0, +∞)), (34)
where y + is the right-end point of F Y . It is worth noticing the fact that we are free to change the marginal of Y without disturbing the convergence in ( 34), but this is not true for the X variable. In the paper it is proven that, in the case of the bivariate standard Gaussian, it is not possible to obtain a non-degenerate expression for µ in (34) if the marginals of X and Y are both transformed into the Pareto type, while an asymptotic form is still reachable if the marginals are transformed into exponentials. We refer to Heffernan and Resnick (2007, Sec. 7) for a thorough discussion of this aspect.

.

Revisiting the Gaussian case

Now that we have a better and more complete understanding on the issues underlying the treatment made in the literature of the Gaussian case, we provide a way to find suitable constants through the application of model ( 32) and formulas (29a) and (29b). We know that, if (X N , Y N ) is a bivariate normal random vector and correlation ρ ∈ (0, 1),

P Y N > ãn y + bn | X N > a n x + b n → n→∞ 1 -Φ(y), (35)
where the constants are given in (29a) and (29b). The following results hold:

Theorem 2. Let (X E , Y E ) be a bivariate random vector with unit exponential marginals and correlation ρ ∈ (0, 1), such that (32) holds, with ζ 1 := ζ 1 (n, y). Then, we have

P Y E > ãn,E y + bn,E | X E > log n → 1 y √ 2π e -(ρy) 2 /2 , as n → ∞, (36)
where

ãn,E = ρ 1 -ρ 2 2 log n - 1 2 log (4π log n) 2 log n (37a) and bn,E = ρ 2 log n - ρ 2 2 log (4π log n) + ρ 2 16 (log (4π log n)) 2 log n . ( 37b
)
Corollary 2. Note that, if we assume y large enough (as in [START_REF] Ledford | Concomitant tail behaviour for extremes[END_REF]), then we can use the Mill's ratio to obtain

P Y E > ãn,E y + bn,E | X E > log n ∼ n→∞ ρ (1 -Φ(ρy)) . ( 38
)
Refer to Appendix A for the proofs.

The limit distribution given in ( 38) is similar to those in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF] and [START_REF] Wadsworth | A new representation for multivariate tail probabilities[END_REF]. The differences among the final formulas are due to the different chosen orders of the Taylor approximations.

The work presented has provided a deep understanding of the impact of the extremal dependence of a bivariate random vector (X, Y) on the joint asymptotic behavior of suitably rescaled maxima over subsets of concomitants. Studying the framework of asymptotically dependent random variables is actually of great interest in practical applications and completes the literature.

Our main result addresses this problem, providing the expression of the asymptotic joint cdf of the two rescaled maxima of concomitants. In particular, we take advantage of the joint-tail model proposed in [START_REF] Ledford | Modelling dependence within joint tail regions[END_REF] to provide the asymptotic form of all the conditional distributions describing Y given the behaviour of X, which contribute to the asymptotic joint cdf of the two maxima of concomitants. We also explain in detail how to choose suitable normalizing constants needed to avoid the degeneracy of such distributions (thus preventing the joint cdf of the maxima of concomitants to be degenerate as well). Our theoretical results are illustrated through two examples of interest in the risk analysis literature: We compute the analytical expressions of the conditional distribution and perform a numerical integration to obtain the joint cdf. We show the goodness of the asymptotic approximation in all the examined cases through a numerical simulation. Furthermore, we reconsider Theorem 2 in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF] and show how some hypotheses are connected to the upper-tail dependence coefficient characterizing the copula of the parent vector (X, Y), thus providing an intuitive interpretation of the formulas therein and a new reformulation of the cited Theorem. More importantly, we revisit the example of the bivariate Gaussian distribution in our case, furnishing a comprehensive discussion on it. We deeply study the reasons why the joint-tail model in [START_REF] Ledford | Modelling dependence within joint tail regions[END_REF] fails in providing Zamanzade, E. and X. Wang (2018). Proportion estimation in ranked set sampling in the presence of tie information. Computational Statistics 33(3), 1349-1366.

-Proof of Lemma 1. Using the notations of Lemma 1, and indicating with D(G γ 1 ), D(G γ 2 ), D(G 0 ) the MDAs of the Fréchet, the Weibull and the Gumbel distributions, respectively, we can write, for x, y ∈ R and x + denoting the upper end-point of the distribution of X,

nP(X > a n x + b n , Y > c n y + d n ) = x + a n x+b n nP(Y > c n y + d n | X = t) f X (t) dt n x + a n x+b n f X (t) dt = n(1 -F X (a n x + b n )) < ∞, ∀n ∈ N.
As shown in the proof of Proposition 2.5 in [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF], it holds

na n f X (a n z + b n )1 x z < a -1 n (x + -b n ) → n→∞ 1 (z x) g(z)
with

g(z) :=      γ 1 z -(γ 1 +1) , z > 0, γ 1 > 0, if F X ∈ D(G γ 1 ), γ 2 (-z) γ 2 -1 , z < 0, γ 2 > 0, if F X ∈ D(G γ 2 ), e -z , z ∈ R, if F X ∈ D(G 0 ).
Note the role of the indicator 1 x z < a -1 n (x + -b n ) , needed to ensure that we are integrating over the support of X. In order to apply the Extended Dominated Convergence Theorem, we need to check that

R na n f X (a n z + x + )1 (z x) dz = n (1 -F X (a n x + b n )) → n→∞ R g(z)1 (z x) dz.
For the Gumbel case, we have

n (1 -F X (a n x + b n )) → e -x and R g(z)1 (z x) dz = e -x .
The Fréchet and Weibull cases follow the same reasoning. Since the two quantities are equal, we can exchange limit and integral, as follows:

lim n→∞ R nP(Y > c n y + d n | X = t) f X (t)1 (t a n x + b n ) dt = R lim n→∞ P(Y > c n y + d n | X = a n z + b n ) • lim n→∞ na n f X (a n z + b n ) • 1 z < a -1 n (x + -b n ) • 1 (z x) dz.
Since the integrand is positive, then the integral tends to 0 if and only if the latter tends to 0, that is, if and only if

P(Y > c n y + d n | X = a n z + b n ) = 1 -F 3 (c n y + d n | a n z + b n ) → n→∞ 0 for all z x (note that lim n→∞ na n f X (a n z + b n ) = 0). Hence the result.
Proof of Lemma 2. In order to make use of ( 15), we write all the conditionals in terms of F and recall the notation C ⊥ ⊥ for the case α + β = 1 and λ U > 0. For what concerns F 2 (y | x), we have

F 2 (y | x) = P(Y y | X < x) = 1 -P(Y > y | X < x) = 1 - FY (y) -F(x, y) F X (x) . Since FY (ny) -F(nx, ny) F X (nx) = e 1 nx 1 -e -1 ny - L(nx, ny) n α+β x α y β ∼ n→∞ 1 n 1 + 1 nx 1 y - L(nx, ny) n α+β-1 x α y β = 1 ny 1 + 1 nx 1 - L(nx, ny) n α+β-1 x α y β-1 ∼ n→∞ 1 y 1 + 1 nx 1 - L(n, n) n α+β-1 r(x, ny)x -α y 1-β ∼ n→∞ 1 ny 1 -λ U r(x, y) y x α 1 C ⊥ ⊥ ,
we obtain

F n 2 (ny | nx) = 1 - FY (ny) F X (nx) - F(nx, ny) F X (nx) n ∼ n→∞ 1 - 1 ny 1 -λ U r(x, y) y x α 1 C ⊥ ⊥ n . Thus F n 2 (ny | nx) → n→∞ e -1 y 1-λ U r(x,y)( y x ) α 1 C ⊥ ⊥ .
For F 3 (y | x), note that

F 3 (y | x) = 1 f X (x) ∂ ∂x F(x, y) = 1 f X (x) ∂ ∂x F X (x) -F(y) + F(x, y) = 1 + 1 f X (x) ∂ ∂x F(x, y) = 1 + x 2-α y -β e 1 x ∂ ∂x L(x, y) - α x L(x, y) .
It follows that

F 3 (ny | nx) ∼ n→∞ 1 + n 2-α-β e 1 nx x 2-α y -β ∂ ∂(nx) L(nx, ny) - α nx L(nx, ny) = 1 + e 1 nx x 2-α y -β n 2-α-β 1 n ∂ ∂x L(nx, ny) - α nx L(nx, ny) ∼ n→∞ 1 + e 1 nx n 1-α-β L(n, n)x 2-α y -β ∂ ∂x r(x, y) - α x r(x, y) → n→∞ 1 + λ U x 2-α y -β ∂ ∂x r(x, y) - α x r(x, y) 1 C ⊥ ⊥ ,
which concludes the proof.

Proof of Theorem 1. The proof is based on Lemma 2, and [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF]. To begin with, we report the finite-sample bivariate cumulative distribution of (V 1 , V 2 ) (Joshi and Nagaraja (1995, Equation (2.1))),

F (V 1 ,V 2 ) (v 1 , v 2 ) = E h v 1 , v 2 , X (n-k) (39)
where the function h is defined as

h(v 1 , v 2 , x) = F k 1 (v 1 | x)F n-k-1 2 (v 2 | x) F 3 (v 2 | x) .
The rational under Equation ( 39) lies in [START_REF] Kaufmann | On conditional distributions of nearest neighbors[END_REF]. Conditioned on the event {X (n-k) = x}, the maxima V 1 and V 2 have the following behaviour: V 1 can be considered as the sample maximum of a random sample of size k from F 1 (• | x), while V 2 behaves like the sample maximum of a sample of size n -k, composed by n -k -1 i.i.d. variables from F 2 (• | x) and a single random variable from F 3 (• | x).

Following [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF], from Equation (39), we can write

F (V 1 ,V 2 ) ( ãn v 1 + bn , c n v 2 + d n ) = F 1 ( ãn v 1 + bn | a n x + b n )F n-k-1 2 (c n v 2 + d n | a n x + b n ) F 3 (c n v 2 + d n | a n x + b n ) a n f X (a n x + b n ) dx.
Since our goal is to compute the joint asymptotic distribution of V 1 and V 2 as n → ∞, we refer to [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF] and [START_REF] Nagaraja | Distribution of the Maximum of Concomitants of Selected Order Statistics[END_REF] for what concerns the swap between the limit operation and the integral. Since X is unit Fréchet distributed, it is clear that the marginal F X is in the domain of attraction of G X (x) = e -1/x and that g X (x) = x -2 e -1/x . We refer to Lemma 1 and Result 1 in [START_REF] Nagaraja | Distribution of the Maximum of Concomitants of Selected Order Statistics[END_REF], for the convergence as n → ∞ of the term a n f X (a

n x + b n ) to g (k) (x) := (-log G X (x)) k k! g X (x)
. Theorem 1 then follows by Condition (18) and Lemma 2.

Proof of Theorem 2. Note that X N > a n x + b n is equivalent to X N > b n , as a n → 0 when n → ∞. Let us now report the model in [START_REF] Wadsworth | A new representation for multivariate tail probabilities[END_REF] in terms of standard normal marginals, namely

P(Y E > ζ 1 (n, y) log n | X E > log n) = P Y N > Φ -1 1 -n -ζ 1 (n,y) X N > Φ -1 1 -n -1 ∼ L 1 (n; ζ 1 (n, y), 1) n 1-κ(ζ 1 (n,y),1) (40)
where

κ(ζ 1 (n, y), 1) = ζ 1 (n, y) + 1 -2ρ ζ 1 (n, y) 1 -ρ 2 , and (41)
L 1 (n; ζ 1 (n, y), 1) = (4π log n) 2ρ 2 -ρ( √ ζ 1 (n,y)+1/ √ ζ 1 (n,y)) 2(1-ρ 2 ) ζ 1 (n, y) 1-ρ/ √ ζ 1 (n,y) 2(1-ρ 2 ) (1 -ρ 2 ) 3/2 ( ζ 1 (n, y) -ρ)(1 -ρ ζ 1 (n, y)) , (42)
whenever ρ 2 < min{ζ 1 (n, y), 1}. The authors explain that the SV function L, which we have in the case of the standard exponential marginals (from Equation (32)), is such that L/L 1 ∼ 1 as n → ∞.

To ease the notation, in the following, we use ζ 1 for ζ 1 (n, y). We want to find some constants ãn,E and bn,E such that {Y E > ζ 1 (n, y) log n} is equivalent to {Y E > ãn,E y + bn,E }, so we set

ζ 1 = ãn,E y + bn,E log n . ( 43
)
By using the first equality in ( 40) and knowing that if

Y N ∼ N(0, 1) then -log (1 -Φ(Y N )) ∼ Exp(1), it follows that ζ 1 = -log 1 -Φ( ãn y + bn log n . ( 44
)
This means that we can rewrite (35) by transforming the marginals into standard exponential and then apply the model (40). Combining ( 43) and ( 44 Note that, up to now, we have not used the constants in (29b), but found some general conditions on the their growth as n increases. By using (29b), we have

ãn,E = ρ 1 -ρ 2 2 log n - 1 2 log (4π log n) 2 log n bn,E = ρ 2 log n - ρ 2 2 log (4π log n) + ρ 2 16 (log (4π log n)) 2 log n .
where we neglected all the terms characterized by a faster decay with respect to (log n) -1 , we obtain

P(Y E > ζ 1 (n, y) log n | X E > log n) ∼ n→∞ 1 y √ 2π e -(ρy) 2 /2 .
Proof of Corollary 2. Notice that

P(Y E > ζ 1 (n, y) log n | X E > log n) ∼ n→∞ ρ ρy √ 2π e -(ρy) 2 /2 = ρ ϕ(ρy) ρy ,
where ϕ denotes the probability density function of the standard normal distribution. Under the assumption that y is large, it is straightforward to show, by using the Mill's ratio, that

P(Y E > ζ 1 (n, y) log n | X E > log n) ∼ n→∞ ρ (1 -Φ(ρy)) ,
which completes the proof.

k

The present section contains some graphical results of the examples considered in Section 3. We recall that each example runs a simulation of 5000 replicates of 100000 bivariate sample (X, Y). As a result, we obtain a bivariate sample for (V 1 , V 2 ) of dimension 5000 (each replicate of the (X, Y) sample gives one couple of (V 1 , V 2 )), for each k. The values of k indicating the dimension of the extreme set are fixed to 1, 50, 100, 250, 500, 1000, 1500, 2000. However, the figures do not show the results for each of the considered values: This choice is made to foster the readability, as a lot of almost similar plots on a single page may be redundant. From a quick look at Figures 4 to 7 for both examples, we can draw the same important conclusions:

(a) The range of the entries relative to V 1 is much wider than the range of the entries corresponding to V 2 ;

(b) The cloud of black points in the surfaces in Figure 4 and Figure 6 becomes closer and closer to the v 1 axis as k increases.

This is a clear effect of the positive dependence between (X, Y). To see this, think about the definition of (V 1 , V 2 ): V 1 is the maximum of the concomitants of the largest k order statistics and V 2 is the maximum of the remaining set. The clear difference between the ranges of the two samples indicates that the highest concomitants contribute to the values of V 1 and the lowest concomitants are responsible of the relatively low values of V 2 . Moreover, as k increases, almost all the values that V 2 takes are lower than 1, while V 1 can assume values which are greater than 100, for example. Explained differently, high values of concomitants generally correspond to high order statistics of X. It has to be said that the plotting range of the surfaces is chosen in order to enhance the visualization purposes, as extremely wide intervals would have been not suitable to observe the results. However, there can be points which fall outside the plotting area: Figure 4e contains more points in the middleright part, while Figure 4f contains more points in the upper-right part. This indicates that when k is 1000, some concomitants high in value do not belong to the extreme set, but when k increases to 2000, the partition to which they belong changes. In other words, V 1 "steals" high values to V 2 when k increases. Plots relative to the marginal distributions support the presented conclusions: There is no noticeable difference among the log cdf shown in Figure 5a and Figure 7a, which means that the sample in V 1 does not change so much, while the upper end point of the distribution of V 2 appears to shrink towards 0 (Figure 5b and Figure 7b). This is again an indication of the "stealing" argument: Since high values of concomitants contributing to V 2 pass into the extreme set as k increases, only low values of concomitants are left to compute V 2 , which becomes smaller and smaller. As previously said, this is the effect of the positive dependence between (X, Y). Finally, the extremely low errors shown in the boxplots in Figure 4g and Figure 6g ensure that the asymptotic approximation is almost perfect. The same conclusion can be made from the plots sowing the marginal errors in Figure 5 and Figure 7. Despite some variability for low values of v 1 and v 2 (which is small anyway), the errors rapidly approach 0 as v 1 and v 2 increase, which is important for practical applications as we are usually interested in studying what happens over high thresholds.

-This section reports the Von Mises conditions, as written in [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF], borrowing his notation.

Proposition 1.15 [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]): Suppose F is absolutely continuous with positive density f in some neighborhood of ∞. (c) Equation (47) holds iff for some z 0 and all x > z 0 , we have

1 -F(x) = c exp - x z 0 t -1 α(t) dt
where lim t→∞ α(t) = α and c is a positive constant.

Proposition 1.16 [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]): Suppose F has finite right end-point x 0 and is absolutely continuous in a left neighborhood of x 0 with positive density f. where lim t↑x 0 δ(t) = α, lim t↑x 0 c(t) = c 0 > 0.

Proposition 1.17 [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF]): Let F be absolutely continuous in a left neighborhood of x 0 with density f. (g(t)/f(t)) dt , z 0 < x < x 0 , (50)

where lim t↑ 0 g(x) = 1 and f is absolutely continuous with density f(x) → 0 as x ↑ x 0 .

(d) Equation (49) or (50) are equivalent to tf ((1/(1 -F)) ← (t)) ∈ RV 0 [START_REF] Sweeting | On Domains of Uniform Local Attraction in Extreme Value Theory[END_REF]).

  y) x 1-α < y β , c(x, y) ∼ y→0 y -β , c(x, y) = o y -β for y → ∞, with ãn = O n 1-α β and bn = O ( ãn ) or o ( ãn ) .

Figure 1 :

 1 Figure 1: Joint tail F(x, y) (gray) with unit Fréchet marginals, as expressed in (24), and its approximation (black and gold) given by Equation (15). The parameters α and β are both set to 0.5 in the two examples.

  Marginal distributions diagnostics.

Figure 2 :

 2 Figure 2: Survival Clayton(θ = 2), with k = 10. Panel (a) shows the empirical univariate log cdf of V 1 (in black), together with the asymptotic theoretical log values (in gold) obtained by numerically integrating the asymptotic joint cdf (19) to obtain the marginal. Panel (b) shows the absolute values of the errors between empirical and theoretical cdf of V 2 . The upper plot of panel (c) shows the empirical univariate cdf of V 1 (in black), together with the asymptotic theoretical values (in gold), in log scale, obtained by numerically integrating the asymptotic joint cdf (19). The lower plot of panel (c) represents the L 1 error between the marginal empirical and theoretical cdf.

  Marginal distributions diagnostics.

Figure 3 :

 3 Figure 3: Logistic(γ = 1/2), with k = 10. Panel (a) shows the empirical univariate log cdf of V 1 (in black), together with the asymptotic theoretical log values (in gold) obtained by numerically integrating the asymptotic joint cdf (19) to obtain the marginal one. Panel (b) shows the absolute values of the errors between the empirical and the theoretical cdf of V 2 . The upper plot of panel (c) shows the empirical univariate cdf of V 1 (in black), together with the asymptotic theoretical values (in gold), in log scale, obtained by numerically integrating the asymptotic joint cdf (19). The lower plot of panel (c) represents the L 1 error between the marginal empirical and theoretical cdf.

  Boxplots of the errors.

Figure 4 :

 4 Figure 4: Survival Clayton(θ = 2). Bivariate cdf computed by numerically integrating equation (19), for various values of k. The black bullets correspond to the bivariate empirical cdf of the simulated samples. Panel g) shows the boxplots of the errors in absolute value between the bivariate asymptotic cdf computed in (V 1,i , V 2,i ) for i = 1, . . . , 5000 and the corresponding empirical bivariate cdf. The gold dots represent the error means.

  V 1 : log marginal cdf comparison and marginal error behaviour with respect to v 1 . Error (absolute value) (b) V 2 : log marginal cdf comparison and marginal error behaviour with respect to v 2 .

Figure 5 :

 5 Figure 5: Survival Clayton(θ = 2). Empirical univariate cdf of V 1 and V 2 , together with the asymptotic theoretical values obtained by numerically integrating the asymptotic joint cdf (19).

  Boxplots of the errors.

Figure 6 :

 6 Figure 6: Logistic(γ = 1/2). Bivariate cdf computed by numerically integrating equation (19), for various values of k. The black bullets correspond to the bivariate empirical cdf of the simulated samples. Panel g) shows the boxplots of the errors in absolute value between the bivariate asymptotic cdf computed in (V 1,i , V 2,i ) for i = 1, . . . , 5000 and the corresponding empirical bivariate cdf. The gold dots represent the error means.

  V 1 : log marginal cdf comparison and marginal error behaviour with respect to v 1 . Error (absolute value) (b) V 2 : log marginal cdf comparison and marginal error behaviour with respect to v 2 .

Figure 7 :

 7 Figure 7: Logistic(γ = 1/2) example. Empirical univariate cdf of V 1 and V 2 , together with the asymptotic theoretical values obtained by numerically integrating the asymptotic joint cdf (19).

  then F ∈ D (Φ α ). We may choose a n to satisfy a n f (a n ) ∼ α/n. (b) If f is nonincreasing and F ∈ D (Φ α ) then (47) holds.

  -x) f(x)/(1 -F(x)) = α (48) then F ∈ D (Ψ α ). (b) If f is nonincreasing and F ∈ D (Ψ α ) then (48) holds.(c) Equation (48) holds iff c(x) can be taken to be constant in some left neighborhood of x 0 in the representation1 -F(x) = c(x) exp -x x 0 -1δ(t)/ (x 0 -t) dt and for x < x 0

  F(s)) ds/(1 -F(x)) b n = (1/(1 -F)) ← (n), a n = f (b n ) . (b) If f is nonincreasing and F ∈ D(Λ) then holds. (c) Equation (49) holds iff 1 -F(x) = c exp -x z 0

1

  Introduction 2 Joint asymptotics for maxima over subsets of concomitants 2.1 The LT joint-tail model . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Discussion on the choice of ãn and bn . . . . . . . . . . . . . . . . 3 Illustration of the results 3.1 Pareto -Lomax marginals with Survival Clayton copula . . . . . 3.2 Unit Fréchet marginals with Logistic Extremal copula . . . . . . . 4 Discussion on the Gaussian case 4.1 Existing literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Revisiting the Gaussian case . . . . . . . . . . . . . . . . . . . . . .

	5 Discussion
	a Appendix -Proofs
	b Appendix -Discussion on k
	c Appendix -The Von Mises conditions

a non-degenerate asymptotic distribution and fill the gap in the extension of the joint-tail model proposed in [START_REF] Wadsworth | A new representation for multivariate tail probabilities[END_REF], by showing how to compute suitable normalization constants in this case. Our results are coherent with those obtained by direct computation in [START_REF] Joshi | Joint distribution of maxima of concomitants of subsets of order statistics[END_REF].

The authors are grateful for funding to Labex MME DII (ANR-11-LABX-0023-01) and ESOMAS department (University of Turin). This paper was mainly written during the first author's visiting period at the ESSEC Business School, Paris, France and during the second author's visiting period at the ESOMAS Department, University of Turin, Turin, Italy. Both authors are grateful to the hosting institutions.

Replacing the constants in (43) it yields

Now, notice that with ζ 1 defined in ( 46), the condition ρ 2 < min{ζ 1 (n, y), 1} is always satisfied for large n and for ρ ∈ (0, 1), as it yields ρ 2 < ρ 2 + O (log n) -1/2 . Hence, since

By using these expressions into Equation ( 40), we have