MSc Bjorn Stemkens
email: b.stemkens@umcutrecht.nl

Rob H Tijssen

Baudouin Denis De Senneville

MD Hanne D Heerkens

MD, PhD Marco Van Vulpen

PhD Jan Jw Lagendijk

PhD Cornelis At Van Den Berg

Dphil

PhD Baudouin D De Senneville

Optimizing 4-dimensional magnetic resonance imaging data sampling for respiratory motion analysis of pancreatic tumors

Keywords: 4D-MRI, radiotherapy, retrospective reconstruction, 4D-CT, non-rigid registration, pancreatic cancer, respiratory motion

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Radiotherapy treatments for pancreatic cancer is challenging due to motion, primarily induced by respiration. Conventionally, margins are used to ensure full dose coverage of the target at the cost of extra radiation to nearby organs at risk (OARs). Prior knowledge about patient specific motion, however, may lead to more conformal radiation plans [START_REF] Keall | The management of respiratory motion in radiation oncology report of AAPM Task Group 76a[END_REF].

Currently, 4D-CT is the method of choice for motion characterization in radiotherapy.

Unfortunately, CT has limited soft tissue contrast and 4D-CT often suffers from volume inconsistencies (e.g. duplicating, incomplete, or overlapping structures) caused by discrepancies between the breathing frequency and the pitch [START_REF] Abdelnour | Phase and amplitude binning for 4D-CT imaging[END_REF][START_REF] Yamamoto | Retrospective analysis of artifacts in four-dimensional CT images of 50 abdominal and thoracic radiotherapy patients[END_REF]. Magnetic Resonance Imaging (MRI), on the other hand, has superior soft tissue contrast, allowing better differentiation between tumor and surrounding tissue. Moreover, the flexibility of image acquisition in MRI could be used to overcome some of the image artifacts typical of 4D-CT.

Previous authors have reported MRI based motion characterization using 2D-cine MRI [START_REF] Feng | Characterization of pancreatic tumor motion using cine MRI: surrogates for tumor position should be used with caution[END_REF][START_REF] Bussels | Respiration-induced movement of the upper abdominal organs: a pitfall for the three-dimensional conformal radiation treatment of pancreatic cancer[END_REF][START_REF] Bjerre | Three-dimensional MRI-linac intra-fraction guidance using multiple orthogonal cine-MRI planes[END_REF][START_REF] Heerkens | MRI-based tumor motion characterization and gating schemes for radiation therapy of pancreatic cancer[END_REF].

However, these methods solely provide tumor motion characteristics, and do not provide information of the OARs. Consequently, retrospective binning methods have gained interest recently. To bin the acquired data into appropriate respiratory phases, a number of respiratory surrogates have been investigated, such as respiratory bellows [START_REF] Feng | Characterization of pancreatic tumor motion using cine MRI: surrogates for tumor position should be used with caution[END_REF][START_REF] Tryggestad | Respiration-based sorting of dynamic MRI to derive representative 4D-MRI for radiotherapy planning[END_REF], or image based body area [START_REF] Cai | Four-dimensional magnetic resonance imaging (4D-MRI) using image-based respiratory surrogate: A feasibility study[END_REF]. Another method is to use a 1D MRI navigator that records the diaphragm position [START_REF] Ehman | Adaptive technique for high-definition MR imaging of moving structures[END_REF].

To construct 4D-MR data, previous studies have reported retrospective stacking of 2D multi-slice images using different respiratory surrogates [START_REF] Tryggestad | Respiration-based sorting of dynamic MRI to derive representative 4D-MRI for radiotherapy planning[END_REF][START_REF] Cai | Four-dimensional magnetic resonance imaging (4D-MRI) using image-based respiratory surrogate: A feasibility study[END_REF][START_REF] Yang | Four-Dimensional Magnetic Resonance Imaging Using Axial Body Area as Respiratory Surrogate: Initial Patient Results[END_REF][START_REF] Siebenthal | 4D MR imaging of respiratory organ motion and its variability[END_REF]. These methods are closely related to 4D-CT and therefore suffer from similar volume inconsistencies. Volumetric 3D acquired MRI scans overcome these problems, since the entire volume is excited each time.

Moreover, signal-to-noise ratio (SNR) is inherently larger for 3D acquisitions. Buerger et al.

[13] used a custom 3D read-out and an extensive reconstruction process ( ) to create multiple 3D volumes, which is undesirable for clinical practice.

The goal of this study is to develop a 4D-MRI method to characterize respiratory induced motion of pancreatic tumors and its surrounding OARs using clinically available scan protocols. For this purpose, two motion surrogate signals that serve as input signal for the binning process are assessed. Moreover, two MRI sampling strategies and three binning methods are explored in terms of efficiency, imaging time and image quality. The final resulting technique is applied to three pancreatic cancer patients for motion characterization using a 3D optical flow (OF) algorithm, because of its short processing time and minimal number of control parameters [START_REF] Noe | Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy[END_REF].

Methods

We performed three different experiments in which the surrogate signal, MRI sampling strategy and binning method are independently optimized. The optimized protocol was used to scan three pancreatic cancer patients. Subsequently 3D motion of ten respiratory phases is calculated for both the tumor and the duodenum. All experiments are conducted on a 1.5T Philips Achieva MR scanner (Philips Healthcare, Best, The Netherlands). Processing is performed using MATLAB (The Mathworks, Natick, MA).

Surrogate signal

Two healthy subjects were scanned using a fast, coronal, 2D transient balanced steady-state gradient echo (balanced Turbo Field Echo, bTFE) sequence with radial read-out ( = , TE/TR = 1.3/2.6 , FOV = 294x294 , voxelsize = 1.5x1.5 , turbo-factor = 128, 500 dynamics, temporal frequency 2.6

). The imaging plane was angulated parallel to the spine to sample along the principal axis of pancreatic motion, thereby reducing through-plane motion. The image acquisition was interleaved with a 1D MRI navigator, which was positioned at the dome of the right hemi-diaphragm and respiratory bellows, provided by the MRI vendor, was placed on the upper abdomen.

To assess the agreement of the two surrogates with the pancreatic motion, we performed a general linear model (GLM) analysis, which is a least-squares fit to the data described by . Here is the pancreatic motion, obtained by a 2D OF algorithm [START_REF] De Senneville | MR-guided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF][START_REF] Roujol | Automatic nonrigid calibration of image registration for real time MR-guided HIFU ablations of mobile organs[END_REF][START_REF] Stemkens | Optical flow analysis on undersampled radial acquisitions for real-time tracking of the pancreas in MR guided radiotherapy[END_REF], the parameter estimates, the surrogate signal and represents the residual error. Spatial residual error ( ) maps were calculated to find the surrogate with the least residuals.

Furthermore, power spectrum analysis (PSA) was performed on the surrogate signals and the pancreatic motion to get insight in the agreement in frequency range of the surrogates with the pancreatic motion.

Binning Simulations

Another important consideration is how to reorder the continuously acquired MR data with respect to the respiratory cycle (i.e. retrospective binning). 4D-CT scanners employ two binning methods. Absolute amplitude and phase binning (see Fig. 1b). In the first method, the ten bins (i.e., respiratory phases) are determined by the minimum and maximum of the entire respiratory signal. The bins are thus fixed throughout the reordering process. The disadvantage of this method is that, when the breathing amplitude varies, certain phases will be sparsely sampled [START_REF] Abdelnour | Phase and amplitude binning for 4D-CT imaging[END_REF][START_REF] Ozhasoglu | Issues in respiratory motion compensation during external-beam radiotherapy[END_REF]. The phase binning method determines two consecutive exhale peaks and defines ten phases with equal time spacing. The disadvantage of this method, however, is that for an irregular breathing pattern, the amplitude variation within each bin can be much larger compared to absolute amplitude binning [START_REF] Lu | A comparison between amplitude sorting and phase-angle sorting using external respiratory measurement for 4D CT[END_REF]. To address these disadvantages, we introduce a third, hybrid, binning method, relative amplitude binning. This method defines the maximum exhale and inhale for each respiratory cycle and normalizes them before defining ten respiratory phases. We further distinguish between in-and exhale to be able to model hysteresis (i.e. varying motion paths for in-and exhale).

As retrospective reordering of continuously acquired data is a stochastic process, the binning method directly influences the dynamics of the reordering and thus the number of (continuously acquired) volumes required to fill each of the bins with sufficient data. To get insight into how quickly the bins are filled for each of the three binning methods, we performed various simulations. The navigator signal, which was acquired in the subsequent 3D in-vivo experiments, was used as input (8:35 minutes). The time to fill a percentage of all k-space lines required for the ten reconstructed volumes was calculated for an increasing number of acquired volumes. Moreover, to evaluate the intra-phase amplitude variation (which should be as small as possible) for the three binning methods, the mean amplitude and standard deviation for each phase was calculated based on the navigator amplitude.

In-vivo experiments

Seven healthy volunteers and three patients with an unresectable pancreas tumor provided written informed consent to be scanned according to institutional rules. A clinically available 3D bTFE was used with turbo direction along kz and a shot length equal to the amount of slices, acquiring all kz lines within a single shot (see Fig. 1a). Prior to each shot, a navigator was acquired.

We compared cartesian and radial in-plane sampling, using the parameters described in Table 1, in terms of image quality and robustness. 25 Dynamics were acquired as a concession between imaging time and filling. One subject was scanned with the cartesian protocol, whereas the other subjects and patients were scanned with the radial protocol. Two subjects, who were scanned with the radial protocol, had a smaller FOV (330x330x64 ), which was sufficient to cover the pancreatic motion path.

Reconstruction was performed offline on an 8-core CPU computer using ReconFrame reconstruction software (Gyrotools, Zurich, CH). Due to the stochastic nature of the filling of respiratory phases, it is likely that certain k-space segments are acquired multiple times during acquisition. When this occurred, complex averaging was performed to increase SNR. Image quality was assessed both visually and using the structural similarity (SSIM) index [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], which calculates the similarity between two images as a ratio between 0 (no similarity) and 1 (perfect similarity).

4D Motion Analysis

Data of three patients were reconstructed using relative amplitude binning. Non-rigid motion was calculated using a multi-threaded optimized 3D OF algorithm [START_REF] Noe | Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy[END_REF]. In short, this algorithm calculates the motion between two imaging frames based on intensity gradients with an additional constraint on motion smoothness to model elastic organ deformation. All respiratory phases were registered to one reference phase (the first phase) as can be seen in Fig. 1c. The tumor and OARs were delineated on the reconstructed reference volume, and the mean motion within each delineated structure was calculated for all respiratory phases.

To assess the quality of the motion estimation by the OF algorithm, all images were registered to the reference phase and visually inspected for residual motion. Moreover, the motion vectors of consecutive phases were added. Since this trajectory comprises a full respiratory cycle, the results should be close to zero. Furthermore, the OF results were compared to results obtained by Elastix [START_REF] Klein | Elastix: a toolbox for intensity-based medical image registration[END_REF]. Elastix is a widely used ITK-based non-rigid registration package, which uses B-spline registration, but requires much more computation time, and can therefore not be used in an online setting.

Results

Surrogate signal

Fig. 2a shows the residuals of the GLM using either the respiratory bellows or the MRI navigator as surrogate signal to model the motion within the abdomen for a representative subject along with an anatomical reference image. It is seen that the residuals within the pancreas are 3 times higher for the bellows compared to when the navigator is used. These maps further demonstrate the large motion heterogeneity over the FOV. Fig. 2b displays the power spectra of the pancreatic motion, and the spectra of the two surrogates. Aside from large respiratory peaks between 0.25 and 0.32 in all spectra, low frequencies ( 0.05 ) are present in the spectra of the pancreatic motion and the navigator, but absent in the bellows spectrum, which may indicate that the bellows are unable to pick up slow frequency position drifts. A diaphragm navigator is therefore used for all subsequent experiments.

Binning Simulations

Fig. 3a illustrates the amplitude variation within each phase for the three binning methods.

Absolute amplitude binning displays the smallest intra-phase variation, with a mean standard deviation over all subjects and phases of 2.0 . Phase binning reveals the largest amplitude variation, with a mean standard deviation of 3.9 , whereas relative amplitude binning has a mean standard deviation of 3.3 . Fig. 3b shows the ratio of the total k-space that is assigned to its appropriate respiratory phase as a function of the acquired number of dynamics. It reveals a rapid increase in k-space filling in the beginning, but an asymptotical approach towards due to the stochastic character of the respiratory filling. The average filling over all subjects after 25 dynamics is for absolute amplitude, for relative amplitude, and for phase binning. Fig. 3c, finally, shows the distribution of the assigned -space data across the ten phases. Missing data are denoted in white. The figure clearly demonstrates that the inhale phases (i.e. phase five/six) are only sparsely sampled with data when absolute amplitude binning is used. Considering the trade-off between low intra-phase amplitude variation and uniform k-space filing, relative amplitude binning will be used for the subsequent experiments.

In-vivo experiments

4D Motion Analysis

Discussion

In this paper, a novel 4D-MRI method was presented using only clinically available scan protocols. The total imaging and processing time was less than 15 minutes, which makes MRI based motion characterization of both the tumor and OARs clinically feasible.

Various 4D motion analyses were performed in three patients. Inter-organ and inter-patient variation was observed in the motion trajectories, which shows the clinical importance of patient specific motion models. Motion heterogeneity was apparent within the tumors, which may be caused by fixation of the tumor to blood vessels, due to ingrowth, which was confirmed by a radiologist.

