G Kariniotakis
email: georges.kariniotakis@ensmp.fr

P Pinson

Evaluation of the MORE-CARE Wind Power Prediction Platform. Performance of the Fuzzy Logic Based Models

Keywords: Wind power, time-series forecasting, numerical weather predictions, on-line software, adaptive fuzzy-neural networks

The paper presents an advanced wind forecasting system that uses on-line SCADA measurements, as well as numerical weather predictions as input to predict the power production of wind parks 48 hours ahead. The prediction tool integrates models based on adaptive fuzzy-neural networks configured either for short-term or long-term forecasting. In each case, the model architecture is selected through non-linear optimization techniques. The forecasting system is integrated within the MORE-CARE EMS software developed in the frame of a European research project. Within this on-line platform, the forecasting module provides forecasts and confidence intervals for the wind farms in a power system, which can be directly used by economic dispatch and unit commitment functions. The platform can run also as a stand-alone application destined only for wind forecasting. Detailed results are presented on the performance of the developed models over a one-year evaluation period on five real wind farms in Ireland, using HIRLAM numerical weather prediction and SCADA data as input.

I. INTRODUCTION

he large-scale integration of wind power in interconnected or isolated power systems emerges the development of appropriate tools to assist the wind farm operators on their daily management task. Short-term forecasts of the wind farms production, up to 48 hours ahead, are necessary for a secure and economic large-scale wind power integration. Wind power prediction tools are useful for end-users such as Independent Power Producers, Transmission and Distribution System Operators, Energy Service Providers, Traders a.o. In a liberalised electricity market environment, such tools enhance the competitiveness of wind power, since they reduce the penalties resulting from the wind resource intermittence. Reduced operational and financial risk for the wind farm developers is a motivating factor for undertaking investments on wind farms.

Wind power forecasting is a far from trivial problem. Wind speed is a non-stationary process both in the mean and variance. Wind power is nonlinear w.r.t. speed with a major difficulty in the area of cut-off speed, where prediction intervals can extend from maximum to zero wind power.

Among the difficulties, one should add the error of numerical weather forecasts, which are often used as input to the models. Often, no adequate information is available online by a data acquisition system (SCADA) to assess the actual operational status of the wind farm (i.e. how many turbines are in operation). The available on-line data can be detailed (i.e. power, speed of each wind turbine) or not (i.e. only total power available). In some situations there is complete lack of data and information from neighbor wind farms has to be assessed.

Research on wind power forecasting is actively pursued by several research centres in Europe. Actually there are two main stateof-the-art approaches; one based on physical or deterministic modelling and a second one based on statistical or timeseries modelling.

The "physical" approach for wind power forecasting is based on a detailed description of the wind park site (orography, roughness, obstacles), the wind turbines (hub height, power curve, thrust curve) and the layout of the wind plant. In [START_REF] Giebel | Stateof-the-Art on Methods and Software Tools for Short-Term Prediction of Wind Energy Production[END_REF][START_REF] Focken | Reduction of wind power prediction error by spatial smoothing effects[END_REF], wind power forecasting platforms based on physical methods are described. The main input is numerical weather predictions (NWP). Model output statistics are developed to account for systematic errors. Weather predictions are however updated only a limited number of times per day by meteorological services. For this reason, the performance of these models is often satisfactory for rather longer (>6 hours ahead) than short-term horizons.

The alternative "time series", or statistical, approach includes typical linear models (ARMA, ARX etc) and non-linear ones (i.e. neural networks, conditional parametric models, etc). These models aim to predict the future by "capturing" temporal and spatial dependencies in the data [START_REF] Giebel | Stateof-the-Art on Methods and Software Tools for Short-Term Prediction of Wind Energy Production[END_REF], [START_REF] Fukuda | The development of a wind velocity prediction method based on a Data-Mining type Auto-Regressive model[END_REF][START_REF] Nielsen | Prediction of regional wind power[END_REF][START_REF] Kariniotakis | Wind power forecasting using advanced neural network models[END_REF][START_REF] Kariniotakis | Advanced Short-Term Forecasting of Wind Power Production[END_REF]. The input to these models can be on-line SCADA data and numerical weather predictions (NWP). For look-ahead times more than ~10 hours (mentioned hereafter as "long-term"), NWPs are indispensable for an acceptable performance, since they represent weather dynamics that cannot be modelled using only recent on-line data. For shorter horizons, up to ~10 hours ahead (mentioned hereafter as "short-term"), time series models can be based exclusively on recent measurements; however even in this case, NWPs as explanatory input improves results. It is noted that the threshold of 10 hours is mentioned as an example rather than a rule, since it depends on the characteristics of a specific wind profile.

The models presented in this paper belong to the time series approach. In previous work of the authors, linear autoregressive models, radial basis functions, wavelet networks, feed forward and recurrent neural networks and finally adaptive fuzzy-neural network models were compared for the task of short-term prediction [START_REF] Kariniotakis | Advanced Short-Term Forecasting of Wind Power Production[END_REF][START_REF] Dutton | Load and Wind Power Forecasting Methods for the Optimal Management of Isolated Power Systems with High Wind Penetration[END_REF][START_REF] Kariniotakis | Evaluation of advanced wind power and load forecasting methods for the optimal management of isolated power systems[END_REF]. Fuzzy neural networks, originally used here for wind forecasting, were found to outperform the other approaches in both tasks of shortterm (0-10 hours) and longer-term (0-48 hours) wind prediction.

The developed models have been implemented in an online prediction module, which is the core of the Armines Wind Power Prediction System (AWPPS) developed by Ecole des Mines de Paris. The AWPPS platform integrates a relational database for data management and a Java Man-Machine Interface (MMI) for higher portability. Communication between the modules is performed using ODBC/JDBC.

T

The core prediction module of AWPPS has been integrated: (i) in the More-Care Energy Management System (EMS) developed in the frame of a European project and installed for on-line operation in the islands of Crete, Madeira and in Ireland (see Figure 8);

(ii) in an industrial SCADA software (see Figure 9). This paper presents results from the implementation of More-Care in Ireland where the aim is the prediction of the output of 11 wind farms of a total capacity of several tens of MWs.

II. THE AWPPS WIND PREDICTION SYSTEM.

The core prediction module of AWPPS provides forecasts for multiple wind farms. It integrates several functionalities such as:

1. short-term models considering only SCADA data as input based on adaptive fuzzy-neural networks (F-NN). Such models provide predictions for 1-10 hours ahead and are destined for small power systems. 2. Longer-term models based on adaptive F-NNs for 0-48 hours ahead using NWPs as input and on-line data if available. 3. Simple alternative models (persistence, moving average techniques, power curve models based on direct conversion of NWPs to power), which can be optionally activated as backup models. 4. Methods to estimate the uncertainty of the predictions. 5. A method to combine forecasts by different models. 6. Scheduled maintenance of the wind turbines, etc.

Adaptive fuzzy-neural networks have been used for both shortterm and long-term wind power prediction. The adaptivity property stands for the capacity of the model to fine-tune its parameters during on-line operation. This is an important requirement for a nonstationary process like wind speed or power. Adaptivity of the model compensates changes in the environment of the application that may happen during the lifetime of a wind farm. Such changes can be changes in the number of wind turbines (extension of the wind farm, maintenance or availability of the machines that is usually not available through SCADA), in the performance of the wind turbines due to aging, changes in the surrounding of the wind park (i.e. vegetation), or changes in the configuration of the model used to produce the NWPs.

The generic F-NN model, described in [START_REF] Kariniotakis | An advanced on-line wind resource prediction system for the optimal management of wind parks[END_REF], can be trained on appropriate input depending on the final use, which can be either short-term or long-term prediction.

A. Short-term models based on on-line data.

Short-term models receive historic values of wind power as input, as well as explanatory data, such as wind speed and direction, to predict wind power. The general form of a simple model with input only past values of power is:

( ) ( ) ( ) ( ) ( ) m t P t P t P f t P - - = + , , 1 , 1 ˆK
The function f(.) stands for the generic fuzzy-neural function. Multi-step ahead forecasts are generated using the model in an iterative way. I.e., in order to produce a forecast for t+2, the forecast for t+1 is fed back as input to the model. This approach presents the drawback that does not permit to iterate explanatory input, since no forecasts can be available for such quantities. To handle this problem, models using the look-ahead time k as an input variable can be considered. An alternative approach is to develop multi-output models, or to tune a different model for each time-step. The implementation of this approach is complex and requires high development effort, which can be prohibitive in case of a large number of wind farms.

The short-term models based on fuzzy-neural networks are found to outperform Persistence up to 20% according to the time-step [START_REF] Dutton | Load and Wind Power Forecasting Methods for the Optimal Management of Isolated Power Systems with High Wind Penetration[END_REF][START_REF] Kariniotakis | Evaluation of advanced wind power and load forecasting methods for the optimal management of isolated power systems[END_REF]. Persistence is a simple approach used as reference to evaluate the performance of advanced models. It assumes that the "wind in the future will be the same as the wind now".

Short-term predictions are adequate for small applications, e.g. islands, for which NWPs are not available. However, in larger systems, where predictions up to 48 hours are required, timeseries models based on meteorological information, as the one presented below, are necessary. Such models, due to the fact that consider NWPs as input, manage to double the improvement w.r.t. Persistence for the first 10 hours (up to 40% versus 20% for short-term models).

B. Longer-term models based on meteorological information.

For "long-term" horizons up to 24-48 hours ahead, it is necessary to include numerical weather predictions (NWP) as explanatory input to the model in order to have an acceptable performance. NWPs include usually wind speed, direction and temperature at 10 m and at several levels defined by atmospheric pressure. NWPs can be provided for the geographical coordinates of the wind farm or for a grid of four points surrounding the farm. In the second case, the spatial resolution of the NWP model is of primary importance. Meteorological models with high resolution are often more accurate but require high computation time to produce forecasts, and as a consequence, they do not update frequently their output (i.e. 1-4 times per day). In contrast, forecasts from low-resolution NWP models are more frequently available.

The developed forecasting tool is able to operate with input from different NWP systems. In the frame of the More-Care project it was tested and gave satisfactory results with input from the SKIRON system for the case of Crete, and also from HIRLAM for the case of Ireland. SKIRON forecasts were provided for a grid of 15x15 km (System B in Fig. 1), while HIRLAM predictions were provided at the level of the wind farm as interpolated values but from a model resolution of 33 km (System A in Fig. 1).

The developed model receives on-line data as well as NWPs as input to predict the wind farms production for the next 48 hours. These forecasts are updated every hour based on the most recent wind power measurements. Wind power data are necessary for the on-line updating procedure, independently if they are used or not as input variables to the model. The updating procedure permits mainly a good performance of the model for the first hours (i.e. 1-6 hours) of the considered horizon. Model configurations that do not update their forecasts based on recent wind power data were found to perform worse than Persistence in look-ahead times up to 6 hours ahead. Finally, the consideration of on-line information, other than wind power (i.e. wind speed or direction), was not found to contribute in the accuracy of the results. The general scheme of the model is shown in Fig. 1.