The calculated patient specific 4D motion trajectories can be used for optimizing the radiation plan for both the GTV and the OARs in abdominal tumors. This can be used to transition from population based internal target volume (ITV) [START_REF] Underberg | Four-dimensional CT scans for treatment planning in stereotactic radiotherapy for stage I lung cancer[END_REF] to patient specific ITV. Moreover, the ability to perform motion characterization of OARs, such as the duodenum, allows for better definition of planning organs at risk volumes (PRVs) [START_REF] Mckenzie | Margins for geometric uncertainty around organs at risk in radiotherapy[END_REF], which will aid reducing toxicity.

Further, with the development of MR-Linac systems (e.g. [START_REF] Lagendijk | MRI/linac integration[END_REF]), the techniques from this study can be used to calculate the daily variation in motion, online, just prior to treatment and used for online tracking models and beam steering. Although the technique is only used for pancreatic tumors in this study, we believe it can be used for other abdominal tumors, which suffer from respiratory motion, e.g. liver and kidney tumors.

We have shown a 1D MRI navigator correlated better with pancreatic motion, than external respiratory bellows. Moreover, synchronization between the acquired image data and the navigator signal is more robust than for the bellows as mentioned previously [START_REF] Tryggestad | Respiration-based sorting of dynamic MRI to derive representative 4D-MRI for radiotherapy planning[END_REF]. It is worth pointing out, however, that the acquisition of the MR navigator does require some additional scanning time ( increase in scanning time).

The proposed relative amplitude binning showed the best trade-off between k-space filling efficiency and small intra-phase amplitude variation. Phase binning also performed well in the simulations, which is probably caused by the relatively stable respiratory signals used for these simulations. However, when breathing becomes more erratic, phase binning will likely result in larger intra-phase amplitude variation.

Finally, it was shown that radial in-plane sampling is superior to cartesian sampling. In radial sampling the center of -space is sampled for every read-out line, which makes the sequence more robust against image artifacts resulting from undersampling or intra-view motion [START_REF] Glover | Projection reconstruction techniques for reduction of motion effects in MRI[END_REF][START_REF] Scheffler | Reduced circular field-of-view imaging[END_REF].

In the reordered data sets this occurs when the bins are not completely filled or when the intra-phase variation is large (Fig. 4).

Limitations

Higher frequency induced motion (e.g. pulsatile motion) will not be visualized by retrospective respiratory binning, which is a limitation of any retrospective method. Furthermore, retrospective binning represents the motion of an average breathing cycle, whereas inter-and intra-cycle variation is not modeled, which may vary substantially (as was seen in the 2D data).

Moreover, the overall accuracy of the presented method is very difficult to determine. The retrospective binning, intra-phase amplitude variation and averaging will modulate k-space in a complex way. Future studies will therefore mainly focus on investigating the overall accuracy of the proposed technique.

Conclusion

A way to estimate 3D motion for ten respiratory phases based on retrospectively binning of k-space data using solely clinically scanning protocols has been presented. Optimized respiratory binning, based on an MR navigator, in combination with robust MR sampling resulted in ten 3D volumes. Using a fast 3D optical flow algorithm, we were able to calculate 4D motion for all abdominal structures within 13 minutes after starting acquisition.

Fig. 4a displays

  Fig. 4a displays the reconstructed images of the first phase (exhale) for cartesian and radial

Fig. 5

 5 Fig.5shows a coronal and sagittal slice with the vectorfield displayed on top in panels a-b

Figure 1

 1 Figure 1 Workflow of the in-vivo experiments. (a) Cartesian and radial trajectories with turbo

Figure 2

 2 Figure 2 (a) Residual errors of the GLM after regression with both surrogates in one

Figure 3

 3 Figure 3 (a) Amplitude variation within the respiratory phases for absolute amplitude, relative

Figure 4

 4 Figure 4 Resulting images for cartesian and radial in-plane sampling. (a) Reconstructions of

Figure 5

 5 Figure 5 Vectorfields displayed on top of (a) coronal and (b) sagittal slice of the first patient


Mikael Escobar-Bach

Yuri Goegebeur

Armelle Guillou

Local robust estimation of the Pickands dependence function

Keywords: MSC 2010 subject classifications: Primary 62G32, 62G05, 62G20; secondary 60F05, 60G70. Conditional Pickands dependence function, robustness, stochastic convergence

We consider the robust estimation of the Pickands dependence function in the random covariate framework. Our estimator is based on local estimation with the minimum density power divergence criterion. We provide the main asymptotic properties, in particular the convergence of the stochastic process, correctly normalized, towards a tight centered Gaussian process. The finite sample performance of our estimator is evaluated with a simulation study involving both uncontaminated and contaminated samples. The method is illustrated on a dataset of air pollution measurements.

Introduction

Modelling dependence among extremes is of primary importance in practical applications where extreme phenomena occur. To this aim, the copula function can be used as a margin-free description of the dependence structure. Indeed, according to the well-known result of [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF], the distribution function of a pair pY p1q , Y p2q q can be represented in terms of the two margins F 1 and F 2 of Y p1q and Y p2q respectively, and a copula function C as follows:

P ¡ Y p1q ¤ y 1 , Y p2q ¤ y 2
© C pF 1 py 1 q, F 2 py 2 qq .

This function C characterizes the dependence between Y p1q and Y p2q and is called an extreme value copula if and only if it admits a representation of the form Cpy 1 , y 2 q exp ¢ logpy 1 y 2 qA ¢ logpy 2 q logpy 1 y 2 q , where A: r0, 1s Ñ r1{2, 1s is the Pickands dependence function, which is convex and satisfies maxtt, 1¡tu ¤ Aptq ¤ 1, see [START_REF] Pickands | Multivariate extreme value distributions[END_REF]. Statistical inference on the bivariate function C is therefore equivalent to the statistical inference on the one-dimensional function A. Estimating this function A has been extensively studied in the literature. We can mention, among others, [START_REF] Capéraà | A nonparametric estimation procedure for bivariate extreme value copulas[END_REF], [START_REF] Fils-Villetard | Projection estimators of Pickands dependence functions[END_REF] or [START_REF] Bücher | New estimators of the Pickands dependence function and a test for extreme-value dependence[END_REF].

In this paper, we extend the above framework to the case where the pair pY p1q , Y p2q q is recorded along with a random covariate X R p . In that context, the copula function together with the marginal distribution functions depend on the covariate X. In the sequel, we denote by C x , F 1 p.|xq and F 2 p.|xq the conditional copula function and the continuous conditional distribution functions of Y p1q and Y p2q given X x. Our model can thus be written as

P ¡ F 1 pY p1q |xq ¤ y 1 , F 2 pY p2q |xq ¤ y 2 § § §X x © C x py 1 , y 2 q, (1.1)
where C x admits a representation of the form C x py 1 , y 2 q exp ¢ logpy 1 y 2 qA ¢ logpy 2 q logpy 1 y 2 q § § §x , ¦ This work was supported by a research grant (VKR023480) from VILLUM FONDEN and an international project for scientific cooperation (PICS-6416). Computation/simulation for the work described in this paper was supported by the DeIC National HPC Centre, SDU.

with Ap.|.q : r0, 1s¢R p Ñ r1{2, 1s is the conditional Pickands dependence function which is again a convex function satisfying maxtt, 1¡tu ¤ Apt|xq ¤ 1 for all x R p . From a practical point of view, the considered family of extreme value distributions has sufficiently large potential for data analysis. Firstly, the family of extreme value distributions is very rich, and includes among others the logistic, the asymmetric logistic, the negative logistic, the Hüsler-Reiss, the t extreme value, and Dirichlet model. Secondly, multivariate extreme value distributions arise naturally as the limiting distributions of properly normalised componentwise maxima, making them a useful approximation to the true, but typically unknown, distribution of these component-wise maxima in practice. We refer to [START_REF] Kotz | Extreme value distributions -Theory and applications[END_REF] and [START_REF] Gudendorf | Extreme-value copulas, in: Copula theory and its applications[END_REF] for further motivation and discussion of this class of distributions and additional examples. As a possible application, we consider modelling extremal dependence between air pollutants, like groundlevel ozone and particulate matter, conditional on location and time; see Section 5 for more details. Moreover, in addition to the covariate context, we consider the case of contamination and we propose a robust estimator of the conditional Pickands dependence function Ap.|xq. To reach this goal, we use the density power divergence method introduced by [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF]. In particular, the density power divergence between two density functions g and h is defined as follows ∆ α pg, hq : Here the density function h is assumed to depend on a parameter vector θ, and if Z 1 , ..., Z n is a sample of independent and identically distributed random variables according to the density function g, then the minimum density power divergence estimator (MDPDE) of θ is the point p θ minimizing the empirical version (up to a constant independent of θ) p ∆ α pθq : 6 9 9 9 8 9 9 9 7

» R h 1 α pyqdy ¡ ¢ 1 1 α 1 n n i1 h α pZ i q, α ¡ 0, ¡ 1 n n i1
log hpZ i q α 0.

We can observe that for α 0 one recovers the log-likelihood function, up to the sign. A large value of α allows us to increase the robustness of the estimator, whereas a smaller value implies more efficiency. This parameter α can thus be selected in order to ensure a trade-off between these two antagonist concepts. The nonparametric estimation of extremal dependence in presence of random covariates is still in its infancy, despite the huge potential of such methods for practical data analysis. [START_REF] Gardes | Nonparametric estimation of the conditional tail copula[END_REF] introduce an estimator for the tail copula based on a random sample from a distribution in the maxdomain of attraction of an extreme value distribution, and provide a finite dimensional convergence result for their estimator, when properly normalised. [START_REF] Portier | On the weak convergence of the empirical conditional copula under a simplifying assumption[END_REF] considered model (1.1) but under the simplifying assumption that the dependence between Y p1q and Y p2q does not depend on the value taken by the covariate, i.e. C x C (see also [START_REF] Gijbels | Estimation of a copula when a covariate affects only marginal distributions[END_REF]. In the present paper we introduce a nonparametric and robust estimator for Ap.|xq which is obtained by an adjustment of the above introduced density power divergence estimation criterion to the situation of local estimation, and we study the asymptotic properties of the obtained estimator in terms of stochastic process convergence. To the best of our knowledge, nonparametric and robust estimation of the conditional Pickands dependence function has not been considered in the literature.

The remainder of the paper is organized as follows. In Section 2, we simplify the situation to the case where the two marginal distributions are known, we propose a robust estimator for Ap.|xq and prove its convergence in terms of a stochastic process. Then, in Section 3, we extend this result to the case of unknown margins. The efficiency and robustness of the estimator are examined with a simulation study, described in Section 4. Finally, in Section 5 we illustrate the practical applicability of the method for modelling extremal dependence between air pollution measurements. Additional simulation results are available in the online supplementary material. All the proofs are postponed to the Appendix.

Case of known margins

We denote by f the density function of the covariate X and by x 0 a reference position such that x 0 IntpS X q, the interior of the support S X of f . In this section, we restrict our interest to the case where the marginals F 1 p.|xq and F 2 p.|xq are known, and we denote by A 0 p.|xq the true conditional Pickands dependence function associated to the pair pY p1q , Y p2q q.

Construction of the estimator

For convenience we reformulate the model (1.1) into standard exponential margins. After applying the transformations r Y pjq ¡ log F j pY pjq |xq, j 1, 2, we obtain the following bivariate survival function

Gpy 1 , y 2 |xq : P ¡ r Y p1q ¡ y 1 , r Y p2q ¡ y 2 § § §X x © exp ¢ ¡py 1 y 2 qA 0 ¢ y 2 y 1 y 2 § § §x
, for all y 1 , y 2 ¡ 0. Let t r0, 1s. Considering the univariate random variable

Z t : min £ r Y p1q 1 ¡ t , r Y p2q t ,
it is clear that PpZ t ¡ z|X xq e ¡zA0pt|xq , dz ¡ 0 and x R p . Consequently, the conditional distribution of Z t given X x is an exponential distribution with parameter A 0 pt|xq.

Let pZ t,i , X i q, i 1, . . . , n, be independent copies of the random pair pZ t , Xq. In the present paper, we will develop a nonparametric robust estimator for A 0 pt|x 0 q by fitting this exponential distribution function locally to the variables Z t,i , i 1, ..., n, by means of the MDPD criterion, adjusted to locally weighted estimation, i.e. we minimize for α ¡ 0 p ∆ α,x0,t paq : 1

n n i1 K h px 0 ¡ X i q 4» V 0 ae ¡az ¨1 α dz ¡ ¢ 1 1 α ae ¡aZt,i ¨αB a α n n i1 K h px 0 ¡ X i q 4 1 1 α ¡ ¢ 1 1 α e ¡αaZt,i B .
(2.1)

Here K h p.q : Kp.{hq{h p where K is a joint density on R p and h h n is a positive non-random sequence satisfying h n Ñ 0 as n Ñ V. The MDPDE p A α,n pt|x 0 q for A 0 pt|x 0 q satisfies the estimating equation

p ∆ p1q α,x0,t p p A α,n pt|x 0 qq 0, (2.2)
where p ∆ pjq α,x0,t p.q denotes the derivative of order j of p ∆ α,x0,t p.q. The minimization of p ∆ α,x0,t is here performed without constraints, which means that p A α,n p.|x 0 q does not automatically satisfy the conditions of the Pickands dependence function. In fact, this is the case for several of the estimators proposed in the literature, see e.g. [START_REF] Pickands | Multivariate extreme value distributions[END_REF], [START_REF] Deheuvels | On the limiting behavior of the Pickands estimator for bivariate extreme-value distributions[END_REF] or [START_REF] Capéraà | A nonparametric estimation procedure for bivariate extreme value copulas[END_REF]. To overcome this, one could follow the idea of [START_REF] Fils-Villetard | Projection estimators of Pickands dependence functions[END_REF], and project the obtained estimator onto the space of Pickands dependence functions.

Our aim in this paper is to show the weak convergence of the stochastic process

3 c nh p ¡ p A α,n pt|x 0 q ¡ A 0 pt|x 0 q © , t r0, 1s A , (2.3)
in the space of all continuous functions on r0, 1s, denoted as Cpr0, 1sq, when n Ñ V.

Our starting point is the estimating equation (2.2). By applying a Taylor series expansion around the true value A 0 pt|x 0 q, we get

0 p ∆ p1q α,x 0 ,t pA0pt|x0qq ¡ p Aα,npt|x0q ¡ A0pt|x0q © p ∆ p2q α,x 0 ,t pA0pt|x0qq 1 2 ¡ p Aα,npt|x0q ¡ A0pt|x0q © 2 p ∆ p3q α,x 0 ,t p r

Apt|x0qq

where r

Apt|x 0 q is a random value between A 0 pt|x 0 q and p A α,n pt|x 0 q. A straightforward rearrangement of the above display gives

c nh p ¡ p Aα,npt|x0q ¡ A0pt|x0q © ¡ c nh p p ∆ p1q α,x 0 ,t pA0pt|x0qq p ∆ p2q α,x 0 ,t pA0pt|x0qq 1 2 p ∆ p3q α,x 0 ,t p r Apt|x0qq ¡ p Aα,npt|x0q ¡ A0pt|x0q © . (2.4)
Consequently, in order to obtain the convergence of the stochastic process (2.3), we need to study the properties of the derivatives p ∆ pjq α,x0,t , j 1, 2, 3. According to Appendix A.5, these can be expressed as a linear combination of a key statistic T n , defined as

T n pK, a, t, λ, β, γ|x 0 q : a γ n n i1 K h px 0 ¡ X i qZ β t,i e ¡λaZt,i , (2.5)
for a r1{2, 1s, t r0, 1s, λ, β ¥ 0 and γ R.

Asymptotic properties of T n

Due to the regression context, we need some Hölder-type conditions on the density function f and on the conditional Pickands dependence function A 0 . Let }.} be some norm on R p , and denote by B x prq the closed ball with respect to }.} centered at x and radius r ¡ 0. Assumption pDq. There exist M f ¡ 0 and η f ¡ 0 such that |fpxq ¡ f pzq| ¤ M f }x ¡ z} η f , for all px, zq S X ¢ S X . Assumption pA 0 q. There exist M A0 ¡ 0 and η A0 ¡ 0 such that |A 0 pt|xq ¡ A 0 pt|zq| ¤ M A0 }x ¡ z} η A 0 , for all px, zq B x0 prq ¢ B x0 prq, r ¡ 0 and t r0, 1s.

Also a usual condition is assumed on the kernel K.

Assumption pK 1 q. K is a bounded density function on R p with support S K included in the unit ball of R p with respect to the norm }.}.

As a preliminary result, in Lemma 2.1 we prove the convergence in probability of the key statistic T n .

Lemma 2.1. Assume that for all t r0, 1s, x Ñ A 0 pt|xq and the density function f are both continuous at x 0 IntpS X q non-empty. Under Assumption pK 1 q, if h Ñ 0 and nh p Ñ V, then for a r1{2, 1s, λ, β ¥ 0, γ R and x 0 such that f px 0 q ¡ 0, we have

T n pK, a, t, λ, β, γ|x 0 q P ÝÑa γ Γpβ 1q A 0 pt|x 0 q pλa A 0 pt|x 0 qq β 1 f px 0 q, as n Ñ V, where Γ is the gamma function defined as Γprq : ³ V 0 t r¡1 e ¡t dt, dr ¡ 0. Now, our interest is in the rate of convergence in Lemma 2.1 when a is replaced by A 0 pt|x 0 q. More precisely, we want to show the weak convergence of the stochastic process

4 c nh p ¢ T n pK, A 0 pt|x 0 q, t, λ, β, γ|x 0 q ¡ Γpβ 1q rA 0 pt|x 0 qs γ¡β pλ 1q β 1 f px 0 q , t r0, 1s B .
To establish such a result, we use empirical processes arguments based on the theory of Vapnik-Červonenkis classes (VC -classes) of functions as formulated in van der [START_REF] Van Der Vaart | Weak convergence and empirical processes, with applications to statistics[END_REF]. This allows us to show the following theorem.

Theorem 2.1. Let γ R and pλ, βq p0, Vq ¢ R or pλ, βq p0, 0q. Under the assumptions of Lemma 2.1 and if pDq and pA 0 q hold with c nh p h minpη f ,η A 0 q Ñ 0, then the process 4 c nh p ¢ T n pK, A 0 pt|x 0 q, t, λ, β, γ|x 0 q ¡ Γpβ 1q rA 0 pt|x 0 qs γ¡β pλ 1q β 1 f px 0 q , t r0, 1s B weakly converges in Cpr0, 1sq towards a tight centered Gaussian process tB t , t r0, 1su with covariance structure given by CovpB t , B s q rA 0 pt|x 0 qA

0 ps|x 0 qs γ }K} 2 2 f px 0 q ¢ 5 » R 2 gpu, vqG t,s pu, v|x 0 qdudv 1 ¡ λ 1 λ δ 0 pβq C ,
for all ps, tq r0, 1s 2 , where δ 0 is the Dirac measure on 0, and gpu, vq :

u β¡1 pβ ¡ λA 0 pt|x 0 quqe ¡λA0pt|x0qu v β¡1 pβ ¡ λA 0 ps|x 0 qvqe ¡λA0ps|x0qv , G t,s pu, v|x 0 q : G ¡ maxpp1 ¡ tqu, p1 ¡ sqvq, maxptu, svq § § §x0 © .
We now derive the limiting distribution of a vector of statistics of the form (2.5), when properly normalized. Let T n be a pm ¢ 1q vector defined as

Tn : pTnpK, A0pt1|x0q, t1, λ1, β1, γ1|x0q, ..., TnpK, A0ptm|x0q, tm, λm, βm, γm|x0qq T , for some positive integer m and let Σ be a pm¢mq covariance matrix with elements pσ j,k q 1¤j,k¤m defined as

σ j,k : rA 0 pt j |x 0 qs γj rA 0 pt k |x 0 qs γ k }K} 2 2 f px 0 q ¢ 5 » R 2 g j,k pu, vqG tj ,t k pu, v|x 0 qdudv δ 0 pβ j q Γpβ k 1q rλ k 1s β k 1 rA 0 pt k |x 0 qs β k δ 0 pβ k q Γpβ j 1q
rλ j 1s βj 1 rA 0 pt j |x 0 qs βj ¡ δ 0 pβ j qδ 0 pβ k q B (2.6) where g j,k pu, vq :

u βj ¡1 rβ j ¡ λ j A 0 pt j |x 0 qus e ¡λjA0ptj|x0qu v β k ¡1 rβ k ¡ λ k A 0 pt k |x 0 qvs e ¡λ k A0pt k |x0qv .
The aim of next theorem is to provide the finite dimensional convergence result which will, together with the tightness, allow us to establish the joint convergence of several processes related to the statistic T n .

Theorem 2.2. Under the assumptions of Lemma 2.1, we have

c nh p pT n ¡ ErT n sq N m p0, Σq ,
where N m denotes a m¡dimensional normal distribution.

We have now all the needed ingredients for proving the asymptotic properties of the MDPDE p A α,n pt|x 0 q.

Asymptotic properties of p

A α,n pt|x 0 q

The first result states the existence and uniform consistency of a sequence of solutions to the estimating equation (2.2).

Theorem 2.3. Let α ¡ 0. Under the assumptions of Theorem 2.1, with probability tending to 1, there exists a sequence

¡ p

A α,n pt|x 0 q © nN of solutions for the estimating equation (2.2) such that

sup tr0,1s § § § p A α,n pt|x 0 q ¡ A 0 pt|x 0 q § § § o P p1q.
Now, we come back to our final goal which is the weak convergence of the stochastic process (2.3).

Theorem 2.4. Let

¡ p

A α,n pt|x 0 q © nN be the consistent sequence defined in Theorem 2.3. Under the assumptions of Theorem 2.1, the process

3 c nh p ¡ p A α,n pt|x 0 q ¡ A 0 pt|x 0 q © , t r0, 1s
A weakly converges in Cpr0, 1sq towards a tight centered Gaussian process tN t , t r0, 1su with covariance structure given by

Cov pN t , N s q }K} 2 2 A 0 pt|x 0 qA 0 ps|x 0 q f px 0 q p1 αq 2 p1 α 2 q 2 v T α Σpt, sqv α ,
where

v α : ¤ ¦ ¥ α 1 α ¡p1 αq 1 α
and Σpt, sq :

¤ ¥ p1 αq 2 1 α 1 1 α Σ 2,2 pt, sq Σ 2,3 pt, sq 1 Σ 2,3 ps, tq Σ 3,3 pt, sq with Σ 2,2 pt, sq : p1 ¡ αqp1 αq α 2 p1 αq 2 A 0 pt|x 0 qA 0 ps|x 0 q » R 2 e ¡αrA 0 pt|x 0 qu A 0 ps|x 0 qvs Gt,spu, v|x 0 qdudv Σ 2,3 pt, sq : 1 ¡ αp1 αq 2 A 0 pt|x 0 qA 0 ps|x 0 q » R 2 p1 ¡ αA 0 ps|x 0 qvqe ¡αrA 0 pt|x 0 qu A 0 ps|x 0 qvs Gt,spu, v|x 0 qdudv Σ 3,3 pt, sq : p1 αq 2 A 0 pt|x 0 qA 0 ps|x 0 q » R 2 p1 ¡ αA 0 pt|x 0 quqp1 ¡ αA 0 ps|x 0 qvqe ¡αrA 0 pt|x 0 qu A 0 ps|x 0 qvs Gt,spu, v|x 0 qdudv.
In particular, for all t r0, 1s, we have

c nh p ¡ p Aα,npt|x 0 q ¡ A 0 pt|x 0 q © N 1 ¢ 0, }K} 2
2 rA 0 pt|x 0 qs 2 f px 0 q p1 αq 2 p1 4α 9α 2 14α 3 13α 4 8α 5 4α 6 q p1 2αq 3 p1 α 2 q 2 , as n Ñ V.

The asymptotic standard deviation is shown as a function of α in Figure S1 of the supplementary material. As is clear from this plot, the asymptotic standard deviation is increasing in α. Note that our results could also be obtained under different assumptions by using the local empirical process results of [START_REF] Stute | Conditional empirical processes[END_REF] and [START_REF] Einmahl | Gaussian approximation of local empirical processes indexed by functions[END_REF], combined with the functional delta method.

Case of unknown margins

In this section, we consider the general framework where both F 1 p.|xq and F 2 p.|xq are unknown conditional distribution functions. We want to mimic what has been done in the previous section and transform to standard exponential margins. To this aim, we consider the triplets

¡ ¡ log ¡ F n,1 pY p1q i |X i q © , ¡ log ¡ F n,2 pY p2q i |X i q © , X i ©
, i 1, ..., n, for suitable estimators F n,j of F j , j 1, 2, and we compute the univariate random variables q Z n,t,i : min

¤ ¥ ¡ log ¡ F n,1 pY p1q i |X i q © 1 ¡ t , ¡ log ¡ F n,2 pY p2q i |X i q © t , i 1, ..., n.
Then, similarly as in Section 2, the statistic

q T n pK, a, t, λ, β, γ|x 0 q : a γ n n i1 K h px 0 ¡ X i q q Z β n,t,i e ¡λa q Zn,t,i , (3.1)
is the cornerstone for the MDPDE, denoted q A α,n pt|x 0 q, which satisfies the estimating equation

q ∆ p1q α,x0,t p q A α,n pt|x 0 qq 0, (3.2)
where q ∆ p1q α,x0,t p.q is the first derivative of q ∆ α,x0,t p.q and q ∆ α,x0,t paq :

a α n n i1 K h px 0 ¡ X i q 4 1 1 α ¡ ¢ 1 1 α e ¡αa q Zn,t,i B .
The final goal is still the same, that is the weak convergence of the stochastic process

3 c nh p ¡ q A α,n pt|x 0 q ¡ A 0 pt|x 0 q © , t r0, 1s A . (3.3)
Again this result relies essentially on the asymptotic properties of the statistic q T n , and so the idea will be to decompose

c nh p ¡ q T n ¡ Er q T n s © pK, a, t, λ, β, γ|x 0 q, into the two terms 3 c nh p pTn ¡ ErTnsq pK, a, t, λ, β, γ|x0q A 3 c nh p ¡ r q Tn ¡ Tns ¡ Er q Tn ¡ Tns © pK, a, t, λ, β, γ|x0q A . (3.4)
The first term can be dealt with using the results of Section 2.2, whereas we have to show that the second term is uniformly negligible. To achieve this objective, let us introduce the following empirical kernel estimator of the unknown conditional distribution functions