The aim of the prediction model is to capture the relations between input (meteorological information, on-line data) and output (future total wind park power). Such mapping includes the following implicit relations:

• Temporal correlations between past and future data of the process (autoregressive aspect of the model).

• Conversion of wind speed (meteorological predictions) from the height or the atmospheric level they are given to the hub height of the wind turbines.

• Spatial projection of the meteorological wind speed forecasts from the NWP grid points (e.g. 15x15 km) to the level of the wind farm ("downscaling").

• Correction of the wind park output for factors affecting the total production (i.e. array effects, effect of wind direction etc). The advantage of a model such as the fuzzy neural network model, compared to models based on the "physical" approach, is that it permits to avoid all the above intermediate modeling steps. Moreover, its adaptive mode can compensate situations like the ones explained in the previous Section.

C. Integration in the More-Care Energy Management System.

The core prediction Module, with the above functionalities ( 1)-( 6) has been integrated in the More-Care EMS. More-Care consists of a number of Modules for functions such as load, wind or hydro forecasting, economic dispatch, unit commitment, fast security assessment etc. The general architecture is shown in Figure 2. Numerical Weather Predictions are entered in the system via FTP connection, email etc. The short-term models, when activated, produce forecasts every 10-20 minutes for the next hours (sliding window scheme). The long-term models produce forecasts also with a sliding window scheme of one hour for the next 48 hours. The sliding window operation is a major difference from existing systems that produce forecasts only when new NWPs arrive (i.e. once or four times per day).

D. Off-line tuning of the prediction models.

Any wind power prediction software is not "plug-and-play" since it is always site-dependent. In order to run with acceptable accuracy when installed to a new site, it is always necessary to devote considerable effort for tuning off-line the prediction model on the characteristics of the local wind profile or for describing the environment of the wind farms. This task requires considerable expertise.

The outstanding problem in time series forecasting is to define the structure of the forecasting model for which optimal accuracy will be obtained. Optimal accuracy is required when predicting new data that is, data that have not been used during the model development. The capacity of the model to predict such "out of sample" data is called "generalization" and is a primary requirement during on-line operation.

Nowadays, models based on artificial intelligence are adopted in several prediction applications. This is because they permit to consider easily available explanatory input when this is available. In the case of wind power prediction, input can be past measurements of wind power, speed or direction, measurements from neighbor sites, numerical weather predictions of wind speed, direction, temperature etc for various levels of the NWP system and for grid points around the wind farm. A model however that uses excessive input will have a very high number of parameters and as a consequence will "learn" not only the useful information included in the data but also their noise ("overfitting"). Then, the capacity of that model to generalize will be low. As a conclusion: one should be cautious of models that use huge quantities of data as input. Models that use selectively and in an intelligent way the available input are expected to be more accurate. This is the well-known principle of "parsimony" in timeseries forecasting.

The off-line tuning of the prediction models considered here is based on an advanced methodology, which permits to select the optimal structure of a model and the most relevant input based on non-linear constrained optimization techniques. The estimation of the parameters is performed using learning algorithms that optimize simultaneously both the error response and the Information content of the model. The off-line methodology is presented in detail in [START_REF] Kariniotakis | Contribution to the development of an advanced control system for the optimal management of autonomous winddiesel systems[END_REF][START_REF] Kariniotakis | An advanced on-line wind resource prediction system for the optimal management of wind parks[END_REF].

III. EVALUATION RESULTS

Evaluation results are presented here for five wind farms in Ireland (WF-A to WF-E) with a total installed power of a few tens of MW. The available time series cover a period of almost two years from which 6600 hours were used for training (learning set), 1000 hours for cross-validation and one year for testing the performance of the models. Hirlam Numerical Weather Predictions (speed, direction and temperature at 10m or 2m and at model Levels 22, 23 and 24) are considered. Their resolution is 33 km. Only the total production of the wind farms is measured and used as input to the models.

The distribution of prediction errors varies as a function of the look-ahead time. As the prediction horizon increases, the standard deviation of the error distribution augments and the percentage of occurrence of low prediction errors diminishes significantly.

Figure 3 ahead). The bar width represents 5% of the wind farm nominal power. These are representative error distributions for this advanced model. As one can notice, for the first look-ahead time, almost 70% of the prediction errors are lower than 7.5% of the wind farm nominal power, though for 24-hour ahead prediction, there are 38% of the prediction errors that are inferior to that threshold. Figure 4 illustrates the performance of the advanced models versus that of Persistence for two wind farms (WF-A and WF-E).

The criteria used are the NMAE (Normalized Mean Absolute Error) and NRMSE (Normalized Root Mean Square Error). Normalization is based on the wind farm nominal capacity. The NRMSE criterion has higher values than the NMAE since larger errors weight more than small errors in contrast to the NMAE where all errors weight equally. In some situations, although NRMSE can provide some improvement w.r.t., Persistence this may not be the case for the NMAE criterion. Thus both criteria need to be used for a thorough evaluation.

The F-NN model always outperforms Persistence whatever the prediction horizon. This is also true for the first 3-4 look-ahead hours thanks to the use of SCADA data as input to the model. The NRMSE error of Persistence ranges between 9-39% of the nominal power, while that of the advanced model varies between 9-22%. The NMAE takes values between 5-16% according to the time step. A detailed analysis of the errors led to the conclusion that a major part of the uncertainty comes indeed from the Hirlam NWPs. Also given the fact that prediction of individual wind farms is performed, no spatial smoothing is present as would be the case in prediction of regional or national power. The level of accuracy increases when one considers the total number of wind farms. Then, average prediction error for the first 24 hours is inferior to 10%. The improvement obtained by the Fuzzy-NN model w.r.t. to Persistence is depicted in Fig. 5 where the improvement is plotted for the five wind farms. In all cases the advanced model: (i) is better than persistence, even for the first hours, (ii) outperforms Persistence up to 50% depending on the time step.

Since the evaluation of the prediction model is done over a oneyear period, it is of particular interest to visualize the monthly performance (Fig. 6). This performance is slightly variable from a month to another; for instance the prediction NMAE error ranges from 9.7% in July, up to 14.5% in May, for the 6-hour ahead horizon. In the results for each month of the year for all time steps the improvement obtained over persistence rises up to 60%. However, for several months, such as in August, this improvement is very low and this penalizes the overall result. This is because during this month low and relatively constant wind regimes are present. As a conclusion to this, further analysis could be based on a normalization of the criteria using average wind production as an alternative to nominal capacity. Finally, Fig. 7 depicts an episode with the wind power predictions for the next 43 hours compared to the real values for WF-E. The 85% confidence intervals are built with the method described in [START_REF] Pinson | On-line assessment of prediction risk for wind power production forecasts[END_REF].

IV. CONCLUSIONS

The paper presents results from the application of the fuzzyneural network prediction model for the case study of Ireland. By using both SCADA data and Hirlam Numerical Weather Predictions, the F-NN model can provide high quality forecasts. The model presents a clear advantage for the first hours [START_REF] Giebel | Stateof-the-Art on Methods and Software Tools for Short-Term Prediction of Wind Energy Production[END_REF][START_REF] Focken | Reduction of wind power prediction error by spatial smoothing effects[END_REF][START_REF] Fukuda | The development of a wind velocity prediction method based on a Data-Mining type Auto-Regressive model[END_REF][START_REF] Nielsen | Prediction of regional wind power[END_REF][START_REF] Kariniotakis | Wind power forecasting using advanced neural network models[END_REF][START_REF] Kariniotakis | Advanced Short-Term Forecasting of Wind Power Production[END_REF] where it outperforms Persistence. Emphasis is given in optimizing the model architecture using nonlinear constrained optimization techniques. By this way the most relevant input is automatically selected and the models obtained are parsimonious with optimal generalization capability.

The level of accuracy is satisfactory for single wind farm prediction. Part of the inaccuracy comes however from the NWP system and is mainly due to phase errors of the Hirlam predictions. On the other hand, the resolution of the available NWPs (33 km) was quite low. In the last months the resolution changed to 14 km and this is expected to improve accuracy in the future. When summing predictions for all wind farms, due to spatial smoothing effect of errors, the accuracy increases; the average error when predicting the total power reduces to less than 10% for the first 24 hours.

The developed methods have been implemented in operational software and installed for on line operation in Ireland as well as in Crete and Madeira in the frame of the European Project More-Care. Figure 8 shows a view of the Man-Machine interface of the More-Care implementation. Alternatively, the core prediction module has been also integrated in an industrial SCADA software (Figure 9) for predicting the output of individual wind farms connected in a power system.

Fig. 1 :

 1 Fig. 1: General scheme of the "long-term" prediction model with examples of two configurations of NWP systems used as input (SKIRON, HIRLAM).

Fig. 2 :

 2 Fig. 2: General architecture of the More-Care EMS system.

Fig. 3 :

 3 Fig. 3: Distributions of prediction errors for two prediction horizons (left: 1-hour ahead prediction error distribution, right: 24hour ahead prediction error distribution). Results are for WF-A.

Fig. 4 :

 4 Fig. 4: Error (NMAE, NRMSE) as a function of the look-ahead time for WF-A and WF-E as obtained by the F-NN model and persistence.

Fig. 5 :

 5 Fig. 5: Performance of the Fuzzy-NN model: improvement w.r.t. Persistence for the 5 wind farms over a one-year evaluation period.

Fig. 6 :

 6 Fig. 6: NMAE performance of the F-NN model and Persistence for WF-E for 6 hours ahead forecasts. Performance is given as a function of the month of the year (Oct. 2001-Aug. 2002).

Fig. 7 :

 7 Fig. 7: Example of wind power prediction by the Fuzzy-NN model with 85% confidence intervals (WF-E).

Fig. 8 :

 8 Fig. 8: User interface of the wind power prediction module of More-Care.

Fig. 9 :

 9 Fig. 9: User interface of an industrial SCADA software integrating the wind power prediction module.

V. ACKNOWLEDGMENTS

The authors gratefully acknowledge ESB National Grid for providing data for the realization of this work. This work was performed in the frame of the Project MORE-CARE (ERK-CT1999-00119) funded in part by the European Commission.


Armelle Guillou
email: armelle.guillou@math.unistra.fr

Simone A Padoan
email: simone.padoan@unibocconi.itstefano.rizzelli@phd.unibocconi.it

Stefano Rizzelli

Inference for asymptotically independent samples of extremes

Keywords: MSC 2010 subject classifications: Primary 62G32, 62G05, 62G20; secondary 60F05, 60G70 Extremal dependence, extreme-value copula, nonparametric estimation, Pickands dependence function

An important topic of the multivariate extreme-value theory is to develop probabilistic models and statistical methods to describe and measure the strength of dependence among extreme observations. The theory is well established for data whose dependence structure is compatible with that of asymptotically dependent models. On the contrary, in many applications data do not comply with asymptotically dependent models and thus new tools are required. This article contributes to the methodological development of such a context, by considering a componentwise maxima approach. First we propose a statistical test based on the classical Pickands dependence function to verify whether asymptotic dependence or independence holds. Then, we present a new Pickands dependence function to describe the extremal dependence under asymptotic independence. Finally, we propose an estimator of the latter, we establish its main asymptotic properties and we illustrate its performance by a simulation study.