F n,j py|xq : °n i1 K c px ¡ X i q1l tY pjq i ¤yu °n i1 K c px ¡ X i q , j 1, 2,
where c : c n is a positive non-random sequence satisfying c n Ñ 0 as n Ñ V. Here we kept the same kernel K as in the previous section, but of course any other kernel function can be used.

Before stating our main results, we need to impose again some assumptions, in particular a Hölder-type condition on each marginal conditional distribution function F j similar to the one imposed on the density function of the covariate.

Assumption pFq. There exist M Fj ¡ 0 and η Fj ¡ 0 such that |F j py|xq ¡ F j py|zq| ¤ M Fj }x ¡ z} η F j , for all y R and all px, zq S X ¢ S X , and j 1, 2. Concerning the kernel K a stronger assumption than pK 1 q is needed. Assumption pK 2 q. K satisfies Assumption pK 1 q, there exists δ, m ¡ 0 such that B 0 pδq S K and Kpuq ¥ m for all u B 0 pδq, and K belongs to the linear span (the set of finite linear combinations) of functions k ¥ 0 satisfying the following property: the subgraph of k, tps, uq : kpsq ¥ uu, can be represented as a finite number of Boolean operations among sets of the form tps, uq : qps, uq ¥ ϕpuqu, where q is a polynomial on R p ¢ R and ϕ is an arbitrary real function.

The latter assumption has already been used in [START_REF] Giné | Rates of strong uniform consistency for multivariate kernel density estimators[END_REF] or [START_REF] Giné | Weighted uniform consistency of kernel density estimators[END_REF]. In particular, it allows us to measure the discrepancy between the conditional distribution function F j and its empirical kernel version F n,j .

Lemma 3.1. Assume that there exists b ¡ 0 such that f pxq ¥ b, dx S X R p , f is bounded, and pK 2 q and pFq hold. Consider a sequence c tending to 0 as n Ñ V such that for some q ¡ 1 | log c| q nc p ÝÑ 0.

Also assume that there exists an ε ¡ 0 such that for n sufficiently large

inf xS X λ ptu B 0 p1q : x ¡ cu S X uq ¡ ε, (3.5)
where λ denotes the Lebesgue measure. Then, for any 0 η minpη F1 , η F2 q, we have

sup py,xqR¢S X |F n,j py|xq ¡ F j py|xq| o P £ max £ | log c| q nc p , c η , for j 1, 2.
Note that the assumption f pxq ¥ b, dx S X , for some b ¡ 0, is similar to the one already used in [START_REF] Gijbels | Estimation of a copula when a covariate affects only marginal distributions[END_REF] and [START_REF] Portier | On the weak convergence of the empirical conditional copula under a simplifying assumption[END_REF].

We are now able to study the second term in (3.4).

Theorem 3.1. Assume that there exists b ¡ 0 such that f pxq ¥ b, dx S X R p , f is bounded, and pK 2 q, pDq and pFq hold together with condition (3.5). Consider two sequences h and c tending to 0, such that for nh p Ñ V and for some q ¡ 1 and any 0 η minpη F1 , η F2 q c nh p r n :

c nh p max £ | log c| q nc p , c η ÝÑ 0,
as n Ñ V. Then, for all γ R and pλ, βq p0, Vq ¢ R or pλ, βq p0, 0q, we have

sup tr0,1s,ar1{2,1s c nh p § § § q T n ¡ T n ¡ E q T n ¡ T n % § § § pK, a, t, λ, β, γ|x 0 q o P p1q.
Finally, the decomposition (3.4) combined with Theorem 3.1 and the results from Section 2.2, yields the desired theoretical result of this paper.

Theorem 3.2. Let α ¡ 0. Under the assumptions of Theorem 3.1 and pA 0 q, with probability tending to 1, there exists a sequence

¡ q

A α,n pt|x 0 q © nN of solutions for the estimating equation (3.2) such that

sup tr0,1s § § § q A α,n pt|x 0 q ¡ A 0 pt|x 0 q § § § o P p1q.
Moreover, for this consistent sequence, if

c nh p h minpη f ,η A 0 q Ñ 0, the process 3 c nh p ¡ q A α,n pt|x 0 q ¡ A 0 pt|x 0 q © , t r0, 1s
A weakly converges in Cpr0, 1sq towards the tight centered Gaussian process tN t , t r0, 1su defined in Theorem 2.4.

Simulation study

Our aim in this section is to illustrate the efficiency of the proposed robust estimator for the conditional Pickands dependence function with a simulation study. We assume that the conditional distribution function of pY p1q , Y p2q q given X x is a mixture model of the form F ε py 1 , y 2 |xq p1 ¡ εqF py 1 , y 2 |xq εF c py 1 , y 2 |xq, where ε r0, 1s represents the fraction of contamination in the dataset. The main distribution F is the logistic distribution given by

F py 1 , y 2 |xq : exp 3 ¡ ¡ y ¡1{x 1 y ¡1{x 2 © x A , for y 1 , y 2 ¥ 0 and A 0 pt|xq ¡ t 1{x p1 ¡ tq 1{x © x ,
where the covariate X is a uniformly distributed random variable on r0, 1s. For this model, complete dependence is obtained in the limit as x Ó 0, whereas independence can be reached for x 1. Note also that the conditional marginal distributions of Y pjq given X, j 1, 2, under this logistic model are unit Fréchet distributions. Moreover, we can check that this model satisfies conditions pDq, pA 0 q and pFq.

Two completely different types of distributions F c will be considered throughout the paper and additional examples will be given in the online supplementary material.

First type of contamination: Given X x, the distribution function F c is

F c py 1 , y 2 |xq 1 2 3 e ¡y ¡1 1 e ¡y ¡1 2 A 1l ty1¥0,y2¥0u .
The mixture based on this distribution F c is illustrated in Figure S2 of the supplementary material. For this mixture, the contaminated points are on the axes.

Second type of contamination: The distribution function F c has completely dependent unit exponential margins. Figure S3 in the supplementary material shows an example of a simulated dataset from this model. This time, the contaminated points are on the diagonal.

To compute the estimator q A α,n , two sequences h and c have to be chosen. Both are determined by a cross validation criterion. Because of the very high computational burden of the cross validations for sample sizes n ¥ 1000, a random selection of size n r : n 1000 from the original observations is obtained, denoted tpY p1q i,r , Y p2q i,r , X i,r qu i1,...,nr , and the cross validations are implemented on these random subsamples. Concerning c, we can use the following cross validation criterion, already used in an extreme value context by [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF]:

c j : arg min r cj C nr i1 nr ķ1 1l ! Y pjq i,r ¤Y pjq k,r ) ¡ r F nr,¡i,j pY pjq k,r |X i,r q & 2 , j 1, 2,
where C is a grid of values of r c j and r F nr,¡i,j py|xq :

°nr k1,k$i K r cj px ¡ X k,r q1l tY pjq k,r ¤yu °nr k1,k$i K r cj px ¡ X k,r q .
Also the bandwidth parameter h is selected using a cross validation criterion. In particular, h : arg min

r hH 1 n r M nr i1 M j1 q A α,n,p¡iq pt j |X i,r q α ¢ 1 1 α ¡ ¢ 1 1 α e ¡α q A α,n,p¡iq ptj|Xi,rq q Zn,t j ,i,r for α ¡ 0 h : arg min r hH 1 n r M nr i1 M j1 ¡ log ¡ q
A 0,n,p¡iq pt j |X i,r q e ¡ q A 0,n,p¡iq ptj|Xi,rq q Zn,t j ,i,r © for α 0 where q A α,n,p¡iq pt|xq denotes the estimator of A 0 pt|xq obtained on all but observation i, q Z n,tj ,i,r is as q Z n,tj ,i but now calculated for pY p1q

i,r , Y p2q i,r , X i,r q, and q A 0,n pt|xq :

°n i1 K r h px ¡ X i q °n i1 K r h px ¡ X i q q Z t,i .
This criterion can be seen as a generalisation of a commonly used cross validation from the context of local likelihood estimation (see e.g. [START_REF] Abegaz | Semiparametric estimation of conditional copulas[END_REF] to the context of local MDPD estimation.

After extensive simulation studies, we have chosen the grids C t0.06, 0.12, . . . , 0.3u and H t0.02, 0.03, . . . , 0.06u. These choices provide a reasonable trade off between stability of the estimates and accuracy of approximation by asymptotic results.

Concerning the kernel, each time we use the bi-quadratic function Kpxq : 15 16 p1 ¡ x 2 q 2 1l r¡1,1s pxq.

As an indicator of efficiency, we compute over a grid the L 2 -error in the estimation of the Pickands dependence function Ap.|xq as a function of x, i.e.

MISEpε, α|xq :

1

N M N i1 M m1 q A piq α,ε,n pt m |xq ¡ A 0 pt m |xq % 2 .
Here q

A piq α,ε,n pt m |xq is our estimator of A 0 pt m |xq obtained with the i¡th sample when the contamination is ε. We set t m m{50, m 1, . . . , 49. Our simulations are based on datasets of sizes n 1000 and n 5000, and the procedure is repeated N 200 times.

Figure 1 represents the MISEpε, α|xq as a function of ε t0, 0.025, 0.05, ..., 0.2u for the two types of contamination (rows 1 and 2, and 3 and 4, respectively). From the left to the right, three positions have been considered: x 0.1, 0.3 and 0.5 for the first type of contamination and x 0.5, 0.7 and 0.9 for the second type. Also, four different values of α have been reported: 0 (black), 0.1 (blue), 0.5 (green) and 1 (red), and two sample sizes: n 1000 and n 5000. Based on these simulations, we can draw the following conclusions:

• As expected, the MISE curves show less variability for n 5000 compared to n 1000;

• If the percentage of contamination ε is very small (in the range 0 -0.025), the MISE indicators are typically very similar, whatever the value of α. This result is a nice feature of our method, because if there is almost no contamination then in principle one does not need a robust procedure, but as is clear from these figures, the MDPDE performs similarly to the maximum likelihood method (corresponding to α 0), which is efficient (but not robust);

• If we increase the percentage of contamination ε, then it is crucial to increase the value of α to 0.5 or 1 in order to have good results. Indeed, for increasing ε a small value of α implies a drastic increase of the MISE;

• The MISE values are almost constant for α 0.5 and 1, whatever the percentage of contamination.

This illustrates again the robustness of our method, since it means that the methodology can handle a quite large percentage of contamination without deterioration of the results;

• For the first type of contamination, the gain in MISE by taking α 0.5 or 1 over α 0 or 0.1 is more important for small x than for large x. In this case the contamination is on the axes (independently), and this is less disturbing for x close to 1, which corresponds also to independence, than for x close to 0, which corresponds to complete dependence. For the second type of contamination one can observe the opposite effect. The gain of taking α 0.5 or 1 over α 0 or 0.1 is more important for

x close to 1 than for x close to zero. Indeed, the perfectly dependent contamination is less disturbing for small x than for large x; • Figure 1 gives us also some indications about the breakdown point of our estimator, which is a common concept in the robust framework. Indeed the breakdown point can be interpreted as the smallest ε where the MISE indicator starts to increase. For α small, in the range 0 -0.1, the breakdown point is very small, say ε around 0.025, while for α 0.5 and 1, one can go to ε 0.15 or a larger value, depending on the type of contamination, which illustrates again the nice robustness property of our method.

Next to the above mentioned MISE indicators, we also used our simulated data to compute the empirical coverage probabilities of 90% confidence intervals based on the limiting distribution given in MISEpε, α|xq as a function of ε t0, 0.025, 0.05, ..., 0.2u with α 0 (black), α 0.1 (blue), α 0.5 (green) and α 1 (red). First type of contamination: rows 1 pn 1000q and 2 pn 5000q, x 0.1, 0.3, 0.5 from the left to the right. Second type of contamination: rows 3 pn 1000q and 4 pn 5000q, x 0.5, 0.7, 0.9 from the left to the right. Theorem 2.4. These are given in Tables 1 and2 for the first and second type of contamination, respectively. From these tables we can see that

• When ε 0 then the empirical coverage probabilities are generally larger than 0.90, meaning that the confidence interval based on the limiting distribution is conservative; • For increasing ε, the coverage probabilities generally decrease when α is small (0 or 0.1), while for larger α, especially for α 1, the coverage probabilities do not seem to be as much affected by the contamination.

Since the main objective of this paper is to estimate the conditional Pickands dependence function Ap.|xq, we also provide in the online supplementary material the boxplots of our estimator q A α,n p.|xq based on 200 replications for the two examples of contamination introduced above and two additional examples. These figures emphasize again the robustness properties of our estimator.

Application to air pollution data

In this section we illustrate the practical applicability of our method on a dataset of air pollution measurements. Extreme temperature and high levels of pollutants like ground-level ozone and particulate matter pose a major threat to human health. We consider the data collected by the United States Environmental Protection Agency (EPA), publicly available at https:{{aqsdr1.epa.gov{aqsweb{aqstmp{airdata{download files.html. The dataset contains daily measurements on, among others, maximum temperature, and ground-level ozone, carbon monoxide and particulate matter concentrations, for the time period 1999 to 2013. These data are collected at stations spread over the U.S. We focus the analysis on the ground-level ozone and particulate matter concentrations. In order to estimate the extremal dependence between these, we calculate the component-wise monthly maximum of daily maximum concentrations, and estimate the Pickands dependence function conditional on the covariates time and location, where the latter is expressed by latitude and longitude. The estimation method was implemented with the same cross validation criteria as in the simulation section, including the same choices for C and H, after standardising the covariates to the interval r0, 1s. As kernel function K ¦ we use the following generalisation of the bi-quadratic kernel

K : K ¦ px 1 , x 2 , x 3 q : 3 ¹ i1 Kpx i q,
where x 1 , x 2 , x 3 , refer to the covariates time, latitude and longitude, respectively, in standardised form.

Note that K ¦ has as support the unit ball with respect to the max-norm on R 3 . We report here only the results for the city of Houston. Similar results can though be obtained for other cities or regions in the U.S.

In the left panel of Figure 2, we show the time plot of the estimates for the conditional extremal coefficient over the observation period. The conditional extremal coefficient is defined as ηpxq 2A 0 p0.5|xq, and is often used as a summary measure of extremal dependence. Its range is r1, 2s, where 1 corresponds with perfect dependence and 2 with independence. The time plot shows a seasonal pattern in the extremal dependence, and moreover the extremal dependence seems to decrease with time. We also observe that the estimates for α 0 and 0.1 are similar, but different from those obtained with α 0.5 and 1 (which are also similar), indicating that the dataset contains contamination with respect to the dependence structure. In order to get a better idea about the extremal dependence we show in the right panel of Figure 2 the estimate of A 0 pt|xq for a particular month (April 2002). This plot shows again estimates which are similar for α 0 and 0.1, but different from those obtained with α 0.5 and 1 (which are similar), confirming our earlier observation that there are observations which are contaminating with respect to the dependence structure.

Appendix A: Proofs of the results

A.1. Proof of Lemma 2.1

Using the fact that for any x R p the conditional distribution function of Z t given X x is an exponential distribution with parameter A 0 pt|xq and since λa A 0 pt|xq ¡ 0, we have Then

E Z β t e ¡λaZt § § §X x % Γpβ 1q A 0 pt|xq pλa A 0 pt|xqq β 1 . (A.
E K h px 0 ¡ XqZ β t e ¡λaZt % E K h px 0 ¡ XqΓpβ 1q A 0 pt|Xq pλa A 0 pt|Xqq β 1 & Γpβ 1q » R p K h px 0 ¡ yq A 0 pt|yq pλa A 0 pt|yqq β 1 f pyqdy Γpβ 1q » S K Kpzq A 0 pt|x 0 ¡ zhq pλa A 0 pt|x 0 ¡ zhqq β 1 f px 0 ¡ hzqdz Γpβ 1q
A 0 pt|x 0 q pλa A 0 pt|x 0 qq β 1 f px 0 qp1 op1qq.

(A.2)

The last transition in the above display follows since z S K and for n large enough, using the continuity of A 0 pt|.q and f at x 0 IntpS X q non-empty, we have boundedness in a neighborhood of x 0 , allowing us to use Lebesgue's dominated convergence theorem. Consequently

ErT n pK, a, t, λ, β, γ|x 0 qs a γ Γpβ 1q

A 0 pt|x 0 q pλa A 0 pt|x 0 qq β 1 f px 0 qp1 op1qq.

Also, similar arguments yield

VarpT n pK, a, t, λ, β, γ|x 0 qq 1 nh p }K} 2 2 A 0 pt|x 0 q a 2γ Γp2β 1q f px 0 q p2λa A 0 pt|x 0 qq 2β 1 p1 op1qq op1q, from which the convergence in probability simply follows.

A.2. Asymptotic covariance matrix of the finite dimensional vector T n

Our aim in this section is to compute the explicit expression of the elements of the covariance matrix Σ pσ j,k q 1¤j,k¤m given in (2.6). In this section we work under the assumptions of Lemma 2.1. According to (A.2), we have

E K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j %
f px 0 qΓpβ j 1q rλ j 1s βj 1 rA 0 pt j |x 0 qs βj p1 op1qq, for 1 ¤ j ¤ m. In order to compute the cross expectation, we need to derive the conditional distribution function of the pair pZ tj , Z t k q given X x. Let u, v ¡ 0

P Zt j ¡ u, Zt k ¡ v|X x ¨ P ¡ Y p1q ¡ max pp1 ¡ tjqu, p1 ¡ t k qvq , Y p2q ¡ max ptju, t k vq § § §X x © G ¡ max pp1 ¡ tjqu, p1 ¡ t k qvq , max ptju, t k vq § § §x © .
Hence, for j, k t1, ..., mu 2 , we have

E K h px0 ¡ XqZ β j t j e ¡λ j A 0 pt j |x 0 qZ t j K h px0 ¡ XqZ β k t k e ¡λ k A 0 pt k |x 0 qZ t k % (A.3) E K 2 h px0 ¡ XqE Z β j t j e ¡λ j A 0 pt j |x 0 qZ t j Z β k t k e ¡λ k A 0 pt k |x 0 qZ t k § § § X %% .
We focus now on the conditional expectation. Using (A.1) and the fact that

z β e ¡aλz ¡ δ 0 pβq » R 1l tz¡uu u β¡1 pβ ¡ aλuqe ¡aλu du, (A.4)
we have

E Z βj tj e ¡λjA0ptj|x0qZt j Z β k t k e ¡λ k A0pt k |x0qZt k § § § X % E ¡ Z βj tj e ¡λjA0ptj|x0qZt j ¡ δ 0 pβ j q © ¡ Z β k t k e ¡λ k A0pt k |x0qZt k ¡ δ 0 pβ k q © § § § X % ¡ δ 0 pβ j qδ 0 pβ k q δ 0 pβ j qE Z β k t k e ¡λ k A0pt k |x0qZt k § § § X % δ 0 pβ k qE Z βj tj e ¡λjA0ptj|x0qZt j § § § X % » R 2 g j,k pu, vqG tj ,t k pu, v|Xqdudv ¡ δ 0 pβ j qδ 0 pβ k q δ 0 pβ j q Γpβ k 1q rλ k 1s β k 1 rA 0 pt k |Xqs β k δ 0 pβ k q Γpβ j 1q
rλ j 1s βj 1 rA 0 pt j |Xqs βj .

(A.5)

Combining the continuity at x 0 and boundedness of the functions f , A 0 pt|.q and Gpu, v|.q, the expression of σ j,k in (2.6) follows.

A.3. Proof of Theorem 2.1

First, remark that to show Theorem 2.1, it is sufficient to look at the weak convergence of the process

3 c nh p pT n pK, A 0 pt|x 0 q, t, λ, β, γ|x 0 q ¡ E rT n pK, A 0 pt|x 0 q, t, λ, β, γ|x 0 qsq , t r0, 1s A , (A.6) since lim nÑV sup tr0,1s c nh p § § § § E rT n pK, A 0 pt|x 0 q, t, λ, β, γ|x 0 qs ¡ Γpβ 1q rA 0 pt|x 0 qs γ¡β pλ 1q β 1 f px 0 q § § § § 0.
Indeed, according to (A.2), we have

§ § § § E rT n pK, A 0 pt|x 0 q, t, λ, β, γ|x 0 qs ¡ Γpβ 1q rA 0 pt|x 0 qs γ¡β pλ 1q β 1 f px 0 q § § § § ¤ Γpβ 1qA γ 0 pt|x 0 q » S K Kpyq § § § § § A 0 pt|x 0 ¡ yhq pλA 0 pt|x 0 q A 0 pt|x 0 ¡ yhqq β 1 f px 0 ¡ hyq ¡ A ¡β 0 pt|x 0 q pλ 1q β 1 f px 0 q § § § § § dy.
Now, using Assumptions pDq and pA 0 q, we deduce that §

§ § § § A0pt|x0 ¡ yhq pλA0pt|x0q A0pt|x0 ¡ yhqq β 1 f px0 ¡ hyq ¡ A ¡β 0 pt|x0q pλ 1q β 1 f px0q § § § § § ¤ A0pt|x0 ¡ yhq pλA0pt|x0q A0pt|x0 ¡ yhqq β 1 |fpx0 ¡ yhq ¡ f px0q| § § § § § A0pt|x0 ¡ yhq pλA0pt|x0q A0pt|x0 ¡ yhqq β 1 ¡ A ¡β 0 pt|x0q pλ 1q β 1 § § § § § f px0q Oph minpη f ,η A 0 q q
for n large enough such that h ¤ 1, with a bound which is uniform in t.

Then, to show the weak convergence of the stochastic process (A.6), we will use Theorem 19.28 in van der Vaart (1998). To apply this result, we need to introduce some notations. Define the covering number N pF, L 2 pQq, τ q as the minimal number of L 2 pQq-balls of radius τ needed to cover the class of functions F and the uniform entropy integral as

Jpδ, F, L 2 q : » δ 0 log sup Q N pF, L 2 pQq, τ }F} Q,2 q dτ,
where Q is the set of all probability measures Q for which 0

}F} 2 Q,2 : ³ F 2 dQ
V and F is an envelope function of the class F.

Let P denote the law of the vector pY p1q , Y p2q , Xq and define the expectation under P , the empirical version and empirical process as follows

P f : » f dP, P n f : 1 n n i1 f ¡ Y p1q i , Y p2q i , X i © , G n f : c npP n ¡ P qf,
for any real-valued measurable function f .

For any γ R and pλ, βq p0, Vq ¢ R or pλ, βq p0, 0q, we introduce our sequence of classes F n as

Fn : tpy1, y2, zq Ñ fn,tpy1, y2, zq, t r0, 1su : 3 py1, y2, zq Ñ c h p K h px0 ¡ zqrA0pt|x0qs γ¡β rA0pt|x0qZtpy1, y2qs β e ¡λA 0 pt|x 0 qZ t py 1 ,y 2 q , t r0, 1s A , where Z t py 1 , y 2 q : min

¢ y 1 1 ¡ t , y 2 t . Remark that Z t Z t ¡ r Y p1q , r Y p2q ©
. Denote now by F n an envelope function of the class F n and for any y S X , define the bivariate function ρ y : r0, 1s 2 Ñ R as

ρ y pt, sq : E ¡ A γ 0 pt|x 0 qZ β t e ¡λA0pt|x0qZt ¡ A γ 0 ps|x 0 qZ β s e ¡λA0ps|x0qZs © 2 § § § § X y & .
Naturally, ρ y defines a semimetric on r0, 1s 2 and since it is bi-continuous, it makes r0, 1s totally bounded. Now, according to Theorem 19.28 in van der Vaart (1998), the weak convergence of the stochastic process (A.6) follows from the four following conditions sup ρx 0 pt,sq¤δn P pf n,t ¡ f n,s q 2 ÝÑ 0 for every δ n × 0, (A.7)

P F 2 n Op1q, (A.8) P F 2 n tF n ¡ ε c
nu ÝÑ 0 for every ε ¡ 0, (A.9)

Jpδ n , F n , L 2 q ÝÑ 0 for every δ n × 0.

(A.10)

We start to prove (A.7). By definition, we have Such expectations have been computed in (A.1) and (A.5). Using the mean value theorem combined with the boundedness of A 0 p.|.q and Assumption pA 0 q, we can easily infer that for all py, y I q B x0 prq ¢ B x0 prq, we have sup pt,sqr0,1s 2 |ρ y pt, sq ¡ ρ y Ipt, sq| ¤ C}y ¡ y I } η A 0 , for some positive constant C. This implies (A.11) and thus (A.7) is established. Now, we move to the proof of (A.8) and (A.9). Since the function x Ñ x β e ¡λx is bounded over R by pβ{λq β e ¡β and A 0 pt|x 0 q r1{2, 1s, F n admits the natural envelope function py 1 , y 2 , zq Ñ F n py 1 , y 2 , zq :

P pf n,t ¡ f n,s q 2 » R p h ¡p K 2 ¢ x 0 ¡ u h ρ u pt, sqf puqdu » S K K 2 puqρ x0¡hu pt, sqf px 0 ¡ huqdu }K} 2 2 f px 0 qρ x0 pt, sq » S K K 2 puqfpx 0 ¡ huqrρ x0¡hu pt, sq ¡ ρ x0 pt, sqsdu ρ x0 pt, sq » S K K 2 puqrfpx 0 ¡ huq ¡ f px 0 qsdu.
c h p K h px 0 ¡ zqM, (A.12)
where M :

¢ β λ β e ¡β maxp1, 2 β¡γ q. Consequently P F 2 n M 2 » R p h ¡p K 2 ¡ x0 ¡ u h © f puqdu M 2 » S K K 2 puqfpx0 ¡ huqdu M 2 }K} 2 2 f px0qp1 op1qq, P F 2 n tFn ¡ ε c nu M 2 » tKpuq¡M ¡1 ε c nh p u K 2 puqfpx0 ¡ huqdu 0,
for all ε ¡ 0 and n sufficiently large, since nh p Ñ V, K satisfies Assumption pK 1 q and f is continuous.

Finally, it remains to prove (A.10). First, we introduce the class of functions W : tpy 1 , y 2 q Ñ A 0 pt|x 0 qZ t py 1 , y 2 q, t r0, 1su and its subgraph σ t in R 2 ¢ R as σ t : tpu, v, wq : A 0 pt|x 0 qZ t pu, vq ¡ wu 4 pu, v, wq :

A 0 pt|x 0 q 1 ¡ t u ¡ w B 4 pu, v, wq : A 0 pt|x 0 q t v ¡ w B .
We can show that tσ t : t r0, 1su is a VC -class of sets. Indeed, if we look more generally, at the collection of sets C : ttpx, yq : δx ¡ yu, δ ¡ 0u in R ¢R and if we define two points px 1 , y 1 q and px 2 , y 2 q such that, without loss of generality,

y 1 x 1 ¤ y 2 x 2
. Then, for any δ ¡ 0, δx 2 ¥ y 2 implies that δx 1 ¥ y 1 . Thus, C cannot shatter the set tpx 1 , y 1 q, px 2 , y 2 qu and by consequence it is a VC -class of sets. Now, the collection of one set R is naturally a VC -class of sets. According to Lemma 2.6.17 (vii) in van der [START_REF] Van Der Vaart | Weak convergence and empirical processes, with applications to statistics[END_REF], C ¢ R is aVC -class of sets as well. Invoking Lemma 2.6.17 (ii), tσ t : t r0, 1su belongs to a VC -class and as such is VC. Define now for all z R φ λ,β pzq : z β e ¡λz .

We can easily check that φ λ,β is of bounded variation. This implies that φ λ,β can be decomposed as the sum of two monotone functions, say φ p1q λ,β and φ p2q λ,β . Thus, according to Lemma 2.6.18 (viii) in van der

Vaart and Wellner (1996), φ p1q

λ,β ¥ W and φ p2q λ,β ¥ W are VC. Now, according to Theorem 2.6.7 in van der [START_REF] Van Der Vaart | Weak convergence and empirical processes, with applications to statistics[END_REF], there exists a universal constant C such that for any j 1, 2 and 0 τ 1

sup Q N pφ pjq λ,β ¥ W, L 2 pQq, τ }W j } Q,2 q ¤ CV j p16eq Vj ¢ 1 τ 2pVj ¡1q
, where V j is the VC -index of φ pjq λ,β ¥ W and W j its envelope function. Now, consider the sequence of class of functions

F n,j : tz Ñ c h p K h px 0 ¡ zqu φ pjq λ,β ¥ W, for j 1, 2
, where denotes the direct product between the two classes involved. Since we only update the previous sets with one single function and only one ball is needed to cover the class tz Ñ c h p K h px 0 ¡zqu whatever the measure Q, we have sup

Q N pF n,j , L 2 pQq, τ }κF n } Q,2 q ¤ CV j p16eq Vj ¢ 1 τ 2pVj ¡1q
, where κ is a suitable constant. Moreover since sup tr0,1s rA 0 pt|x 0 qs γ¡β maxp1, 2 β¡γ q, for any 0 τ 1, the minimal number of balls of radius τ maxp1, 2 β¡γ q needed to cover the interval