Introduction

Multivariate extreme-value theory provides the mathematical foundation for performing real data analysis of rare events. To characterize the joint upper tail of a multivariate distribution, two different approaches can be used: either by considering the componentwise maxima, or all the observations above a high threshold. A description of these methodologies can be found for instance in Coles (2001, Ch. 8), Beirlant et al. (2004, Ch. 8-9), de Haan and Ferreira (2006, Ch. 6) and Resnick (2007, Ch. 6), among others. Unfortunately the flexibility of the dependence structures provided by the classical theory of multivariate extreme-values may not be sufficient for statistical modelling (see e.g. Ledford andTawn, 1996, 1997). To solve this issue, different coefficients of tail dependence or probabilistic models have been introduced. They allow to govern/describe the strength of the extremal dependence. In this paper, we are particularly interested in the notion of asymptotic independence which is common in real data analysis. This concept can be defined as follows.

Let Y be a multivariate random vector of dimension d, with distribution function F and marginals F j , 1 ≤ j ≤ d. We say that F is in the max-domain of attraction of a multivariate extreme-value distribution G, if there exist sequences of constants a n > 0 and b n ∈ R d such that lim n→∞ F n (a n y + b n ) = G(y), for all y ∈ R d . Under this condition, a particular case arises when G is equal to the product of its marginal distributions. In this setting, we say that Y satisfies the property of asymptotic independence (or tail independence) which is equivalent to say that the elements of Y are asymptotically independent in the upper tail, i.e. lim u→1 Pr(F j (Y j ) > u|F i (Y i ) > u) = 0 for all 1 ≤ i j ≤ d. On the contrary, if the above limit is positive, then the elements of Y are called asymptotically dependent. The classical theory expects asymptotic dependence and independence as the only two possible scenarios, conceiving extremes as independent in the second case. Many efforts have been made to characterize a residual tail dependence in the data (if there is any) by offering new dependence coefficients or probabilistic and statistical models under asymptotic independence, see [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF], Coles (2001, Ch 8.4), [START_REF] Resnick | Hidden regular variation, second order regular variation and asymptotic independence[END_REF], [START_REF] Maulik | Characterizations and examples of hidden regular variation[END_REF], Ramos andLedford (2009, 2011), [START_REF] Wadsworth | A new representation for multivariate tail probabilities[END_REF] and [START_REF] Wadsworth | Modelling across extremal dependence classes[END_REF], to name a few. If we restrict our framework to the dimension d = 2, several statistical tests for checking asymptotic independence or tail independence have been proposed, among them, [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF], [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], [START_REF] Hüsler | Testing asymptotic independence in bivariate extremes[END_REF] and Falk et al. (2010, Ch. 6.5) and the references therein. However, the extension to dimensions higher than 2 are still in its infancy. Recent proposals are based on the kth largest order statistics of the sample. Although these approaches are simple to implement, the performance of the resulting tests depends strongly on the choice of k, see e.g. [START_REF] Kiriliouk | Nonparametric estimation of extremal dependence In Extreme Value Modeling and Risk Analysis[END_REF].

In this paper, we propose in Section 2 an alternative approach to test asymptotic independence for an arbitrary dimension d ≥ 2, based on the componentwise maxima. We illustrate the performance of our proposal up to dimension d = 4. Then, using again the componentwise maxima approach and in particular the framework proposed by [START_REF] Ramos | An alternative point process framework for modeling multivariate extreme values Comm[END_REF], we introduce, in Section 3, a new dependence function similar to the well-known Pickands dependence function which allows us to measure the residual dependence under asymptotic independence. Finally, we estimate this new dependence function and we establish the main asymptotic properties of the estimator. By means of a simulation study, its good performance is highlighted. A discussion on the methodological assumptions and possible extensions of our work ends the paper. All the proofs are postponed to the appendix.

Throughout the paper, the following notations are used. For any arbitrary dimension d and f : X ⊂ R d → R, we set f ∞ = sup x∈X f (x). We denote by ∞ (X) the space of all bounded real-valued functions on X. The symbol " " stands for convergence in distribution of random vectors, but also for weak convergence of bounded real-valued functions in ∞ (X), the difference will be clear from the context.

A test for asymptotic independence

A d-dimensional random vector X = (X 1 , . . . , X d ) follows the law of a multivariate extreme-value distribution if the one-dimensional marginal distributions, G j (x) = Pr(X j ≤ x) for all x ∈ R, j = 1, . . . , d, are Generalized Extreme-Value (GEV) distributions, and the joint distribution takes the form

G(x) = C G 1 (x 1 ), . . . , G d (x d ) , x ∈ R d ,
where C is an extreme-value copula, i.e., Haan and Ferreira, 2006, Ch. 1, 6, for details). Consider the map L : [0, ∞) d → [0, ∞), defined by L(z) := V(1/z) with z = 1/y for y ∈ (0, ∞] d . The function L is known as the stable tail dependence function. As it is a homogeneous function of order one, i.e. L(az) = aL(z) for all a > 0, we have