0, maxp1, 2 β¡γ q $ is r1{2τ s. Hence sup Q N trA 0 pt|x 0 qs γ¡β , t r0, 1su, L 2 pQq, τ maxp1, 2 β¡γ q ¨ 1 2τ ¤ 3 2 ¢ 1 τ 2 . Consequently, we have sup Q N ptrA0pt|x0qs γ¡β , t r0, 1su Fn,j, L2pQq, τ maxp1, 2 β¡γ q}κFn}Q,2q ¤ 3C 2 Vjp16eq V j ¢ 1 τ 2V j .
Finally, since our class of interest F n is included in the class of functions r F n : trA 0 pt|x 0 qs γ¡β , t r0, 1su F n,1 trA 0 pt|x 0 qs γ¡β , t r0, 1su F n,2 , with envelope function 2 maxp1, 2 β¡γ qκF n , using Lemma 16 in [START_REF] Nolan | U -processes: rates of convergence[END_REF], we have sup

Q N pF n , L 2 pQq, 2τ maxp1, 2 β¡γ q}κF n } Q,2 q ¤ sup Q N p r F n , L 2 pQq, 2τ maxp1, 2 β¡γ q}κF n } Q,2 q ¤ 9C 2 4 V 1 V 2 p16eq V1 V2 ¢ 4 τ 2pV1 V2q : L ¢ 1 τ V .
Thus, (A.10) is established since for any sequence δ n × 0 and n large enough, we have

Jpδ n , F n , L 2 q ¤ » δn 0 logpr2κ maxp1, 2 β¡γ qs V Lq ¡ V logpτ qdτ op1q.
This achieves the proof of Theorem 2.1 since the covariance structure follows from (2.6).

A.4. Proof of Theorem 2.2

To prove this theorem, we will make use of the Cramér-Wold device (see, e.g., Severini, 2005, p. 337), according to which it is sufficient to show that

Λ n : ξ T c nh p pT n ¡ ErT n sq N 1 0, ξ T Σξ ¨,
for all ξ R m . A straightforward rearrangement of the terms leads to

Λn 1 n n i1 c nh p 5 m j1 ξjrA0ptj|x0qs γ j K h px0 ¡ XiqZ β j t j ,i e ¡λ j A 0 pt j |x 0 qZ t j ,i ¡E m j1 ξjrA0ptj|x0qs γ j K h px0 ¡ XiqZ β j t j ,i e ¡λ j A 0 pt j |x 0 qZ t j ,i 'C : 1 n n i1
Wi.

Since W 1 , ..., W n are independent and identically distributed random variables, VarpΛ n q VarpW 1 q

n with VarpW 1 q nh p m j1 m ķ1 ξ j ξ k C j,k , where C j,k : E pA 0 pt j |x 0 qq γj pA 0 pt k |x 0 qq γ k K 2 h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j Z β k t k e ¡λ k A0pt k |x0qZt k % ¡E pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j % E pA 0 pt k |x 0 qq γ k K h px 0 ¡ XqZ β k t k e ¡λ k A0pt k |x0qZt k % .
According to the computations in Appendix A.2, V arpΛ n q ξ T Σξp1 op1qq. Hence, to ensure the convergence in distribution of Λ n to a normal random variable, we have to verify the Lyapounov condition for triangular arrays of random variables (Billingsley, 1995, p. 362). In the present context this simplifies to verifying 1 n 2 Ep|W 1 | 3 q Ñ 0. We have

Ep|W 1 | 3 q ¤ n 3{2 h 3p{2 6 8 7 E ! £ m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j 3 ( ) 3E ! £ m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j 2 ( ) ¢E m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j ' 4 E £ m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j ' 3 D F E .
A similar treatment as for (A.3) yields for all postive integer q

E £ m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j q ' E £ E £ m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j q § § § § § X ' : E rK q h px 0 ¡ XqQpXqs
where the explicit expression of QpXq can be obtained similarly as for (A.5). Hence

E £ m j1 |ξ j |pA 0 pt j |x 0 qq γj K h px 0 ¡ XqZ βj tj e ¡λjA0ptj|x0qZt j q ' 1 h qp » R p K q ¢ x 0 ¡ u h Qpuqf puqdu ph p q 1¡q » S K K q pzqQpx 0 ¡ zhqf px 0 ¡ zhqdz Opph p q 1¡q q
by continuity and boundedness of the functions. Consequently

1 n 2 Ep|W 1 | 3 q O ¡ p c nh p q ¡1 © op1q.
A.5. The derivatives of p ∆ α,x0,t and their asymptotic properties

Straightforward computations for a r1{2, 1s, α ¡ 0, give p ∆ p1q α,x 0 ,t paq αa ¡1 p ∆α,x 0 ,tpaq a α p1 αq 1 Proof of Corollary A.1. As usual, it is sufficient to show the finite dimensional convergence and the tightness of the process. Using Theorem 2.2 we directly solve the finite dimensional convergence issue.

n n i1 K h px0 ¡ XiqZt,ie ¡αaZ t,i , p ∆ p2q α,x 0 ,t paq αa ¡1 p ∆ p1q α,x 0 ,t paq ¡ αa ¡2 p ∆α,x 0 ,tpaq αpα 1qa α¡1 1 n n i1 K h px0 ¡ Xiqp1 ¡ aZt,iqZt,ie ¡αaZ t,i , p ∆ p3q α,x 0 ,t paq α ¡ 2a ¡3 p ∆α,x 0 ,tpaq a ¡1 p ∆ p2q α,x 0 ,t paq ¡ 2a ¡2 p ∆ p1q α,x 0 ,t paq © pα ¡ 1qαpα 1q a α¡2 n n i1 K h px0 ¡ Xiqp1 ¡ aZt,iqZt,ie ¡αaZ t,i ¡αpα 1q a α¡1 n n i1 K h px0 ¡ Xiqpαp1 ¡ aZt,
Next, Theorem 2.1 combined with (A.6) implies tightness for any process t Ñ c nh p pT n ¡ErT n sqpK, A 0 pt|x 0 q, t, λ, β, γ|x 0 q and similarly as in Lemma 1 in [START_REF] Bai | Multivariate limit theorems in the context of long-range dependence[END_REF], we have tightness for any multivariate process with similar coordinates. Corollary A.1 then follows.

A.6. Proof of Theorem 2.3

To prove the theorem we will adjust the arguments used to prove existence and consistency of solutions of the likelihood estimating equation, see e.g. Theorem 3.7 and Theorem 5.1 in Chapter 6 of [START_REF] Lehmann | Theory of point estimation[END_REF], to the MDPD framework. Let ζ, b ¡ 0, Cp.|.q : r0, 1s ¢ S X Ñ r1{2 ¡ ζ, 1 ζs and dt r0, 1s, rptq : |A 0 pt|x 0 q ¡ Cpt|x 0 q|. Define in addition the b-level of r as T b : tt r0, 1s, rptq ¡ bu . We firstly show that for any b ¡ 0

P ¡ dt T b , p ∆ α,x0,t pA 0 pt|x 0 qq p ∆ α,x0,t pCpt|x 0 qq © Ñ 1, (A.13)
as n Ñ V, for any function Cp.|x 0 q different from but close enough to A 0 p.|x 0 q. By applying a Taylor series expansion, we have p ∆α,x 0 ,tpCpt|x0qq ¡ p ∆α,x 0 ,tpA0pt|x0qq pCpt|x 0 q ¡ A 0 pt|x 0 qq p ∆ p1q

α,x 0 ,t pA 0 pt|x 0 qq 1 2 pCpt|x 0 q ¡ A 0 pt|x 0 qq 2 p ∆ p2q

α,x 0 ,t pA 0 pt|x 0 qq 1 6 pCpt|x 0 q ¡ A 0 pt|x 0 qq 3 p ∆ p3q

α,x 0 ,t p r Cpt|x 0 qq, where r

Cpt|x 0 q is an intermediate value between Cpt|x 0 q and A 0 pt|x 0 q. According to Appendix A.5, as

n Ñ V sup tr0,1s § § § p ∆ p1q α,x0,t pA 0 pt|x 0 qq § § § sup tr0,1s § § § p ∆ p1q α,x0,t pA 0 pt|x 0 qq ¡ p1q α,x0,t pA 0 pt|x 0 qq § § § P ÝÑ0.
This convergence implies, that for all 0 ε ¤ b

2 P ¡ dt T b , rptq| p ∆ p1q α,x 0 ,t pA 0 pt|x 0 qq| ¤ r 3 ptq © ¥ P £ dt T b , | p ∆ p1q α,x 0 ,t pA 0 pt|x 0 qq| ¤ r 2 ptq, sup tr0,1s § § § p ∆ p1q α,x 0 ,t pA 0 pt|x 0 q § § § ¤ ε P £ sup tr0,1s § § § p ∆ p1q
α,x 0 ,t pA 0 pt|x 0 q § § § ¤ ε

ÝÑ1,

as n Ñ V. α,x0,t pA 0 pt|x 0 qq ¡ δ 1 r 2 ptq, with probability tending to 1.

Finally, since x Ñ x λ e ¡x is bounded dλ ¥ 1 on R and by Lemma 2.1

T n pK, a, t, 0, 0, 0|x 0 q 1 n n i1 K h px 0 ¡ X i q P ÝÑfpx 0 q, as n Ñ V, we have for any ε ¡ 0, n ¡1 °n i1 K h px 0 ¡ X i q ¤ f px 0 q ε with probability tending to 1. This implies that sup rptq δ 1 {p1 M {6q, (A.13) follows.

To complete the proof we adjust the line of argumentation of Theorem 3.7 in Chapter 6 of [START_REF] Lehmann | Theory of point estimation[END_REF]. Take 0 δ ζ and define the event S n pδq :

3 dt r0, 1s, p ∆ α,x0,t pA 0 pt|x 0 qq p ∆ α,x0,t pA 0 pt|x 0 q ¨δq A .

For υ S n pδq, since p ∆ α,x0,t paq is differentiable with respect to a, there exists r A α,n,δ pt|x 0 q pA 0 pt|x 0 q ¡ δ, A 0 pt|x 0 q δq where p ∆ α,x0,t paq achieves a local minimum, so p ∆ p1q

α,x0,t p r A α,n,δ pt|x 0 qq 0.

By (A.13), PpS n pδqq Ñ 1 for any small enough δ, and hence there exists a sequence δ n Ó 0, such that PpS n pδ n qq Ñ 1, as n Ñ V. Now, let p A α,n pt|x 0 q : r A α,n,δn pt|x 0 q if υ S n pδ n q and arbitrary otherwise. Since υ S n pδ n q implies p ∆ p1q α,x0,t p p A α,n pt|x 0 qq 0, we have that

P ¡ p ∆ p1q α,x0,t p p A α,n pt|x 0 qq 0 © ¥ P pS n pδ n qq Ñ 1,
as n Ñ V, which establishes the existence part. Note that the measurability of the local minimum can be verified in the same way as it is done in the framework of maximum likelihood estimation (see e.g. Serfling, 1980, p. 147).

Concerning now the uniform consistency of the solution sequence, note that for any ε ¡ 0 and n large enough such that δ n ¤ ε, we have

P £ sup tr0,1s § § § p Aα,npt|x0q ¡ A0pt|x0q § § § ¤ ε ¥ P £ sup tr0,1s § § § p Aα,npt|x0q ¡ A0pt|x0q § § § ¤ δn ¥ P pSnpδnqq Ñ 1,
as n Ñ V, whence the uniform consistency of the estimator sequence.

A.7. Proof of Theorem 2.4

The starting point is (2.4). According to Corollary A. n pt|x 0 q, where T p3q

n pt|x 0 q : ¤ ¥ T n pK, A 0 pt|x 0 q, t, 0, 0, α ¡ 1|x 0 q T n pK, A 0 pt|x 0 q, t, α, 0, α ¡ 1|x 0 q T n pK, A 0 pt|x 0 q, t, α, 1, α|x 0 q .

A.8. Proof of Lemma 3.1

We use the following decomposition

F n,j py|xq ¡ F j py|xq 1 p f n pxq 5 1 n n i1 K c px ¡ X i q1l tY pjq i ¤yu ¡ E K c px ¡ Xq1l tY pjq ¤yu $ C ¡ 1 p f n pxq 5 1 n n i1 K c px ¡ X i qE 1l tY pjq i ¤yu § § §Xi % ¡ E K c px ¡ Xq1l tY pjq ¤yu $ C 1 p f n pxq 5 1 n n i1 K c px ¡ X i q rF j py|X i q ¡ F j py|xqs C : 1 p f n pxq tT 1 py|xq ¡ T 2 py|xq T 3 py|xqu , where p f n pxq : 1 n n i1 K c px ¡ X i q
denotes the kernel density estimator of f .

We start by showing that, for some q ¡ 1,

sup py,xqR¢S X § § § § § n ¡1 n i1 Kcpx ¡ Xiq1l tY pjq i ¤yu ¡ E Kcpx ¡ Xq1l tY pjq ¤yu % § § § § § o P £ | log c| q nc p , (A.15) sup xS X § § § § § n ¡1 n i1 Kcpx ¡ Xiq ¡ E rKcpx ¡ Xqs § § § § § o P £ | log c| q nc p , (A.16) sup py,xqR¢S X § § § § § n ¡1 n i1 Kcpx ¡ XiqE 1l tY pjq i ¤yu § § §Xi % ¡ E Kcpx ¡ Xq1l tY pjq ¤yu % § § § § § o P £ | log c| q nc p . (A.17)
To this aim, let us introduce the class

G : 4 pu, vq Ñ K ¢ x ¡ v d 1l tu¤yu ; y R, x S X , d ¡ 0 B 4 K ¢ x ¡ . d ; x S X , d ¡ 0 B 4 1l t.¤yu ; y R B : G 1 G 2 .
Under Assumption pK 2 q, G 1 is a uniformly bounded VC -class of measurable functions (see e.g. [START_REF] Giné | Rates of strong uniform consistency for multivariate kernel density estimators[END_REF]. Next, since the collection of all cells tp¡V, as, a Ru is a VC -class of sets, it follows that G 2 is also a uniformly bounded VC -class of measurable functions. Now, using the fact that the covering number of the direct product of two VC -classes is bounded by the product of the respective covering numbers,

G n : 4 pu, vq Ñ K ¢ x ¡ v c 1l tu¤yu ; y R, x S X , c c n ¡ 0 B ,
admits the same bound for the covering number as G, that is

N pG n , L 2 pQq, τ }K} V q ¤ C V G p16eq V G ¢ 1 τ 2pV G ¡1q : ¢ A G τ ν G ,
where C is a universal constant, τ p0, 1q and V G is the VC -index of G (see Theorem 2.6.7 in van der [START_REF] Van Der Vaart | Weak convergence and empirical processes, with applications to statistics[END_REF]. Now, according to Proposition 2.1 in [START_REF] Giné | On consistency of kernel density estimators for randomly censored data: rates holding uniformly over adaptive intervals[END_REF] (see also Theorem 2.1 in [START_REF] Giné | Rates of strong uniform consistency for multivariate kernel density estimators[END_REF] for σ 2 ¥ sup gGn Varpgq, U ¥ }K} V and 0 σ ¤ U , there exists

a universal constant B such that E sup py,xqR¢S X § § § § § n ¡1 n i1 K c px ¡ X i q1l tY pjq i ¤yu ¡ E K c px ¡ Xq1l tY pjq ¤yu $ § § § § § ' ¤ rnc p s ¡1 B U ν G log ¢ U A G σ d ν G nσ 2 log ¢ U A G σ ' . Since Var ¢ K ¢ x ¡ X c 1l tY pjq ¤yu ¤ c p » S K K 2 puq f px ¡ cuqdu ¤ c p }f} V }K} 2 2 ,
the choices σ 2 σ 2 n : c p }f} V }K} 2 2 and U }K} V imply that σ 2 n ¤ U 2 for n large enough. This yields (A.15). Similar arguments can be used in order to show (A.16). Also (A.17) can be shown similarly, though with a refinement as used in [START_REF] Portier | On the weak convergence of the empirical conditional copula under a simplifying assumption[END_REF], p. 23.