C(u) = exp -V (-log u 1 ) -1 , . . . , (-log u d ) -1 , u ∈ (0, 1] d , with V : (0, ∞] d → [0, ∞) (see de
L(z) = (z 1 + • • • + z d ) A(t), z ∈ [0, ∞) d , with t j = z j /(z 1 + • • • + z d ) for j = 2, . . . , d, t 1 = 1 -t 2 -• • • -t d ,
d : v 1 + • • • + v d = 1 .
The function A is well-known as the Pickands dependence function (see [START_REF] Pickands | Multivariate extreme value distributions, with a discussion Bull[END_REF], and is often used to quantify the dependence among the elements of X. Indeed, A satisfies the constraint 1/d ≤ max(t 1 , . A(t) ≤ 1 for all t ∈ S d , with lower and upper bounds corresponding to the complete dependence and independence cases, respectively (see Falk et al., 2010, Ch. 4, for details). Thus, estimating this Pickands dependence function is crucial for analysing multivariate extremes, and it has been an extensively discussed topic in the literature, see [START_REF] Klüppelberg | Bivariate extreme value distributions based on polynomial dependence functions Math[END_REF], [START_REF] Zhang | Nonparametric estimation of the dependence function for a multivariate extreme value distribution[END_REF], [START_REF] Gudendorf | Nonparametric estimation of an extreme-value copula in arbitrary dimensions[END_REF]Segers (2011), Capéraà et al. (1997), [START_REF] Berghaus | Minimum distance estimators of the Pickands dependence function and related tests of multivariate extreme-value dependence[END_REF] or [START_REF] Vettori | A comparison of dependence function estimators in multivariate extremes[END_REF], among others.

2.1.

A slightly modified version of the Pickands dependence estimator proposed by [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF] This estimator is based on the madogram concept, a notion borrowed from geostatistics in order to capture the spatial structure. Starting from independent and identically distributed (i.i.d.) copies X 1 , ..., X n , of X, our estimator is defined as

A n (t) := ν n (t) + c(t) 1 -ν n (t) -c(t) (2.1)
where

ν n (t) := 1 n n i=1         d j=1 G (1) n, j (X i, j ) 1/t j - 1 d d j=1 G (1) n, j (X i, j ) 1/t j         (2.2) c(t) := 1 d d j=1 t j 1 + t j with G (a) n, j (X i, j ) := G n, j (X i, j )        1 + a a 1 n n k=1 G 1/a n, j (X k, j )        -a , j = 1, . . . , d, for a > 0,
and the empirical distribution functions denoted by

G n, j (x) := 1 n n i=1
1l {X i, j ≤x} , j = 1, . . . , d.

By convention, here u 1/0 = 0 for 0 < u < 1. The estimator (2.2) is a slightly modified version of that proposed in [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF], with G (1) n, j in place of G n, j which ensures that the new Pickands estimator A n now satisfies A n (e j ) = 1 for all j = 1, . . . , d, where e j = (0, . . . , 0, 1, 0, . . . , 0) is the jth canonical unit vector (see Appendix A.1). This is a necessary condition that a function needs to satisfy in order to be a valid Pickands dependence function (see e.g. [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF]. Although, as established in Appendix A.1, our modified estimator shares the same asymptotic properties as the estimator discussed in [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF], our modification greatly improves the latter for finite samples.

Construction of our statistical test

Using our estimator for A, we want now to construct a statistical test to check asymptotic independence in dimensions higher than or equal to two. To this aim, we consider the following system of hypotheses

H 0 : A(t) = 1, ∀ t ∈ S d H 1 : A(t) < 1, for some t ∈ S d .
Note that H 0 means that all the components of X are asymptotically independent, whereas under H 1 some elements of X are asymptotically dependent. Assuming that the extreme-value copula C has continuous partial derivatives over the sets {u ∈ [0, 1] d : 0 < u i < 1}, by Theorem 2.4 in [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF] and according to Appendix A.1, under H 0 we have that where A is a centered Gaussian process on [0, 1] d with continuous sample paths and covariance function equal to

√ n A n (t) -1 t∈S d -4 1 0 A(v
Cov(A(u), A(v)) = d j=1 u j ∧ v j - d j=1         u j ∧ v j i j u i v i         + (d -1) d j=1 u j v j .
As a consequence, by the continuous mapping theorem (see e.g. van der Vaart, 2000, Ch. 2.1) it follows that

S n := sup t∈S d √ n A n (t) -1 S := sup t∈S d 4 1 0 A(v t 1 , . . . , v t d ) dv , as n → ∞.
This convergence result can be used as the cornerstone of our test. Denoting by Q S (α), α ∈ (0, 1), the (1α)-quantile function for the distribution of the random variable S , H 0 can be rejected with approximately α% significance level whenever s n , the observed value of S n , exceeds Q S (α). Unfortunately, there is no closed form for the function Q S (α), however an approximation can still be computed with a Monte Carlo simulation as follows.

Note that for any u, v ∈ [0, 1] and t, w ∈ S d , the covariance function of the Gaussian process A in (2.3), evaluated at the indexes

u t , v w ∈ [0, 1] d , is equal to Cov(A(v t 1 , . . . , v t d ), A(u w 1 , . . . , u w d )) = d j=1 (v t j ∧ u w j ) - d j=1 (v t j ∧ u w j )v 1-t j u 1-w j + (d -1)uv.
(2.4) Thus, for any fixed α ∈ (0, 1), an approximation of the quantile Q S (α) can be obtained by adhering to the following four steps:

1. Divide the unit interval (0, 1) and the simplex S d in p and m equally spaced points, where p and m are positive integers. Let v 1 , . . . , v m and t 1 , . . . , t p be the two sequences of points partitioning (0, 1) and S d , respectively. ∈ [0, 1] d , with r = 1, . . . , m and k = 1, . . . , p, on which the process A is simulated. 2. Sample n * realizations

x i (v t 1,1 1 , . . . , v t 1,d 1 ), . . . , x i (v t p,1 m , . . . , v t p,d m ), i = 1, . . . , n * , of a zero-mean Gaussian process at v t k,1 r , . . . , v t k,d
r , for r = 1, . . . , m and k = 1, . . . , p, with a (mp × mp) variance-covariance matrix defined through the covariance function in (2.4). 3. Simulate samples that approximately follow the distribution of the random variable S , the integral and the sup in S being approximated by a sum and the max for sufficiently large values of m and p. This leads to the realizations

s i = max 1≤k≤p 4 1 m m r=1 x i (v t k ,1 r , . . . , v t k ,d r ) , i = 1, . . . , n * .
4. An approximation of the quantile Q S (α), denoted by Q S (α), can then be obtained by computing the sample quantile of the realizations s 1 , . . . , s n * for sufficiently large n * .

Numerical results

We illustrate the performance of our statistical test through a simulation study. Precisely, we estimate some values of the significance level α and the power 1 -β of the test by computing the empirical proportion of simulated samples under the null hypothesis and the alternative hypothesis that rejected the null hypothesis, respectively. For simplicity we focus on the significance levels α = 0.05 and 0.01. The study consists of five experiments.

First experiment: In order to perform the first experiment, in a first step we compute the approximated quantile Q S (α), for a given α, following our algorithm. The quality of approximation relies on the values of the indexes m, p and n * . Clearly, the larger their values are, the more accurate the approximation is. We set n * = 500 000. We consider increasing values of m and p and for each combination we compute Q S . We stop the search of a better value for these indexes when the value of Q S (α) does not increase anymore, up to the second decimal. The calculation of Q S requires a considerable computational effort, therefore we derive its values only for a dimension d = 2, 3, 4 of the vector X.

In a second step, we compute the rejection rates. We focus on the multivariate logistic extreme-value model introduced by [START_REF] Tawn | Modelling multivariate extreme value distributions[END_REF], with dependence parameter ψ ∈ (0, 1], ψ = 1 corresponding to independent components of X, whereas complete dependence can be reached when ψ → 0. We consider 20 equally spaced values of ψ in (0, 1]. For each of them, we simulate n independent observations from a logistic extreme-value distribution with unit Fréchet margins. Then we estimate the Pickands dependence function by (2.1) and we compute s n . We repeat this task 5000 times and we compute the proportion of times that s n > Q S (α). This experiment is repeated for different values of the samples size n and different dimension d of X. The middle part of Table 1 reports the estimated values of the significance levels α in the case where ψ = 1. We see that accurate estimates of α are already obtained with the sample size n = 50, indicating a good performance of our statistical test. Figure 1 displays the estimated powers of the test. In the first and second rows the results obtained with α = 0.05 and α = 0.01 are reported, respectively. The panels from left to right illustrate the results for the dimensions 2, 3 and 4. Once again, the test shows a good performance already with the sample size n = 50. Indeed in the case d = 2 we see that the power of the test reaches 1 with mild dependence levels, i.e. ψ = 0.5. This figure also outlines that the power of the test improves as the dimension of X increases and that, as expected, for any fixed dimension d = 2, 3, 4, it also improves as the sample size increases.

Second experiment: We repeat the second step of the first experiment approximating the test statistic's distribution, under the null hypothesis, via a Monte Carlo approach. Precisely, we simulate n values from d independent univariate Fréchet distributions, then we estimate the Pickands dependence function by (2.1) and we compute s n . We repeat this task 5000 times and we compute the empirical quantile, for a given α, denoted by Q S n (α). The right-hand side of Table 1 reports the values of them for different values of n and d. We see that the empirical quantiles rapidly approach the asymptotic quantiles, as the sample size increases. Already with the sample size n = 50 the two types of quantiles are very close. The third and fourth lines of Figure 1 display the comparison between the estimates of 1 -β obtained with the two types of quantiles, but also the estimates of α since ψ = 1 corresponds to independent components and thus the proportion of rejection reported by the figures represents an approximation of α in that case. Since the conclusions are the same for both values α = 0.01 and α = 0.05, only the latter are reported. It follows that the inferential results obtained with the empirical quantiles are very close to those obtained with the asymptotic quantiles, already with the sample size n = 50.

Third experiment: We repeat the second experiment using the Genest-Rémillard (GR) rank-based statistical test [START_REF] Genest | Test of independence and randomness based on the empirical copula process[END_REF] and our proposed test with the Capéraà-Fougères-Genest (CFG) estimator [START_REF] Capéraà | A nonparametric estimation procedure for bivariate extreme value copulas[END_REF][START_REF] Zhang | Nonparametric estimation of the dependence function for a multivariate extreme value distribution[END_REF] of the Pickands dependence function in place of (2.1). Figure 2 shows the estimated powers obtained with the GR test and our test (with both the CFG and the Madogram-based (2.1) estimators). For brevity, we show the results for α = 0.05 and the sample sizes n = 25 and n = 50. We see that our test always outperforms the GR test, with the best results provided with the CFG estimator in the case of a dimension equal to 2, whereas, in higher dimensions, similar results can be reached with either the Madogram-based estimator or the CFG.

Fourth experiment: We repeat the third experiment by sampling from three alternative distributions. In the first case, we draw samples from a three-dimensional random vector with a pair that follows the logistic extreme-value distribution and where the last variable is independent from the other two. In the second case, we consider a four-dimensional random vector with two pairs that follow the logistic extreme-value distribution and where the components of one pair are independent from each component of the second pair. In the last case, we consider a four-dimensional random vector, where one pair with independent components and each of these independent from the components of the other pair that follows the logistic extreme-value distribution. The results are collected in the third (n = 25) and fourth row (n = 50) of Figure 2. In these cases we see that our test loses power and provides inferential results very similar to those provided by the GR test, however the latter outperforms our test in the case of the largest number of independent variables. Fifth experiment: We consider the multivariate inverted symmetric logistic model (see e.g. [START_REF] Ledford | Modelling dependence within joint tail regions J. Roy[END_REF][START_REF] Wadsworth | Modelling across extremal dependence classes[END_REF], with dependence parameter ψ ∈ (0, 1], ψ = 1 corresponding to exact independence of the components of X, whereas asymptotic dependence is reached as ψ → 0. This time, we consider 10 equally spaced values of ψ in (0, 1]. For each of them, we simulate 366 values (for similarity with annual maxima) from an inverted logistic distribution with exponential margins. Then, we compute the normalized componentwise maxima and we repeat this procedure in order to obtain n normalized maxima from which we estimate the Pickands dependence function using (2.1) and we calculate s n . We repeat this task 5000 times and we compute the proportion of times that s n > Q S n (0.05). This procedure has been done for different values of d and n and the results are summarized in Table 2. With d = 2, the rejection rates are close to 0.05 whenever ψ is larger than 0.5. Otherwise, the rejection rate is greater than 0.05 and it reaches 1 when ψ approaches 0. In these cases, it can be observed that the normalized maxima show quite a strong dependence, which indeed seems that of an asymptotically dependent model rather than an asymptotically independent one. The strength of the dependence is reduced whenever the normalized maxima are computed on sequences larger than 366, resulting in improvements in the performances of our test. The test performance deteriorates as the dimension of X increases.

In conclusion this study highlights the good performance of our statistical test in detecting the exact independence of sample maxima. However, our test is also useful to detect if multivariate data are asymptotically independent as long as there is a weak dependence within the case of asymptotic independence. On the contrary, in the cases of strong dependence within asymptotic independence our test fails in detecting the data as asymptotically independent. Clearly, these are the most naturally difficult cases to detect and more specific tools are required.

imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=2 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α n=25 n=50 n=100 n=200 =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α n=25 n=50 n=100 n=200 =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α n=25 n=50 n=100 n=200 =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=2 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α n=25 n=50 n=100 n=200 =0.01 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α n=25 n=50 n=100 n=200 =0.01 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α n=25 n=50 n=100 n=200 =0.01 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=2, n=25 ψ q q q q q q q q q q q q q q q q q q q q 1 -β α Empirical quantile Asymptotic quantile =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3, n=25 ψ q q q q q q q q q q q q q q q q q q q q 1 -β α Empirical quantile Asymptotic quantile =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=25 ψ q q q q q q q q q q q q q q q q q q q q 1 -β α Empirical quantile Asymptotic quantile =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=2, n=50 ψ q q q q q q q q q q q q q q q q q q q q 1 -β α Empirical quantile Asymptotic quantile =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3, n=50 ψ q q q q q q q q q q q q q q q q q q q q 1 -β α Empirical quantile Asymptotic quantile =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=50 ψ q q q q q q q q q q q q q q q q q q q q 1 -β α Empirical quantile Asymptotic quantile =0.05 imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=2, n=25 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3, n=25 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=25 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=2, n=50 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3, n=50 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=50 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3, n=25 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=25 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=25 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=3, n=50 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=50 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05 q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 d=4, n=50 ψ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 -β α Madogram CFG GR =0.05

Asymptotic independence for componentwise maxima

Being able to test asymptotic independence versus asymptotic dependence is obviously important. However, since asymptotic independence often arises in applications, it is also crucial to develop some general models that accommodate both situations. In this section, we consider the framework of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] (see also [START_REF] Ledford | Modelling dependence within joint tail regions J. Roy[END_REF]. More precisely, if Y is a d-dimensional random vector with common unit Fréchet margins, i.e. Pr(Y j ≤ y) = e -1/y for every y > 0 and 1 ≤ j ≤ d, this theory relies on the joint survival function of Y which is assumed to be multivariate regularly varying with index -1/η, where η ∈ (0, 1], i. for all y ∈ (0, ∞] d . The function g here is homogeneous of order 0, i.e. such that g(a y 1 , . . . , a y d ) = g(y 1 , . . . , y d ) for any a > 0. This framework implies that the conditional joint survival function admits the following limit representation for every y ≥ 1, the vector of ones,

lim u→∞ Pr(Y > uy|Y > u1) = lim u→∞ Pr(Y > uy) Pr(Y > u1) = η S d d j=1 w j y j 1/η dH η (w), (3.1)
where H η is a non-negative measure satisfying the condition

η S d d j=1 w 1/η j dH η (w) = 1. (3.2)
This measure H η is a particular case of the hidden angular measure introduced by Resnick (2002) (see also [START_REF] Maulik | Characterizations and examples of hidden regular variation[END_REF] when η < 1 and it is a rescaled version of the classical angular measure when η = 1, see [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] for details. According to [START_REF] Ramos | An alternative point process framework for modeling multivariate extreme values Comm[END_REF] we assume that H η is a finite measure. We recall that η is the so-called coefficient of tail dependence, which measures the level of dependence within the asymptotic independence framework. Specifically, η = 1 corresponds to the case of asymptotic dependence, whereas η < 1 corresponds to the case of asymptotic independence. More precisely, when the coefficient η falls in the following sets: (1/d, 1), {1/d} or (0, 1/d), then we say that among the variables there is a positive association, independence or negative association, respectively, within asymptotic independence (see e.g. [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF].

A η-Pickands dependence function

Consider now, n i.i.d. copies Y 1 , . . . , Y n of Y and for a small ε > 0, define M n,ε = (M n,1,ε , . . . , M n,d,ε ) as the vector of componentwise maxima, precisely

M n, j,ε = i∈I n (ε) Y i, j , j = 1, . . . , d, with I n (ε) := {1 ≤ i ≤ n : Y i > 1ε}.
Let b n be a sequence of normalizing constants defined by the equation nPr(Y > b n ) = 1. Then, differently from the classical theory (e.g. de Haan and Ferreira, 2006, Ch. 6), here the limiting distribution for the normalized vector of componentwise maxima M n,ε is obtained as the law of a multivariate η-extreme-value distribution, if the one-dimensional marginal distributions are G η, j (y) = exp(-σ η, j y -1/η ), for all y > 0, j = 1, . . . , d, and the joint distribution takes the form

G η (y) := lim ε→0 lim n→∞ Pr(M n,b n ε ≤ b n y), y ∈ (0, ∞] d , ( 3
G η (y) = C η G η,1 (y 1 ), . . . , G η,d (y d ) , y ∈ (0, ∞] d , (3.4)
where C η is an η-extreme-value copula, i.e.

C η (u) = exp -V η σ η,1 -log u 1 η , . . . , σ η,d -log u d η u ∈ (0, 1] d with V η : (0, ∞] d → [0, ∞) a homogeneous function of order -1/η and σ η, j := V η (∞, . . . , ∞, 1, ∞, . . . , ∞) = η S d w 1/η j dH η (w). (3.5) Introduce now L η (z) := V η ((σ σ σ η /z) η
), for all z = σ σ σ η /y 1/η . This function is called the η-stable tail dependence function and using the homogeneity property, it can be rewritten as

L η (z) = (z 1 + • • • + z d )A η (t), z ∈ [0, ∞) d ,
where

t j = z j /(z 1 + • • • + z d ) for j = 2, . . . , d, t 1 = 1 -t 2 -• • • -t d .
Here, the function A η is called the η-Pickands dependence function and it satisfies the following properties.

Proposition 3.1. The η-Pickands dependence function A η satisfies:

1. For all η ∈ (0, 1], A η (e j ) = 1, j = 1, ..., d; 2. A 1 (t) = A(t), for all t ∈ S d ; 3. For every η ∈ (0, 1] and t ∈ S d , 1/d ≤ max(t 1 , . . . , t d ) ≤ A η (t) ≤ 1.

4. A η (t) is convex, i.e. A η (at 1 + (1a)t 2 ) ≤ aA η (t 1 ) + (1a)A η (t 2 ), for all a ∈ [0, 1] and t 1 , t 2 ∈ S d .

Similarly to the classical literature, a η-madogram function can be defined as the expected distance between the maximum and the mean of the variables G 1/ηt 1 η,1 (Z 1 ), . . . , G 1/ηt d η,d (Z d ), that is,

ν η (t) = E         d j=1 G 1/ηt j η, j (Z j ) - 1 d d j=1 G 1/ηt j η, j (Z j )         .
This function can also be linked to the η-Pickands dependence function as follows.

Proposition 3.2. Any random vector Z with a η-extreme-value distribution admits a η-Pickands dependence function A η which satisfies

A η (t) = 1 η ν η (t) + c η (t) 1 -ν η (t) -c η (t) (3.6)
for all t ∈ S d , where

c η (t) = 1 d d j=1 t j t j + 1/η . (3.7)
This η-Pickands dependence function can be used to represent the level of dependence among the elements of Z, and thus in the next section, we estimate this function and derive the main asymptotic properties of the estimator.

imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

An estimator of the η-Pickands dependence function

Let Z 1 , . . . , Z n be i.i.d. copies of Z with distribution G η and define

H n (y) = 1 n n i=1 1l {Z i ≤y} , y ∈ (0, ∞] d
and its associated empirical process

H n (y) = √ n(H n (y) -G η (y)), y ∈ (0, ∞] d .
In order to estimate the η-Pickands dependence function we first assume that we have at our disposal an estimator η n for η satisfying the condition:

Condition 1. Let η n be an estimator of η satisfying:

(i) η n → η a.s. as n → ∞;

(ii) One of the following holds true

(a) √ n( η n -η) = n -1/2 n i=1 ρ(Z i ) + o p (1), where ρ : (0, ∞] d → R is a measurable function such that Eρ(Z) = 0 and Eρ 2 (Z) < ∞; (b) √ n( η n -η) = χ(H n ) + o p (1), where χ : ∞ ((0, ∞] d ) → R is a bounded linear functional.
In the spirit of (2.1) in Section 2, we propose the following estimator for A η :

A η n ,n (t) := 1 η n ν η n ,n (t) + c η n ,n (t) 1 -ν η n ,n (t) -c η n ,n (t)
where

ν η n ,n (t) := 1 n n i=1         d j=1 H ( η n ) n, j (Z i, j ) 1/ η n t j - 1 d d j=1 H ( η n ) n, j Z i, j 1/ η n t j         c η n ,n (t) := 1 n d n i=1 d j=1 H ( η n ) n, j Z i, j 1/ η n t j (3.8)
with

H (a) n, j (Z i, j ) = H n, j (Z i, j )        1 + a a 1 n n k=1 H 1/a n, j (Z k, j )        -a , j = 1, . . . , d, for a > 0,
and the empirical distribution functions denoted by

H n, j (x) := 1 n n i=1
1l {Z i, j ≤x} , j = 1, . . . , d.

Note that (3.8) comes from the fact that c η defined in Proposition 3.2 can be viewed as

c η (t) = E         1 d d j=1 G η, j (Z j ) 1/ηt j        
and thus in (3.8) we use the empirical counterpart. Another option would have been to replace η by an estimator in (3.7).

We are now able to state our main result on the convergence of a rescaled version of A η n ,n .

imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

Theorem 3.1. Let Z 1 , . . . , Z n be i.i.d. copies of Z with distribution G η , and A η n ,n be our proposed estimator for A η .

Under Condition 1(i), we have A η n ,n -A η ∞ → 0 a.s. as n → ∞.

(3.9)

Under Conditions 1(ii), we have in ∞ (S d ), as n → ∞, √ n( A η n ,n (t) -A η (t)) t∈S d - (1 + ηA η (t)) 2 η 1 0 A η (v ηt 1 , . . . , v ηt d ) dv t∈S d , (3.10)
where A η is a stochastic process defined as

A η (u) := B η (u) - d j=1 Ċη, j (u)B η (1, . . . , 1, u j , 1, . . . , 1), u ∈ [0, 1] d ,
with C η an η-extreme-value copula such that its partial derivative Ċη, j (u) := ∂C η /∂u j (u) exists and is continuous on {u ∈ [0, 1] d : 0 < u j < 1}, for all j = 1, . . . , d, and B η a C η -Brownian bridge, i.e. a zero-mean Gaussian process on [0, 1] d with continuous sample paths and covariance function equal to

Cov B η (u), B η (v) = C η (u ∧ v) -C η (u)C η (v), u , v ∈ [0, 1] d .

Examples of estimators satisfying Condition 1

Our η-Pickands dependence function requires an estimator of η which satisfies Condition 1. Below, two examples of such estimators are proposed.

Example 1. Let Z * = max(Z 1 , . . . , Z d ), where Z follows the distribution (3.4). Then, for any y > 0, the distribution of Z * is G η (y) := G η (y, . . . , y). This distribution can be seen as a two-parameter Fréchet family of distributions. Let η n be the Maximum Likelihood (ML) estimator. By Propositions 3.1 and 3.3 in [START_REF] Bücher | On the maximum likelihood estimator for the generalized extreme-value distribution[END_REF], it follows that the ML estimator satisfies Conditions 1(i) and 1(ii)(a). and on u ∈ (0, 1) the two functions

γ(u) := µ 1,2 u(-log u) -µ 1,1 u(-log u) 2 ϕ(u) := 1 η V η η (1, . . . , 1) u (-log u) 1+η . Then, √ n( η n -η) = - 2 µ 2 1,1 1 0 H n (Q η (u), . . . , Q η (u)) γ(u) ϕ(u) du + o(1) a.s.. (3.11)
Consequently, as n → ∞

η n → η a.s. √ n( η n -η) - 2 µ 2 1,1 1 0 H(Q η (u), . . . , Q η (u)) γ(u) ϕ(u) du
where H is a tight centered Gaussian process on (0, ∞] d , with covariance function

Cov(H(z), H(y)) = G η (z ∧ y) -G η (z)G η (y), z, y ∈ (0, ∞] d .
imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

Simulation

The performance of our estimator A η n ,n is illustrated in a simulation study with two different experiments.

First experiment: We consider the bivariate η-asymmetric logistic dependence model introduced by Ramos and [START_REF] Ramos | An alternative point process framework for modeling multivariate extreme values Comm[END_REF]. Such a dependence structure is characterised by the following features. The hidden spectral measure H η has density

h η (w) = η -ψ ψη 2 N        (ρw) -1/ψ + 1 -w -1/ψ        ψ/η-2
{w(1w)} -(1+1/ψ) , w ∈ (0, 1), (3.12)

where N := -1/η + 1/η -( -1/η + 1/η ) ψ/η and ψ ∈ (0, 1], > 0, η ∈ (0, 1] are dependence parameters. This satisfies the condition (3.2), i.e. in case d = 2,

η -1 = 1/2 0 ω 1/η dH η (w) + 1 1/2 (1 -w) 1/η dH η (w).
The associated limiting distribution in (3.3) takes the form (3.13) where the degenerate case arises when H η is infinite. In the sequel we focus on the case ψ < η and for simplicity we consider = 1. Distribution (3.13) is the attractor for the distribution of normalized componentwise maxima obtained from a random vector whose survival function is ), with u being a high threshold and λ the joint threshold exceedance probability (see Ramos and Ledford, 2009 for details). The survival function (3.14) satisfies (3.1) and it is an asymptotically independent joint probability model for any η ∈ (0, 1), where the strength of the dependence, within asymptotic independence, increases for decreasing values of the parameter ψ. We call (3.13) and (3.14) the η-asymmetric logistic distribution and survival function, respectively. We simulate n values from the η-asymmetric logistic distribution and we estimate the η-Pickands dependence function with A η n ,n . We repeat these steps 1000 times and we compute a Monte Carlo approximation of the Mean Integrated Squared Error (MISE), i.e.,

G η (y 1 , y 2 ) =          exp -N -1 ( y 1 ) -1/ψ + y 2 -1/ψ ψ/η for ψ < η 0 for ψ ≥ η,
Pr(Y 1 > y 1 , Y 2 > y 2 ) = λu 1/η N          ( y 1 ) -1/η + y 2 -1/η -        ( y 1 ) -1/ψ + y 2 -1/ψ        ψ/η          (3.14) where (y 1 , y 2 ) ∈ [u, ∞) × [u, ∞

MISE(

A η n ,n , A η ) = E S d A η n ,n (t) -A η (t) 2 dt .
This study is done for different values of the sample size n and different values of the dependence parameter ψ. The results are summarized in Table 3. For each value of ψ, between the second and the fifth column the mean of the estimates for η obtained with the GPWM (first row) and ML (second row) estimator are reported for increasing sample size. In parentheses is the standard deviation. Between the sixth and ninth columns the approximated MISE is reported. Accurate estimates are obtained with all the dependence levels. GPWM and ML estimators provide similar results, although those of the former seem slightly better. According to the MISE, the better performances are obtained with stronger dependence strengths. For every dependence level the accuracy of estimates increases with increasing sample size.

Second experiment: We show the performance of the estimator A η n ,n under a more realistic scenario. We simulate n × 366 independent observations from a distribution whose survival function is given in (3.14). To do this we use the algorithm described in Theorem 1.1 and Appendix B of [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF]. The simulation procedure in [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF] relies on the condition Pr(Y > y) = ληu 1/η (0,1) min (w/y 1 , (1w)/y 2 ) 1/η dH η (w), y > u1, imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018 where u and λ are as in the first experiment. This condition implies that for every y > u,

1 -Pr(Y 1 > y, Y 2 ≤ u) = exp(-1/y) + λη (0,1) min w y/u , 1 -w 1/η dH η (w), 1 -Pr(Y 1 ≤ u, Y 2 > y) = exp(-1/y) + λη (0,1) min w, 1 -w y/u 1/η dH η (w).
The values of u and λ must be selected in such a way that both functions of y above are monotonically increasing.

When the density of H η is given by (3.12) with = 1 and η = 0.7, the monotonicity conditions are satisfied for every ψ < η by setting u = 10 and λ = 1exp(-0.1) -0.02. With simulated data we compute b 366 , that is the empirical (1 -1/366)-quantile of the minimum between pairs of all observations. For each block of 366 observations we compute the componentwise maxima using ε = Q(0.07)/ b 366 , where Q(0.07) is the ninety-third percentile of a unit Fréchet distribution, i.e. by retaining only the pairs that are both greater than ε b 366 . We standardize the maxima by dividing them by b 366 . With the n normalized maxima we estimate the η-Pickands dependence function by A η n ,n . We repeat these steps 1000 times and we compute an approximation of the MISE. Table 4 collects the results.

We see that the estimates of η and A η are similar to those obtained in Table 3, indicating a good performance of our estimator. We mention that in each block of 366 observations the componentwise maxima are computed, imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018 after the truncation, on average on approximately 17 pairs. Although maxima are obtained with a small number of observations, the estimation results suggest that they are enough to obtain accurate estimates. Estimates are less accurate for ψ = 0.6 for the following reason. The simulation method for maxima produces observations that are approximately drawn from the non-degenerate distribution G η in (3.13), provided that ψ < η, since G η is a degenerate distribution for ψ ≥ η. Furthermore, for this example it can be empirically verified that G η provides a very accurate approximation for the distribution of the simulated maxima when ψ < 0.6. Instead, whenever ψ is close to η (a case that resembles the degenerate case), e.g. 0.6 ≤ ψ < 0.7 = η, the quality of the approximation deteriorates. In this case the mismatch between G η and the distribution of simulated maxima is no longer negligible, thus affecting the estimation results.

Finally, note that the asymptotic properties of our estimator established in Theorem 3.1 are no longer valid in this experiment, although our estimator still performs well. Indeed, they should be re-established under the assumption that the data belong to the domain of attraction of G η . In that case the proofs are much more technical and thus are outside the scope of the present paper.

Discussion

The framework for modelling the dependence within asymptotic independence, based on componentwise maxima, relies on the assumption that H η is a finite measure (see Section 3). Before applying our estimation method, it is desirable to check somehow whether such an assumption is not violated by the available data. To this end, we propose a diagnostic tool. We motivate it on the basis of the following discussion. For simplicity we focus on the bivariate case although our proposal is easily extendable to higher dimensions. Let Y be a two-dimensional random vector defined as in Section 3. Define

F j (y) := lim s→∞ Pr(Y j > sy |Y > s1), j = 1, 2, s > 0, y > 1, = η y -1/η y/(1+y) 0 w 1/η dH η (w) + 1 y/(1+y) (1 -w) 1/η dH η (w) ≤ 2ηMy -1/η ,
where M = H η ((0, 1)) < ∞, and

F min (y) := lim s→∞ Pr(min(Y 1 , Y 2 ) > sy |Y > s1).
Then, it follows that 1 ≤ 2ηM min y -1/η / F 1 (y), y -1/η / F 2 (y) = 2ηM min F min (y)/ F 1 (y), F min (y)/ F 2 (y) .

Consequently marginal survival functions F j , j = 1, 2, heavier than F min , i.e. F min (y)/ F j (y) → 0 as y → ∞, j = 1, 2, provide empirical evidence against the hypothesis that H η is finite. On the contrary, evidence in favor of a finite H η is provided by the conditions F min (y)/ F j (y) → c j as y → ∞, where c j , j = 1, 2, are positive constants. On this basis, we suggest implementing the plot r j (y) := F min (y)

F j (y) , 1 ≤ y ≤ m * /
Pr(Y 1 > sy 1 , Y 2 > sy 2 |Y > s1) ≈ η 1 0 min w y 1 , 1 -w y 2 1/η dH η (w) = -log G η,1 (y 1 ) + log G η,2 (y 2 ) 1 -A η log G η,1 (y 1 ) log G η,1 (y 1 ) + log G η,2 (y 2 ) , log G η,2 (y 2 ) log G η,1 (y 1 ) + log G η,2 (y 2 )
holds for a large threshold s and y 1 , y 2 > 1. Set s = b n , where b n is the empirical (1 -1/n)-quantile of the sequence min(y i,1 , y i,1 ), i = 1, . . . , n, with y i,1 , y i,2 that are independent realizations of Y, see (3.3). Then, the above probability can be approximated by log H n,1 (y 1 ) + log H n,2 (y 2 ) 1 -A η n ,n log H n,1 (y 1 ) log H n,1 (y 1 ) + log H n,2 (y 2 )

, log H n,2 (y 2 ) log H n,1 (y 1 ) + log H n,2 (y 2 ) , (4.1)

where H n, j , j = 1, 2, are the empirical distribution functions, see Section 3.2. We illustrate the extrapolation of the probability of high thresholds exceedances with two examples. We simulate 500 × 366 independent realizations from two distributions with an η-asymmetric logistic survival function and parameters ψ = 0.1, ψ = 0.4, respectively, while u = 10, λ = 1e -0.1 -0.02, = 1, η = 0.7 are the same for both the cases. Then, we obtain the sample of maxima, using b 366 , ε = Q(0.07)/ b 366 (see Section 3.1 and the second experiment of Section 3.4 for details) and we estimate the Pickands dependence function with A η 500 ,500 , where η 500 is the GPWM estimator of η. For y 1 , y 2 ∈ [1, 10], we extrapolate the probability of joint high thresholds exceedances by applying (4.1). Figure 4 displays the estimated probabilities for the two models. The left and right panels report the results for the cases ψ = 0.1 and ψ = 0.4, respectively. To go further with this idea, a topic of interest would be to establish the asymptotic properties of the estimator defined in (4.1). This is outside the scope of the present paper but it will lead to further investigations.

imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

       2 n n k=1 G n, j (X k, j )        -1 = n n + 1 G n, j (X i, j ).
Then, with simple adjustments of the proof of Theorem 2.4 in [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF], the weak convergence of A n and its almost sure consistency follow.

A.2. Proof of Proposition 3.1

Our definition of L η combining with (6.3) in [START_REF] Ramos | An alternative point process framework for modeling multivariate extreme values Comm[END_REF] entails

A η (t) = η S d max        t 1 w 1/η 1 σ η,1 , . . . , t d w 1/η d σ η,d        dH η (w), t ∈ S d .
Then, Property 1 follows by the definition of σ η, j given in (3.5). When η = 1, according to Section 3, we have

lim n→∞ Pr(Y > nx) Pr(Y > n1) = S d d j=1
w j x j dH 1 (w). Now, this limit can also be rephrased with the classical theory (see e.g. de Haan and Ferreira, 2006, Ch. 6), where

lim n→∞ Pr(Y > nx) Pr(Y > n1) = d S d d j=1 w j x j dH(w) R(1, . . . , 1)
, with H and R defined in pages 218 and 225 in de [START_REF] De Haan | Extreme Value Theory: An Introduction Springer[END_REF]. Therefore, Property 2 follows from the relations d -1 R(1, . . . , 1)dH 1 (w) = dH(w), w ∈ S d and σ 1, j = 1/R(1, . . . , 1), j = 1, . . . , d.

For every t ∈ S d we have from which the upper bound in Property 3 follows. To derive the lower bound, it is sufficient to remark that for every t ∈ S d , we have

η S d max        t 1 w 1/η 1 σ η,1 , . . . , t d w 1/η d σ η,d        dH η (w) ≤ η S d d j=1         t j w 1/η j σ η, j         dH η (w) = 1, imsart-
η S d max        t 1 w 1/η 1 σ η,1 , . . . , t d w 1/η
We will prove Theorem 3.1 with H ( η n ) n, j in ν η n ,n and c η n ,n replaced by H n, j . Indeed, this slight modification has no impact on the convergences (3.9) and (3.10) since

H ( η n ) n, j (Z i, j ) = H n, j (Z i, j ) 1 + 1 + η n η n O 1 n -η n =: H n, j (Z i, j ) e -1 n ,
and the terms in (A.2) and (A.3) can be slightly changed by replacing in the integrals v • t j by v • t j e n , j = 1, . . . , d, without any impact. In view of this remark, we pursue the proof of Theorem 3.1 with M(•, t) and M n (•, t) defined in (A.2) and (A.3) without taking care of the adjustment with e n . We start to prove (3.10). To this aim, note that from (A.1) we have

M(η, t) = ηA η (t) 1 + ηA η (t)
and thus the following decomposition holds

√ n( A η n ,n (t) -A η (t)) = √ n       1 η n M n ( η n , t) 1 -M n ( η n , t) - 1 η M(η, t) 1 -M(η, t)       = √ n η n       M n ( η n , t) 1 -M n ( η n , t) - M(η, t) 1 -M(η, t)       + M(η, t) 1 -M(η, t) √ n 1 η n - 1 η =: L n (t) + R n (t),
for all t ∈ S d . We derive a tractable expression for L n by means of the following three results.

Lemma A.1. We have the following decomposition

√ n( M n ( η n , t) -M(η, t)) = √ n( M n (η, t) -M(η, t)) + √ n(M( η n , t) -M(η, t)) + o p (1).
Proof. The proof uses arguments from van der [START_REF] Van Der Vaart | Empirical processes indexed by estimated functions IMS Lecture Notes[END_REF]. Since

√ n( M n ( η n , t) -M(η, t)) = √ n( M n ( η n , t) -M( η n , t)) - √ n( M n (η, t) -M(η, t)) + √ n( M n (η, t) -M(η, t)) + √ n(M( η n , t) -M(η, t)), it remains to show that √ n( M n ( η n , t) -M( η n , t)) - √ n( M n (η, t) -M(η, t)) ∞ = o p (1). (A.4)
By Condition 1(ii) we have that √ n( η n -η) is asymptotically tight. Thus, for every ε > 0, there exists a compact set

K ≡ K ε ⊆ R such that lim inf n→∞ Pr( √ n( η n -η) ∈ K) > 1 -ε.
Furthermore, by the compactness of K, there exist δ > 0, p := p(δ) ∈ N and {h 1 , . . . ,

h p } ⊆ K such that K ⊆ ∪ 1≤s≤p (h s -δ, h s + δ). Therefore, { √ n( η n -η) ∈ K} ⊆        √ n( η n -η) ⊆ p s=1 (h s -δ, h s + δ)        = p s=1 η n ∈ η + n -1/2 (h s -δ), η + n -1/2 (h s + δ) .
imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

Consequently, it follows that, with probability at least 1

-ε, √ n( M n ( η n , t) -M( η n , t)) - √ n( M n (η, t) -M(η, t)) ∞ ≤ sup t∈S d max 1≤s≤p sup |h-h s |<δ | √ n( M n (η n,h , t) -M(η n,h , t)) - √ n( M n (η, t) -M(η, t))| ≤ sup t∈S d max 1≤s≤p | √ n( M n (η n,h s , t) -M(η n,h s , t)) - √ n( M n (η, t) -M(η, t))| + sup t∈S d max 1≤s≤p sup |h-h s |<δ | √ n( M n (η n,h s , t) -M(η n,h s , t)) - √ n( M n (η n,h , t) -M(η n,h , t))| =: I n,1 + I n,2
where η n,• := η + n -1/2 •. Thus to show (A.4) it is sufficient to prove that both I n,1 and I n,2 tends to 0 in probability, as n → ∞. Using (A.2) and (A.3) we obtain

I n,1 = sup t∈S d max 1≤s≤p 1 0 C n (v η n,hs t 1 , . . . , v η n,hs t d ) -C n (v ηt 1 , . . . , v ηt d ) dv ≤ sup t∈S d max 1≤s≤p sup v∈(0,1) C n (v η n,hs t 1 , . . . , v η n,hs t d ) -C n (v ηt 1 , . . . , v ηt d )
and

I n,2 = sup t∈S d max 1≤s≤p sup |h-h s |<δ 1 0 C n (v η n,hs t 1 , . . . , v η n,hs t d ) -C n (v η n,h t 1 , . . . , v η n,h t d ) dv ≤ sup t∈S d max 1≤s≤p sup |h-h s |<δ sup v∈(0,1) C n (v η n,hs t 1 , . . . , v η n,hs t d ) -C n (v η n,h t 1 , . . . , v η n,h t d ) .
Now, for every v ∈ (0, 1) and small > 0, the map ϕ :

(0, 1) → ∞ ([η -, η + ]) : v → ϕ(v), defined by (ϕ(v))(x) = v x
, induces continuously differentiable functions on [η-, η+ ] for every v ∈ (0, 1). The first derivative of such functions is ( φ(v))(x) = v x log v, which is bounded above by ξ v = v η-| log v|. Therefore, (ϕ(v))(x) is a Lipschitz function and it satisfies the condition

|(ϕ(v))(x) -(ϕ(v))(y)| ≤ ξ v |x -y|, ∀ x, y ∈ [η -, η + ].
Furthermore, there exists a positive constant ξ such that sup v∈(0,1) ξ v < ξ, and thus for n sufficiently large ensuring that η n,h , η n,h s ∈ [η -, η + ], we have:

|v η n,hs t j -v ηt j | ≤ ξ |η -η n,h s | = ξ n -1/2 |h s | → 0 |v η n,hs t j -v η n,h t j | ≤ ξ |η n,h s -η n,h | = ξ n -1/2 |h s -h| ≤ ξ δn -1/2 → 0, as n → ∞,
|v η n,hs t j -v η n,h t j | → 0, n → ∞. (A.6)
Since the first partial derivative of C η exists and is continuous on {u ∈ [0, 1] d : 0 < u j < 1} for all j = 1, . . . , d, [START_REF] Fermanian | Weak convergence of empirical copula processes[END_REF][START_REF] Segers | Asymptotics of empirical copula processes under non-restrictive smoothness assumptions[END_REF]. Therefore the sequence C n is asymptotically uniformly equicontinuous in probability (see Theorem 1.5.7 in [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF]. Combining this result with (A.5) and (A.6) we can conclude that I n,1 and I n,2 tends to 0 in probability, as n → ∞. Therefore (A.4) is established and thus Lemma A.1 follows.

C n A η in ∞ ([0, 1] d ) as n → ∞ (see e.g.
Lemma A.2. We have Proof. Let

√ n(M( η n , t) -M(η, t)) = A η (t) (ηA η (t) + 1) 2 √ n( η n -η) + o p (1
ϕ : ((0, ∞), | • |) → ( ∞ (S d ), • ∞ ) : a → M(a, •)
be the map defined by

M(a, •) = a A η (•) 1 + a A η (•)
.

Its Hadamard derivative at η ∈ (0, 1] is

h → ( φη (h)) = h A η (ηA η + 1) 2 .
Indeed, for every n ↓ 0 and

h n → h ∈ (0, ∞), as n → ∞, such that η + n h n ∈ (0, ∞), we have lim n→∞ sup t∈S d (ϕ(η + n h n ))(t) -(ϕ(η))(t) n -( φη (h))(t) = lim n→∞ sup t∈S d 1 n (η + n h n )A η (t) (η + n h n )A η (t) + 1 - ηA η (t) ηA η (t) + 1 - hA η (t) (ηA η (t) + 1) 2 = lim n→∞ sup t∈S d A η (t) ηA η (t) + 1 h n (η + n h n )A η (t) + 1 - h ηA η (t) + 1 ≤ lim n→∞ d 2 |h n -h| + |hh n | n (d + η)(d + η + n h n ) = 0.
Lemma A.2 now follows from Theorem 20.8 in [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF] and under our Condition 1(ii).

Lemma A.3. We have

√ n       M n ( η n , t) 1 -M n ( η n , t) - M(η, t) 1 -M(η, t)       = (1 + ηA η (t)) 2 √ n( M n ( η n , t) -M(η, t)) + o p (1).
Proof. The proof of this lemma is based on an application of the functional delta method after proving the Hadamard differentiability of the map ϕ :

∞ (S d ) → ∞ (S d ) : f → f /(1 -f ), with f in ∞ (S d
), and the existence of the weak limit of

√ n(M( η n , •) -M(η, •)) in ∞ (S d ). First, we start showing that the Hadamard derivative of ϕ at M := M(η, •) is h → ( φM (h)) = h (1 -M) 2 , with h in ∞ (S d ). Indeed, for every sequence n ↓ 0 and h n → h as n → ∞, such that M + n h n in ∞ (S d ), we have lim n→∞ sup t∈S d (ϕ(M + n h n ))(t) -(ϕ(M))(t) n -( φM (h))(t) = lim n→∞ sup t∈S d 1 n M(η, t) + n h n (t) 1 -M(η, t) -n h n (t) - M(η, t) 1 -M(η, t) - h(t) (1 -M(η, t)) 2 = lim n→∞ sup t∈S d (1 + ηA η (t)) 2 h n (t) -h(t) + h(t) n h n (t)(1 + ηA η (t)) 1 -n h n (t)(1 + ηA η (t)) ≤ lim n→∞ (1 + η) 2 h n -h ∞ + n h n h ∞ (1 + η) 1 -n h n ∞ (1 + η) = 0.
Then, combining Lemmas A. (Z i, j )} dv,

W i,t = A η (t) (1 + ηA η (t)) 2 ρ(Z i ).
Note that the new expression for T n,1 is obtained by applying Fubini's theorem. The pair (T n,1 , T n,2 ) is asymptotically tight and so to show that its weak limit exists, it remains to prove that all its finite dimensional distributions converge. This can be done by applying the central limit theorem since, for all k = 1, 2, . . ., the i.i.d. random vectors W i,t 1 , . . . , W i,t k , W i,t 1 , . . . , W i,t k , i = 1, . . . , n, have finite second order moments under the assumptions of our Theorem 3.1 (see Nelsen, 2006, Theorem 2.2.7). This achieves the proof of Lemma A.3.

We come back now to the proof of Theorem 3.1. Combining the three previous lemmas with the definition of M(η, t), we have

L n + R n = (1 + ηA η (t)) 2 η n √ n( M n (η, t) -M(η, t)) + A η (t) η n √ n( η n -η) + ηA η (t) √ n 1 η n - 1 η + o p (1) = (1 + ηA η (t)) 2 η n √ n( M n (η, t) -M(η, t)) + o p (1) = - (1 + ηA η (t)) 2 η 1 0 C n (v ηt 1 , . . . , v ηt d )dv + o p (1).
As in the proof of Lemma A.1, using again the convergence C n A η in ∞ ([0, 1] d ) as n → ∞, (3.10) follows from the continuous mapping theorem and Slutsky's lemma. It remains now to prove (3.9). Note that

A η n ,n -A η ∞ = sup t∈S d 1 η n M n ( η n , t) 1 -M n ( η n , t) - 1 η M(η, t) 1 -M(η, t) = sup
≤ | η n -η| 1 -M(η, t) ∞ + η n M(η, t) -M n ( η n , t) ∞ + | η n -η| ≤ 2| η n -η| + η n M( η n , t) -M n ( η n , t) ∞ + η n M(η, t) -M( η n , t) ∞ ≤ 2| η n -η| + η n C n -C η ∞ + η n k 1 0 v η n t 1 -v ηt 1 , . . . , v η n t d -v ηt d ∞ dv.
Each term on the right-hand side of this inequality tend to 0 a.s. under our assumptions and according to similar arguments to those used in Lemma A.1 for the last term. Thus (3.9) is established and the proof of Theorem 3.1 is now complete.

A.5. Proof of Theorem 3.2

According to [START_REF] Guillou | Madogram and asymptotic independence among maxima[END_REF], η can be rewritten as

η = 2 1 - µ 1,2 µ 1,1 .
A natural estimator can thus be obtained by replacing Q η (u) by the empirical version G ← n (u) where G n (u) := G n (u, . . . , u). This entails Consequently, we can decompose the left-hand side of (3.11) as

η n = 2 1 - µ 1,2 µ 1,
√ n( η n -η) = 2 √ n µ 1,2 µ 1,1 - µ 1,2 µ 1,1 = 2 1 0 Q n (u)γ(u)du n -1/2 µ 1,1 1 0 Q n (u)u(-log u)du + µ 2 1,1 =: 2 N n D n with Q n (u) := √ n(Q n (u) -Q η (u)).
We start to study the numerator N n . To this aim, we define the empirical and quantile processes: We can easily check that G η satisfies the conditions of Theorem 3 in [START_REF] Csörgő | Strong approximations of the quantile process[END_REF] Assembling N n and D n , we deduce that

√ n( η n -η) = - 2 µ 2 1,1 1 0 H n (Q η (u), . . . , Q η (u)) γ(u) ϕ(u) du + o(1) a.s.,
where we used the fact that H n (u) = H n (Q η (u), . . . , Q η (u)). Thus (3.11) is established. The other statements of the theorem are direct consequences.

  and A is the restriction of L into the d-dimensional unit simplex, S d := (v 1 , . . . , v d ) ∈ [0, 1]

Fig 1 .

 1 Fig 1. Estimated power functions. Points report the empirical proportion of simulated samples under H 1 that rejected H 0 as a function of ψ.

Fig 2 .

 2 Fig 2. Estimated power functions. Points are the empirical proportion of simulated samples under H 1 that rejected H 0 as a function of ψ.

  e. Pr(Y > y) = τ(y)(y 1 • • • • • y d ) -1/dη with τ a slowly varying function satisfying lim r→∞ τ(r y 1 , . . . , r y d ) τ(r, . . . , r) = g(y)

Example 2 .

 2 Let η n be the Generalized Probability Weighted Moment (GPWM) estimator of η introduced by Guillou et al. (2014). The next theorem shows that the GPWM estimator admits a stochastic representation implying that Condition 1(ii)(b) is satisfied. The almost sure consistency of η n is a direct consequence. Theorem 3.2. Let η n be the GPWM estimator. For a, b two integers and Q η (u) := G ← η (u), introduce the parameter µ a,b := 1 0 Q η (u)u a (-log u) b du

  751(0.125) 0.766(0.094) 0.775(0.070) 0.781(0.050) 0.7105 0.4298 0.3018 0.2333 0.898(0.201) 0.849(0.121) 0.831(0.082) 0.820(0.053) 0.7293 0.4387 0.3096 0.2381

Fig 4 .

 4 Fig 4. Estimated probabilities of joint high thresholds exceedances.

For

  any u ∈ [0, 1] d , C n (u) = H n (G ← η,1 (u 1 ), . . . , G ← η,d (u d )), so both terms can be expressed as continuous transformations of the empirical process H n . Therefore, the weak convergence of T n,1 + T n,2 follows from the continuous mapping theorem. A similar reasoning can be obtained if Condition 1(ii)(b) is replaced by Condition 1(ii)(a). In that case, we have the following decomposition√ n( M n ( η n , •) -M(η, •)) =: T n,1 + T n,2 + o p (1,t -E(W i,t )), T n,2 (t) :(v t 1 η , . . . , v t d η )1l {v>G 1/ηt j η

  t∈S d 1 η n η{1 -M n ( η n , t)}{1 -M(η, t)} × sup t∈S d η{1 -M(η, t)} M n ( η n , t) -η n {1 -M n ( η n , t)}M(η, t) =: T n,1 × T n,2 .imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018Since η n → η a.s., for a small ε > 0 and large n, we have almost surely that η (v 1+ε , . . . , v 1+ε )dv < ∞.Now, using the Lipschitz property of order k > 0 of C η , we haveT n,2 ≤ η{1 -M(η, t)} -η n {1 -M( η n , t)} ∞ + {1 -M n ( η n , t)}{1 -M(η, t)} ∞ | η n -η|

  u)u a (-log u) b du.

  n (u)u), u ∈ (0, 1), where for i.i.d. copies U 1 , . . . , U n of U = G η (max(Z 1 , . . . , Z d )), we denote i ≤ u), u ∈ (0, 1), and as before Q n := G ← n . Let Ġη (y) and Gη (y) be the first and second derivatives of G η (y) with respect to y > 0. The function defined in Theorem 3.2 is then equal toϕ(u) = Ġη (Q η (u)), u ∈ (0, 1). imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

  ) Q n (u) -Q n (u)| = o(1) a.s. (A.7)and by Bahadur-Kiefer theorem (see e.g.[START_REF] Einmahl | A short and elementary proof of the main Bahadur-Kiefer theorem[END_REF] we have sup u∈(0,1)| Q n (u) + H n (u)| = o(1) a.s.log u) ϕ(u) du + µ 2 1,1 + o(1) = µ 2 1,1 + o(1).

Table 1

 1 Estimated significance levels α. From left to right: the dimension of X, the true significance level, the approximate asymptotic (1 -α)-quantile, the empirical proportion of simulated samples under H 0 that rejected the null hypothesis and the empirical (1 -α)-quantile. Here ψ = 1.

	imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

Table 2

 2 Estimated significance levels α. From left to right: the dimension of X, the sample size and the empirical proportion of simulated samples under H 0 that rejected the null hypothesis for different values of ψ.

	ψ

Table 3

 3 Estimates (standard deviation) of η and MISE for the η-Pickands dependence function, based on a bivariate η-asymmetric logistic dependence model with η = 0.7. The first line corresponds to the GPWM method, whereas the second line is the ML method.

			η n				MISE( A ηn,n , A η )
			n				n
	ψ	25	50	100	200	25	50	100	200
	0.1 0.661(0.115) 0.678(0.084) 0.690(0.062) 0.695(0.043) 0.0111 0.0037 0.0013 0.0005
		0.800(0.201) 0.763(0.128) 0.741(0.085) 0.728(0.055) 0.0110 0.0036 0.0013 0.0005
	0.2 0.667(0.116) 0.679(0.084) 0.688(0.062) 0.692(0.044) 0.0480 0.0195 0.0088 0.0041
		0.807(0.204) 0.761(0.128) 0.740(0.088) 0.724(0.057) 0.0457 0.0187 0.0086 0.0040
	0.3 0.665(0.116) 0.680(0.087) 0.692(0.064) 0.696(0.046) 0.1176 0.0542 0.0262 0.0133
		0.811(0.211) 0.768(0.130) 0.745(0.087) 0.730(0.059) 0.1133 0.0527 0.0256 0.0131
	0.4 0.673(0.114) 0.687(0.088) 0.694(0.062) 0.697(0.045) 0.2177 0.1021 0.0523 0.0260
		0.810(0.204) 0.770(0.129) 0.744(0.084) 0.729(0.057) 0.2118 0.1000 0.0514 0.0257
	0.5 0.670(0.113) 0.684(0.085) 0.692(0.062) 0.695(0.044) 0.3602 0.1795 0.0952 0.0481
		0.805(0.201) 0.766(0.129) 0.742(0.085) 0.728(0.057) 0.3531 0.1765 0.0940 0.0476
	0.6 0.670(0.115) 0.685(0.085) 0.691(0.062) 0.696(0.046) 0.4566 0.2252 0.1192 0.0556
		0.822(0.206) 0.778(0.127) 0.751(0.089) 0.734(0.063) 0.4758 0.2275 0.1156 0.0585

Table 4

 4 Estimates (standard deviation) of η and MISE for the η-Pickands dependence function, based on componentwise maxima with approximate bivariate η-asymmetric logistic model with η = 0.7. The first line corresponds to the GPWM method, whereas the second line is the ML method.

  for every t ∈ S d , indexes s ∈ {1, . . . , p}, j ∈ {1, . . . , d} and for every |hh s | < δ. These results imply that sup

		t∈S d	max 1≤s≤p	sup v∈(0,1)	max 1≤ j≤d	|v η n,hs t j -v ηt j | → 0, n → ∞	(A.5)
	and
	sup t∈S d	max 1≤s≤p	sup |h-h s |<δ	max 1≤ j≤d	sup v∈(0,1)

  (v t 1 η , . . . , v t d η ) -Ċη, j (v t 1 η , . . . , v t d η )C n (1, . . . , 1, v t j η , 1, . . . , 1) + η A η ) 2 χ(H n ).

	where for all t ∈ S d
	T n,1 (t) := -	0	1	        C n d j=1	        dv
	and			T n,2 :=	(1	A η

1, A.2 with Proposition 3.1 in

[START_REF] Segers | Asymptotics of empirical copula processes under non-restrictive smoothness assumptions[END_REF]

, under Condition 1(ii)(b) we have that √ n( M n ( η n , •) -M(η, •)) = T n,1 (•) + T n,2 (•) + o p (

1), imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

Acknowledgements

We sincerely thank the editor, associate editor and the referees for their helpful comments and suggestions that led to substantial improvement of the paper. In addition, we are also gratefull to Carlo Baldassi for his help in parallelizing the code of the first simulation study, to Alexandra Ramos for her suggestions on the second simulation study and to Aad van der Vaart for having kindly replied to our technical questions.

The first author is supported by a research grant (VKR023480) from VILLUM FONDEN 1 imsart-generic ver.

y i, j are observations larger than a positive threshold s, m i = min(y i,1 , y i,2 ), n s is the number of m i exceeding the threshold s and m * is the (n s -1)th order statistic of the sample m 1 , . . . , m n s . When H η is finite, for j = 1, 2, r j (y) approaches a positive constant as y → m * /s, whereas in the infinite case, it decreases toward zero.

We illustrate the diagnostic tool with some examples. We draw samples of 500 × 366 values from six different models that satisfy equations (3.1) and (3.2). We consider three models with an η-asymmetric logistic (η-AL) survival function and parameters ψ = 0.1, ψ = 0.4 and ψ = 0.6, respectively, while u = 10, λ = 1e -0.1 -0.02, = 1, η = 0.7 are the same for all the three cases. In these examples H η is finite and both F min and F j , j = 1, 2, behave approximately as y -1/η for large values of y. Figure 3 displays in the first line the plots of r j , obtained using s = Q(0.07), where Q(0.07) is the ninety-third percentile of unit-Fréchet distribution. As expected, for large values of y, r j (y) stays away from zero and it approaches 1 when ψ = 0.1, 0.4 and a smaller constant when ψ = 0.6 (value close to ψ = 0.7 with which H η is infinite). We also consider a η-asymmetric logistic model with u = 8, λ = 1e -1/8 -0.04, = 1, η = 0.7 and ψ = 0.8, a bivariate standard Gaussian distribution with ρ = 0.5 and a bivariate inverted symmetric logistic (SL) model with ψ = 0.8 (see Section 2.3). In the latter two models the marginal distributions of the data are transformed into uni-Fréchet. In these three cases H η is infinite. Furthermore, for large y, F min behaves approximately as y -1/η with η = 0.7, η = (1 + ρ)/2 = 0.75 and η = 2 -ψ 0.57, respectively. While, F j , j = 1, 2 behaves approximately as y -1/k with k = ψ = 0.8, k = 1 + ρ = 1.5 and k = 2 1-ψ 1.15, respectively. For these three examples the diagnostic plots are displayed in the second line of Figure 3. As expected, r j (y) goes to zero for large values of y.

The procedure for inferring A η discussed in Section 3.2, when possible, provides useful means to extrapolate the probability of joint high thresholds exceedances as we describe next. For simplicity we focus on the bivariate case. By (3.1), (3.5) and the definition of the η-Pickands dependence function we have that the approximation imsart-generic ver. 2014/10/16 file: EtaPickandsHalv3.tex date: February 6, 2018

Appendix A: Proofs A.1. Some properties of A n Note that,

Therefore, A n (e j ) = 1 for all j = 1, . . . , d. The distribution function of the i.i.d. random variables X 1, j , . . . , X n, j , j = 1, . . . , d, being continuous, almost surely there are no ties and thus

For all η ∈ (0, 1] and t ∈ S d , set

By convention u 1/ηt = 0 when t = 0 and u ∈ [0, 1]. By Lemma A.1 in [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF] we have

.

The result (3.6) follows by solving the above equality for A η .

A.4. Proof of Theorem 3.1

We start with some notation. Let

, where C n is the empirical copula defined as