As for p

f n pxq we use (A.16), and obtain p f n pxq Ep p f n pxqq o P p | log c| q {pnc p qq, where the o P term is uniform in x S X . By using the assumptions on K, f , and (3.5) we derive, for n sufficiently large, the following uniform lower bound

E ¡ p f n pxq © » tzS K :x¡czS X u Kpzqf px ¡ zcqdz ¥ b » tzB0pδq:x¡czS X u Kpzqdz ¥ b m λ ptz B 0 pδq : x ¡ cz S X uq ¥ b m δ p ε. Whence tT 1 py|xq ¡ T 2 py|xqu{ p f n pxq o P p | log c| q {pnc p qq, uniformly in py, xq R ¢ S X .
Concerning T 3 py|xq we obtain for x S X the following direct bound |T 3 py|xq| p

f n pxq ¤ 1 p f n pxq 6 8 7 1 n i:}x¡Xi}¤c K c px ¡ X i q |F j py|X i q ¡ F j py|xq| D F E ¤ M Fj c η F j .
Combining the above results establishes the lemma.

A.9. Proof of Theorem 3.1

Let

I n : tg θ,δ,n : θ Θ, δ Hu where for θ : pt, aq Θ : r0, 1s ¢ r1{2, 1s, and δ H : Remark now that Proof of Assertion 2. The idea is to apply Lemma 2.2 in van der [START_REF] Van Der Vaart | Empirical processes indexed by estimated functions, Asymptotics: Particles, Processes and Inverse Problems[END_REF]. To this aim, first observe that the class E n pθ, bq admits an envelope function E n of the same form as F n in (A.12), for some suitable constant M ¡ 0. Thus E n satisfies the conditions (A.8) and (A.9), with F n replaced by One can check that the proof of Theorems 2.3 and 2.4 are mainly due to the asymptotic properties of p ∆ pjq α,x0,t , j 1, 2 and 3. Thus, if we are able to prove that the two key statistics T n and q

2 δ pδ 1 , δ 2 q; δ : R ¢ R ¢ S X Ñ R 2 @ , g θ,δ,n py 1 , y 2 , uq : c h p K h px 0 ¡ uqq θ,δ py 1 , y 2 , uq : c h p K h px 0 ¡ uqa γ rZ θ,
T n are sufficiently close enough, in the sense that sup tr0,1s,ar1{2,1s c nh p § § § q

T n ¡ T n § § § pK, a, t, λ, β, γ|x 0 q o P p1q, T n ¡ T n § § § % pK, a, t, λ, β, γ|x 0 q op1q, (A.24) then we can swap p ∆ pjq α,x0,t by q ∆ pjq α,x0,t , j 1, 2 and 3. According to Theorem 3.1, (A.23) is a direct consequence of (A.24). So it remains to prove (A.24). Note that

c nh p E § § § q T n ¡ T n § § § % pK, a, t, λ, β, γ|x 0 q c n E § § § § § 1 n n i1 c h p K h px 0 ¡ X i qa γ q Z β n,

  Fig 1.MISEpε, α|xq as a function of ε t0, 0.025, 0.05, ..., 0.2u with α 0 (black), α 0.1 (blue), α 0.5 (green) and α 1 (red). First type of contamination: rows 1 pn 1000q and 2 pn 5000q, x 0.1, 0.3, 0.5 from the left to the right. Second type of contamination: rows 3 pn 1000q and 4 pn 5000q, x 0.5, 0.7, 0.9 from the left to the right.

Fig 2 .

 2 Fig 2. Air pollution data: time plot of the estimate for the conditional extremal coefficient (left) and estimate for the conditional Pickands dependence function in April 2002 (right), α 0 (black), α 0.1 (blue), α 0.5 (green) and α 1 (red).

3 e

 3 F 1 r n δ ¦ 1,θ,bn § § , F 1 q, maxp § § F 1 r n δ ¦ 1,θ,bn § § , F 1 q $ A 2 e ¡ts minp § § F 2 r n δ ¦ 2,θ,bn § § , F 2 q, maxp § § F 2 r n δ ¦ 2,θ,bn § § , F 2 q $@ ¡p1¡tqs rF 1 ¡ r n b n , F 1 r n b n s A 2 e ¡ts rF 2 ¡ r n b n , F 2 r n b n s @ : A n,1 psq A n,2 psq.Since for any subsets A and B we have 1l tABu ¤ 1l tAu 1l tBu , we can deduce thatP ¡ s rminpZ θ,δ0 rnδ ¦ θ,bn , Z θ,δ0 q, maxpZ θ,δ0 rnδ ¦ θ,bn , Z θ,δ0 qs|X x 0 ¡ hu © ¤ P pA n,1 psq|X x 0 ¡ huq P pA n,2 psq|X x 0 ¡ huq rv¡rnbn,v rnbnsu dv ¤ 2r n b n 2r n b n 4r n b n . |q θ,δ0 rnδ ¦ θ,bn ¡ q θ,δ0 | § § §X x 0 ¡ hu % ¤ 4 c nh p r n b n sup ar1{2,1s » V 0 a γ |β ¡ λas|s β¡1 e ¡λas ds.This achieves the proof of Assertion 1 since K is bounded, sup ar1{2,1s ³ V 0 a γ |β ¡ λas|s β¡1 e ¡λas ds V, c nh p r n Ñ 0 and b n Ñ 0.

  Jpd n , tG n pθ, bq : θ Θu , L 2 q ÝÑ 0 for all d n × 0.(A.22) We start to show (A.21). SinceP G 2 n pθ, bq » S K K 2 puqE ¡ |q θ,δ0 rnδ ¦ θ,b ¡ q θ,δ0 | 2 § § §X x 0 ¡ hu © f px 0 ¡ huqdu,and (A.19), (A.21) follows from the proof of Assertion 1. Now, to deal with the uniform entropy integral, we can adjust the lines of proof of Theorem 2.1 by considering the classes of functions defined on R¢ R ¢ S X φ pjq λ,β ¥ W ¥ Ψ, j 1, 2,where Ψ is either the functionpy 1 , y 2 , uq Ñ p¡ logpF 1 py 1 |uqq, ¡ logpF 2 py 2 |uqq or py 1 , y 2 , uq Ñ ¡ log § § F 1 py 1 |uq r n δ ¦ 1,θ,b py 1 , y 2 , uq § § ¨, ¡ log § § F 2 py 2 |uq r n δ ¦ 2,θ,b py 1 , y 2 , uq § § ¨ẅhichare VC -classes. This allows us to prove that there exist positive constants C and V such thatsup Q N ptG n pθ, bq : θ Θu , L 2 pQq, τ }E n } Q,2 q ¤ C ¢ 1 τ V , from which (A.22) follows. This achieves the proof of Theorem 3.1.A.10. Proof of Theorem 3.2

c

  nP G n pθ, bq op1qby Assertion 1 since it is clear from its proof that b n Ñ 0 can be replaced by any fixed value b in(A.20) without changing the conclusion. This achieves the proof of Theorem 3.2.

Table 1

 1 First type of contamination -coverage probabilities of 90% confidence intervals.

				t 0.3			t 0.5			t 0.7
	x 0.1	α	0	0.1	0.5	1	0	0.1	0.5	1	0	0.1	0.5	1
		ε 0.0 0.95 0.95 0.96 0.98 0.95 0.96 0.96 0.97 0.95 0.96 0.96 0.98
	n 1000 ε 0.1 0.69 0.77 0.92 0.93 0.50 0.64 0.88 0.91 0.67 0.78 0.92 0.93
		ε 0.2 0.26 0.41 0.83 0.91 0.11 0.17 0.64 0.79 0.24 0.42 0.81 0.90
		ε 0.0 0.98 0.97 0.96 0.97 0.96 0.97 0.97 0.98 0.95 0.96 0.97 0.97
	n 5000 ε 0.1 0.14 0.28 0.84 0.90 0.06 0.12 0.69 0.82 0.13 0.29 0.84 0.90
		ε 0.2 0.02 0.05 0.55 0.83 0.00 0.01 0.16 0.51 0.01 0.03 0.56 0.80
				t 0.3			t 0.5			t 0.7
	x 0.3	α	0	0.1	0.5	1	0	0.1	0.5	1	0	0.1	0.5	1
		ε 0.0 0.97 0.98 0.93 0.93 0.96 0.96 0.95 0.93 0.94 0.95 0.95 0.94
	n 1000 ε 0.1 0.70 0.80 0.93 0.96 0.71 0.80 0.93 0.94 0.70 0.82 0.95 0.95
		ε 0.2 0.29 0.45 0.86 0.94 0.20 0.40 0.78 0.85 0.27 0.47 0.84 0.92
		ε 0.0 0.96 0.95 0.94 0.94 0.94 0.94 0.93 0.94 0.96 0.95 0.94 0.94
	n 5000 ε 0.1 0.14 0.33 0.87 0.93 0.15 0.35 0.81 0.86 0.17 0.41 0.85 0.92
		ε 0.2 0.01 0.04 0.57 0.83 0.01 0.02 0.41 0.63 0.01 0.03 0.55 0.83
				t 0.3			t 0.5			t 0.7
	x 0.5	α	0	0.1	0.5	1	0	0.1	0.5	1	0	0.1	0.5	1
		ε 0.0 0.97 0.99 0.98 0.97 0.97 0.97 0.97 0.95 0.96 0.97 0.97 0.98
	n 1000 ε 0.1 0.77 0.85 0.91 0.95 0.76 0.84 0.91 0.94 0.79 0.83 0.95 0.94
		ε 0.2 0.51 0.69 0.93 0.95 0.53 0.66 0.89 0.92 0.47 0.65 0.92 0.96
		ε 0.0 0.93 0.95 0.95 0.96 0.90 0.92 0.91 0.94 0.94 0.96 0.94 0.94
	n 5000 ε 0.1 0.30 0.53 0.91 0.95 0.39 0.56 0.90 0.91 0.31 0.54 0.93 0.94
		ε 0.2 0.04 0.10 0.70 0.87 0.06 0.11 0.66 0.80 0.05 0.08 0.69 0.86

Table 2

 2 Second type of contamination -coverage probabilities of 90% confidence intervals.

				t 0.3			t 0.5			t 0.7
	x 0.5	α	0	0.1	0.5	1	0	0.1	0.5	1	0	0.1	0.5	1
		ε 0.0 0.96 0.96 0.97 0.96 0.94 0.93 0.96 0.96 0.96 0.97 0.96 0.96
	n 1000 ε 0.1 0.97 0.96 0.96 0.95 0.91 0.93 0.94 0.93 0.96 0.97 0.98 0.96
		ε 0.2 0.99 0.99 0.96 0.96 0.80 0.88 0.92 0.94 0.99 0.99 0.96 0.95
		ε 0.0 0.93 0.94 0.95 0.96 0.93 0.93 0.95 0.97 0.95 0.96 0.95 0.96
	n 5000 ε 0.1 0.93 0.98 0.96 0.96 0.48 0.69 0.94 0.96 0.95 0.99 0.99 0.98
		ε 0.2 0.92 0.98 0.96 0.95 0.20 0.33 0.88 0.96 0.87 0.97 0.96 0.94
				t 0.3			t 0.5			t 0.7
	x 0.7	α	0	0.1	0.5	1	0	0.1	0.5	1	0	0.1	0.5	1
		ε 0.0 0.97 0.97 1.00 1.00 0.94 0.94 0.96 0.99 0.94 0.96 1.00 1.00
	n 1000 ε 0.1 0.98 0.98 1.00 1.00 0.76 0.84 0.96 0.98 0.96 0.95 1.00 1.00
		ε 0.2 0.98 0.98 1.00 1.00 0.55 0.70 0.94 0.99 0.98 0.98 1.00 1.00
		ε 0.0 0.92 0.93 0.94 0.94 0.92 0.94 0.95 0.95 0.94 0.95 0.95 0.96
	n 5000 ε 0.1 0.81 0.89 0.95 0.96 0.24 0.54 0.94 0.95 0.79 0.88 0.95 0.95
		ε 0.2 0.57 0.67 0.90 0.94 0.06 0.11 0.77 0.90 0.60 0.69 0.91 0.93
				t 0.3			t 0.5			t 0.7
	x 0.9	α	0	0.1	0.5	1	0	0.1	0.5	1	0	0.1	0.5	1
		ε 0.0 1.00 1.00 1.00 1.00 0.99 0.99 1.00 1.00 0.98 0.98 0.99 0.99
	n 1000 ε 0.1 0.91 0.95 0.99 0.99 0.64 0.81 0.97 0.99 0.91 0.93 0.99 0.99
		ε 0.2 0.89 0.92 0.98 1.00 0.40 0.63 0.96 0.98 0.89 0.93 0.99 1.00
		ε 0.0 0.98 0.98 0.99 0.98 0.98 0.97 0.98 0.98 0.97 0.97 0.98 0.98
	n 5000 ε 0.1 0.61 0.77 0.97 0.98 0.15 0.39 0.95 0.97 0.60 0.82 0.98 0.99
		ε 0.2 0.26 0.44 0.92 0.95 0.01 0.09 0.74 0.91 0.24 0.42 0.91 0.96

  By the Assumptions pDq, pK 1 q and since ρ x0 is bounded, it remains to show that sup ρx 0 pt,sq¤δn |ρ x0¡hu pt, sq ¡ ρ x0 pt, sq| Ñ 0.

		(A.11)
	Recall that
	ρ y pt, sq rA 0 pt|x 0 qs 2γ E ¡2 rA 0 pt|x 0 qA 0 ps|x 0 qs γ Z 2β t e ¡2λA0pt|x0qZt E Z β t e ¡λA0pt|x0qZt Z β § § §X y % rA 0 ps|x 0 qs 2γ E s e ¡λA0ps|x0qZs § § §X y Z 2β s e ¡2λA0ps|x0qZs % .	§ § §X y	%

  iq 1qZ 2 t,i e ¡αaZ t,i . pA 0 pt|x 0 qq : 1 α 2 p1 αq 2 rA 0 pt|x 0 qs α¡2 f px 0 q. pA 0 pt|x 0 qq, j t1, 2u, to its limit is also useful to study (2.4) and thus to reach our final goal. The aim of the next corollary is to provide such a rate.

	The convergence in probability of the two first derivatives of p ∆ α,x0,t is therefore a direct application of Lemma 2.1, which yields as n Ñ V
	p ∆ p1q α,x0,t pA 0 pt|x 0 qq P ÝÑ p1q α,x0,t pA 0 pt|x 0 qq : 0,
	p ∆ p2q α,x0,t pA 0 pt|x 0 qq P ÝÑ p2q α,x0,t Now the rate of convergence of p ∆ pjq α,x0,t Corollary A.1. Under the assumptions of Theorem 2.1, then for any j t1, 2u, the process 3 c nh p ¡ p ∆ pjq α,x0,t pA 0 pt|x 0 qq ¡ pjq α,x0,t pA 0 pt|x 0 qq © , t r0, 1s A
	weakly converges in Cpr0, 1sq towards a tight centered Gaussian process. In particular, we have
	sup	§ § § p ∆ pjq α,x0,t pA 0 pt|x 0 qq ¡ pjq α,x0,t pA 0 pt|x 0 qq	§ § § o P p1q.
	tr0,1s

  n Ñ V, where the right-hand side of the inequality is positive for rptq δ 1 {p1 M {6q.

	with probability tending to 1.
	Overall, we have shown that
	P	¢	dt T b , p ∆ α,x0,t pCpt|x 0 qq ¡ p ∆ α,x0,t pA 0 pt|x 0 qq ¡ δ 1 r 2 ptq ¡	¢	1	M 6	r 3 ptq	ÝÑ1,
									Thus, setting
			sup
			tr0,1s
				§ § § p ∆ p3q α,x0,t paq § § § : M V			(A.14)
			ar1{2¡ζ,1 ζs,tr0,1s
	with probability tending to 1. We can therefore conclude that
			dt r0, 1s, r 3 ptq 6	§ § § p ∆ p3q α,x0,t p r Cpt|x 0 qq	§ § § ¤ M 6	r 3 ptq,

as

  Apt|x 0 qqp p A α,n pt|x 0 q ¡ A 0 pt|x 0 qq

	Concerning the covariance structure, it follows from Theorem 2.2 and the fact that
				p ∆ p1q α,x0,t pA 0 pt|x 0 qq v T α T p3q
	1, verges, as n Ñ V, towards a tight centered Gaussian process and 3 c nh p p ∆ p1q α,x0,t pA 0 pt|x 0 qq, t r0, 1s	A	weakly con-
		3	p ∆ p2q α,x0,t pA 0 pt|x 0 qq , t r0, 1s	A	P ÝÑ	3 p2q α,x0,t pA 0 pt|x 0 qq, t r0, 1s	A	.
	Combining these results with (A.14), we have, as n Ñ V,
	5	p ∆ p2q α,x0,t pA 0 pt|x 0 qq 1 2	p ∆ p3q α,x0,t p r				& ¡1	, t r0, 1s	C
				P ÝÑ	4 p2q α,x0,t pA 0 pt|x 0 qq	% ¡1	, t r0, 1s	B	.

.

  δ py 1 , y 2 , uqs β exp p¡λaZ θ,δ py 1 , y 2 , uqq with Z θ,δ py 1 , y 2 , uq : min ¢ ¡ log p|δ 1 py 1 , y 2 , uq|q 1 ¡ t , ¡ log p|δ 2 py 1 , y 2 , uq|q t For convenience, denote δ n : pF n,1 , F n,2 q and δ 0 : pF 1 , F 2 q. According to Lemma 3.1, r ¡1 n |δ n ¡ δ 0 | converges in probability towards the null function H 0 : t0u in H, endowed with the norm }δ} H : }δ 1 } V }δ 2 } V for any δ H. In order to apply Theorem 2.3 in van der Vaart and Wellner (2007), we have now to show Assertion 1: sup θΘ c nP G n pθ, b n q ÝÑ 0 for every b n Ñ 0 and Assertion 2: sup θΘ |G n G n pθ, bq| P ÝÑ 0, for every b ¡ 0, where G n pθ, bq is the minimal envelope function for the class E n pθ, bq : tg θ,δ0 rnδ,n ¡ g θ,δ0,n : δ H, }δ} H ¤ bu , θ,δ0 rnδ ¡ q θ,δ0 |py 1 , y 2 , uq, where B : tpx, yq R 2 : |x| |y| ¤ bu. Since B is compact and δ Ñ q θ,δ py 1 , y 2 , uq is continuous, (A.18) reaches its supremum on at least one position δ ¦ θ,b py 1 , y 2 , uq pδ ¦ 1,θ,b py 1 , y 2 , uq, δ ¦ 2,θ,b py 1 , y 2 , uqq in B. Thus, according to Theorem 18.19 in Aliprantis and Border (2006), one can find a measurable function δ ¦ θ,b bounded by b in H such that G n pθ, bq |g θ,δ0 rnδ ¦ θ,b ,n ¡ g θ,δ0,n |. Proof of Assertion 1. For any positive sequence b n Ñ 0, we have c nP G n pθ, b n q cNote that for any pδ, δ I q H ¢ H, using (A.4) |q θ,δ ¡ q θ,δ I| ¤ a γ λas|s β¡1 e ¡λas 1l tsrminpZ θ,δ ,Z θ,δ Iq,maxpZθ,δ,Z θ,δ Iqsu ds.|q θ,δ 0 rnδ ¦ θ,bn ¡ q θ,δ 0 | § § §X x0 ¡ hu Z θ,δ 0 q, maxpZ θ,δ 0 rnδ ¦

	Consequently
								%
	E
	¤ a γ	» V 0	|β ¡ λas|s β¡1 e ¡λas	P	¡	s rminpZ θ,δ 0 rnδ ¦ θ,bn	θ,bn	, Z θ,δ 0 qs	§ § §X x0 ¡ hu	©	ds.
	i.e.
			G n pθ, bq : sup }δ} H ¤b c h p K h px 0 ¡ ¤q sup |g θ,δ0 rnδ,n ¡ g θ,δ0,n | }δ} H ¤b |q θ,δ0 rnδ ¡ q θ,δ0 |.
								sup pδ1py1,y2,uq,δ2py1,y2,uqqB
			nh p	» S K	KpuqE	|q θ,δ0 rnδ ¦ θ,bn ¡ q θ,δ0 | § § §X x 0 ¡ hu %	f px 0 ¡ huqdu.
											(A.19)

(A.18) Now, remark that dpy 1 , y 2 , uq R ¢ R ¢ S X sup }δ} H ¤b |q θ,δ0 rnδ ¡ q θ,δ0 |py 1 , y 2 , uq |q » V 0 |β ¡ ,

  t,i e ¡λa q K h px 0 ¡ X i qa γ Z β t,i e ¡λaZt,i δ n,n pY p1q , Y p2q , Xq ¡ g θ,δ0,n pY p1q , Y p2q , Xq § § § % ¤ cn P G n pθ, bq, since δ n δ 0 r n Bp0, bq where Bp0, bq : tδ : }δ} H ¤ bu. This implies that

			Zn,t,i
	¤ c n E §gθ,sup § § c nh p E § § § q T n ¡ T n § § § % pK, a, t, λ, β, γ|x 0 q ¤	¡	c h p % § % § § sup
	tr0,1s,ar1{2,1s	tr0,1s,ar1{2,1s

Acknowledgements

The authors sincerely thank the editor, associate editor and the referees for their helpful comments and suggestions that led to substantial improvement of the paper.

Supplementary material

Supplementary material to "Local robust estimation of the Pickands dependence function". This document contains additional simulation results.