Farah Gabsi
email: farah.gabsi@univ-lorraine.fr

Frédéric Hamelin
email: frederic.hamelin@univ-lorraine.fr

Nathalie Sauer
email: nathalie.sauer@univ-lorraine.fr

Building hygrothermal modeling by nodal method

Keywords: buildings, hygrothermal modeling, thermal comfort, nodal analysis

.

I. INTRODUCTION

The park of existing buildings in France has 2/3 of dwellings built before 1975. The renewal rate of this park is less than 1% per year. It therefore represents a significant potential for achieving the objectives of "Factor 4", validated by the Grenelle Environnement in 2007, which aims to divide by 4 greenhouse gas emissions by 2050. The impact of mass transfers on the real energy behavior of buildings has been highlighted for several years [START_REF] Woloszyn | Tools for performance simulation of heat, air and moisture conditions of whole buildings[END_REF]. The transfer equations are known, but very few works [START_REF] Peuhkuri | Moisture dynamics in building envelopes[END_REF], [START_REF] Xu | A humidity integrated building thermal model[END_REF] integrate these models in order to implement energy efficient control laws that ensure hygrothermal comfort [START_REF] Djongyang | Thermal comfort: A review paper[END_REF] within buildings. Yet, moisture directly impacts comfort in buildings since its rate reflects both the health of buildings and the quality of indoor air. The main causes of development of moisture [START_REF] Moon | Evaluation of hygrothermal models for mold growth avoidance prediction[END_REF] are related to defects in HVAC systems (direct infiltration of moist air, negative pressure gradient across the envelope, inadequate ventilation, etc.) or inappropriate use of the building (high occupancy density, bad air conditioning settings, etc.).

In this paper, a complete methodology of modeling the hygrothermal behavior of a multi-zone office building is proposed. It uses the Simscape TM library in MATLAB/Simulink R environment. The application target is the "Eco-Safe" platform (CRAN Nancy, France) that fits in a global context of smart home [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF].

II. THERMAL COMFORT

The notion of comfort is closely related to one-off adjustments in body regulation, both global and local. The perception of cold and heat depends on many parameters. However, it is possible to statistically characterize average This work has financial support from the Contrat de Plan Etat-Région (CPER) Lorraine 2015-2020, project "Matériaux, Energie, Procédés". 978-1-5386-4291-7/18/$31.00 © 2018 IEEE comfort zones. If we consider a volume of air i (called afterwards zone i ) bounded by n walls, the most determining parameters of atmosphere are the ambient temperature of the air T i , the temperature of the n walls T k i (1 k n), the relative humidity of the air φ i and the flow rate of air [START_REF]ANSI/ASHRAE Standard 55, Thermal Environmental Conditions for Human Occupancy[END_REF]. The objective of this paragraph is to specify the requirements of hygrothermal comfort in current situations (for resting persons).

A. Operative temperature

The operational temperature T o i in zone i is a simplified measure of human thermal comfort derived from air temperature, mean radiant temperature and air speed. When the air speed is less than 0.2 m/s, (as is typical in buildings), T o i can be approached by the arithmetic mean of air temperature in zone i (T i ) and mean radiant temperature T mr i [START_REF]ANSI/ASHRAE Standard 55, Thermal Environmental Conditions for Human Occupancy[END_REF]:

T o i ≈ T mr i + T i 2 (1)
where T mr i is defined by:

T mr i = k S k i × T k i k S k i (2)
A model of adaptive thermal comfort was proposed by [START_REF] Mccartney | Developing an adaptive control algorithm for europe[END_REF] based on objective investigations carried out in buildings in situ. The model consists of a simple linear regression of the comfort temperature T comf (k) depending on the filtered outside air temperature:

T comf (k) = 0.049 T rm (k)+9.2 if T rm (k) ≤ 283.15 0.206 T rm (k)-34.85 if T rm (k) > 283.15 (3)
with T rm (k) such that:

T rm (k) = 0.8T rm (k -1) + 0.2T dm (k -1) (4)
and T dm (k -1) the daily mean outdoor temperature T e on day (k -1).

A first criterion of thermal comfort is to ensure that operational temperature T o i (1) belongs to the comfort temperature interval defined by T comf ± 2 K. Nevertheless, a source of thermal discomfort remains if the difference between the ambient temperature T i and the wall surface temperature T k i is too big. So, the comfort objectives have to include |T i -T k i | < 8 K for glass areas and 5 K for opaque walls.

B. Humidity

Moisture may condense on surface (Σ k ) bounding zone i when the relative humidity of the air φ i is important and when the difference of temperature between T i and that of the walls T k i is too big. φ i depends on the vapor pressure P v i and the saturation vapor pressure P vs i in zone i according to the relation:

φ i (%) = 100 P v i P vs i = 100 ρ v i ρ vs i ( 5
)
where ρ v i and ρ vs i are respectively the densities of water vapor in zone i associated with vapor pressure P v i and P vs i . Saturation vapor pressure P vs i , which does not decrease with the barometric pressure, corresponds to the maximum partial pressure that can reach P v i without generating condensation.

It can be estimated from Tetens' formula [START_REF]ANSI/ASHRAE Standard 55, Thermal Environmental Conditions for Human Occupancy[END_REF]:

P vs i = 610, 78 × exp 17.269(T i -273.15) T i -35.85 (6)
To avoid the surface condensation on wall (Σ k ), it is therefore necessary that the temperature T k i be higher than the dewpoint temperature T d i of this zone i , which is the temperature from which the air becomes saturated with water vapor and that water condenses. For a relative humidity equal to φ i and an internal temperature T i , T d i is such that:

610, 78 × exp 17.269(T d i -273.15) T d i -35.85 = φ i 100 P vs i (7)
In addition, just as the difference in temperature between two faces of a (Σ k ) creates a heat flow, the difference in humidity revealed by the vapor pressure creates a gradient in the material and therefore a "flow" of moisture. To avoid any condensation phenomenon in the heart of the wall, it is imperative that at all points of this wall the vapor pressure is lower than the saturation vapor pressure P vsk . Furthermore, air quality standards recommend for France a humidity rate φ i (%) between 40 and 60 %. So, overall, we will have to ensure in terms of humidity:

     40% < φ i (%) < 60% P vk i < P vs i ⇒ T k i > T d i P vk < P vsk (8)

III. HYGROTHERMAL MODELING FOR BUILDINGS

The different models of building components are developed in this section. The approximations made are justified by taking as an example the city of Nancy (France, ≈ 230 m altitude).

A. Moisture and heat balances into one zone 1) Sensible heat balance: Taking into account the thermal capacity of moist air, i.e that of dry air and water vapor, the energy balance in a i area delimited by n walls consists in equaling the sum of all incoming and outgoing heat fluxes

S k i Area of surface (Σ k ) in contact with zone i m 2 V i Volume of zone i V k Volume of material (Σ k ) m 3 δ pk Vapor permeability of material (Σ k ) kg/(m•s•Pa) ( e k δ pk = R k is the moisture resistance) λ k Thermal conductivity of material (Σ k ) W/(m•K) ( e k λ k = R k is the thermal resistance) ρ da i Density of dry air in zone i kg/m 3
with the temporal variation of enthalpy H i (t) of moist air (See table I):

Ḣ i (t) = (ρ da i (t)C p da + ρ v i (t)C pv ) V i Ṫ i (t) = C p da P da i (t) R da T i (t) 1 + 1.15 P v i (t) P da i (t) V i Ṫ i (t) = j q j i (t) (9)
knowing that ρ da i (t) =

P da i (t) R da T i (t) , ρ v i (t) = P v i (t) R v T i (t) .
This sensible heat balance is a nonlinear differential equation. However, as in general the pressure at sea level is between 980 and 1040 hPa (usual meteorological conditions in France) and that the temperature in a i zone is such that 288 < T i (t) < 308 K, (9) can be approached for the city of Nancy by the following linear differential equation with a maximum error of ±9% on the estimation of the term weighting Ṫ i (t):

1.15 × C p da V i Ṫ i (t) = j q j i (t) (10)
2) Moisture balance: Into zone i , the density of water vapor ρ v i (t) does not vary as the temperature or pressure changes, as long as moisture is not added or removed. It follows:

V i ρv i (t) = j g j i (t) (11)
3) Convective transfer: Since the area i is assumed to be in contact with n walls (Σ k i ), each surface side S k i , of interior temperature T k i (t) and of water vapor density equal to ρ v i (t), is the place of a convective heat / moisture flow with the air of the i zone according to the relations [START_REF] Webb | Standard nomenclature for mass transfer processes[END_REF]:

           g 1 i (t) = - n k=1 P i (t) n a h m S k i (ρ v i (t)-ρ vk i (t)) q 1 i (t) = - n k=1 h c S k i (T i (t)-T k i (t))+h v g 1 i (t) (12)
where h m , h c , h v and n a are defined in table I and

P i (t) n a ≈ 2 500.
Temperatures T k i and densities of water vapor ρ vk i (t) in relation ( 12) can be determined by means of moisture and heat balances in walls surrounding zone i .

B. Moisture and heat balances in walls / windows

Let us consider a homogeneous and isotropic material (Σ k ) of thickness e k delimited by two parallel surfaces (Σ k i ) and (Σ k j ), this environment being able to represent for example a wall, a window, a door, etc. It is assumed that the temperature and the partial pressure of water vapor are identical at every point of each surface (one-dimensional conduction through the material). Let (T k i (t), P vk i (t), ρ vk i (t)) and (T k j (t), P vk j (t), ρ vk j (t)) be the temperature, the partial vapor pressure and the vapor density characterizing each of the two surfaces (Σ k i ) and (Σ k j ) respectively. It follows from this hypothesis that the mass and thermal flows are carried out perpendicularly to these surfaces. The remainder of this section presents the mass and sensible heat balances performed on the equidistant central surface of e k /2 of (Σ k i ) and (Σ k j ). Let T k (t), P vk (t) and ρ vk (t) be the temperature, the vapor pressure and the density of vapor in (Σ k ). The model of the coupled heat and moisture transport is based on the well-known laws of Fourier, Fick and Darcy.

1) Moisture balance: In porous unsaturated materials, moisture transfer occurs in two forms: liquid and water vapor. The transfer of liquid water is induced by a capillary pressure gradient but is not considered in this work (building supposed healthy). The transport of water vapor is a diffusional process under a vapor pressure gradient, governed by Fick's law (

g v S k = -δ pk grad(P v )).
The mass balance at (Σ k ) is given by (See table I):

V k ρvk (t) = -2S k ρ k C mk × l = i , j
δ pk e k (ρ vk (t)-ρ vk l (t)) [START_REF] Boyer | Thermal building simulation and computer generation of nodal models[END_REF] with:

ρ vk (t) = ρ k C mk P vk (t) ρ vk l (t) = ρ k C mk P vk l (t) (14)
2) Heat balance: The transfer of heat in the porous environment occurs under three forms : a flow of conduction under a gradient of temperature expressed by the law of Fourier ( q cond S k = -λ grad(T )), an advection of sensitive heat through the flows of liquid and steam (classically neglected [START_REF] De | Simultaneous transfer of heat and moisture in porous media[END_REF]) and a transfer of latent heat of change of phase driven by the steam ( q v S k = -h v δ pk grad(P v )). The balance equation of the energy conservation is (See table I):

ρ k C pk V k Ṫk (t) = -2S k × l = i , j λ k e k (T k (t)-T k l (t)) - 2S k σh v ρ k C mk × l = i , j δ pk e k (ρ vk (t)-ρ vk l (t)) (15)
Remark 1: For a wall (Σ k ), if the Biot numbers defined for the transfers of heat and transfer of mass respectively by:

Bi = h c e k λ k ; Bi m = h m e k δ pk (16)
are small in front of 1 (Bi < 0.1 and Bi m < 0.1 by considering the coefficients h c and h m for both surfaces (Σ k i ) and (Σ k j )) then the resistance to conduction within (Σ k ) is much lower than the resistance to convection across the air boundary layer. In this case, the assumption of uniform temperature/partial vapor distribution in (Σ k ) is reasonable and Equations ( 15) and ( 13) become

ρ vk (t) = ρ vk i (t) = ρ vk j (t) T k (t) = T k i (t) = T k j (t)
C. Mass and heat exchanges within a zone 1) linked to the air movements: We consider air exchange and denote in m 3 /s the effective flow rate of outdoor air, including infiltration and mechanical ventilation. Assuming that the air is renewed from outside at a temperature T e (t) due to a constant flow Q a , a mass algebraic flow g 2 i (t) / heat q 2 i (t) exists and is defined by:

g 2 i (t) = -Q a (ρ v i (t) -ρ v e (t)) q 2 i (t) = -1.15 × Q a C p da (T i (t) -T e (t)) (17)
Let us assume that a CMV is used to circulate the air between rooms i and j . It is assumed that the air taken in the 'source' room is replaced by the outside air through the openings at the level of the windows. On the other hand, air received by the 'destination' room replaces a part of its air that is equally displaced through the openings at the level of the windows. Two air masses of different psychrometric characteristics are thus mixed. Let (ρ v i , T i ) (t on ) and (ρ v j , T j ) (t on ) be these characteristics within each zone i and j at the moment t on the CMV is switched on. If we take the example that the air goes from zone i to zone j , the humidity / heat balance is:

g 3 j (t) = -Q cmv (ρ v j (t) -ρ v i (t)) q 3 j (t) = -1.15 × Q cmv C p da (T j (t) -T i (t)) ( 18
)
g 3 i (t) = -Q cmv (ρ v i (t) -ρ v e (t)) q 3 i (t) = -1.15 × Q cmv C p da (T i (t) -T e (t)) (19)
with Q cmv the air flow of the CMV in m 3 /s.

2) linked to heating/cooling systems and air conditioning: Air treatment systems in buildings do not just cool/heat but also treat, filter and humidify indoor air. In this paper, we assume that the temperature T h , T c (in K) and the vapor density ρ v,h , ρ v,c of the forced-air (in kg/m 3 ) in heating / cooling mode are known. The focus is on reversible heat pumps (HP). A heat pump is a thermodynamic device that allows capturing heat from outdoor air and transmitting it inside the building. If it is reversible, the HP is capable of providing air conditioning in the summer. Essentially, and without taking into account the transitional aspects, the equation translating this energy flow in zone i is in heating / cooling mode:

g 4 i (t) = -Q HP ρ v i (t) -ρ v,h/c (t) q 4 i (t) = -1.15 × Q HP C p da T i (t) -T h/c ( 20
)
with Q HP the air flow of the HP in m 3 /s.

3) linked to occupancy: Table II specifies the energy and the water vapor released by one person depending on his activity and the temperature of zone i . 4) linked to the absorbed shortwave solar radiations: The sun provides heat that greatly influences the temperature within the platform rooms. To estimate the solar energy transmitted, it is necessary to follow the positioning of the sun and to take into account direct solar and diffuse radiation. In order to understand the calculations to quantify the global radiation, it is important to define not only the various solar components, but also the different angles used for calculations [START_REF] Calow | Blackwell's concise encyclopedia of environmental management[END_REF]. These angles are the local latitude (Φ), the hour angle (ω) in the local solar time, the solar declination (δ s ) and the solar azimuth (φ s ). The energy related to direct solar radiation is such that:

q 6 i (t) = 1260 S wk exp -1 2.3 (α s (t)+3) sin(α s (t)) × cos (β(t)) (21)
where S wk is the surface of the windows of (Σ k ), α s (t) is the solar elevation angle and β(t) is the angle of incidence between the suns rays and the face of windows. As for the diffused radiation energy, it is calculated by the relationship:

q 7 i (t) = S wk 62.5 ×(1 + cos(i)) sin α s (t) 0.4 +108 ×(1 -cos(i)) sin α s (t) 1.22 (22)
IV. HYGROTHERMAL MODELING BY NODAL METHOD -APPLICATION TO THE "ECO-SAFE" PLATFORM By this type of approach, each node represents volumes of material or components, while inter-node connections represent the flow paths of thermal and mass fluxes (walls, doors, windows, etc.) [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF], [START_REF] Boyer | Thermal building simulation and computer generation of nodal models[END_REF]. The graphic representation of the model has the same appearance as an electrical diagram, hence the analogue method name where:

• the temperature and density of water vapour are represented by electrical potentials; • energy flows and mass transfers are represented by current intensities; • thermal and moisture resistances are represented by electrical resistances; • the thermal and moisture retention capacities are represented by capacitors.

For example illustration, the "Eco-Safe" platform shown in Figure 1 is considered. The electrical circuits translating the energy and mass models within this platform are represented at the top of Figure 2. For each zone, a global model is defined by assembling the different hygrothermal blocks. It consists of four walls, a double glazed window, a floor and a ceiling. The model of this part is illustrated for room 1 at the bottom right corner of Figure 2. The consideration of heterogeneous walls (e.g. window and wall) result in as many circuits RC paralleled. An example of a wall model considering a window with double glazing is represented at the bottom left corner of Figure 2. For the sake of simplicity and due to the low impact of term h v g 1 i (t) in ( 12) during winter conditions, the latter has been neglected. Parameters R and C of the simulator were determined from the nature and characteristics of the different materials used in the wall dividers, walls and windows of the platform. They have been initialized taking into account theoretical values of certain parameters specified in table III and according to the following analogies:

• for zone i :

       R T,k i = 1 h c S k i R H,k i = 1 2 500 h m S k i ; C T, i = 1.15 C p da V i C H, i = V i ( 23
)
• for wall / windows / ceil / floor k:

     R T,k = e k 2λ k R H,k = e k 2δ k ; C T,k = ρ k C pk e k C H,k = ρ k C mk e k (24)
They have been adjusted from the temperature curves actually observed on the system. According to figures 3 and 4, it appears that the profiles of actual and estimated temperature values are found to be similar. The simulated curves were generated from the entries T e (t) and ρ v e (t) (Figure 5), which are the data issued from the weather station over the period 3/5/2018 -3/9/2018.

V. CONCLUSION

In this paper, a complete methodology has been established for modeling the hygrothermal behavior of a multi-zone office building. The application target is the "Eco-Safe" platform (CRAN Nancy, France). This study takes into account not only all the walls/windows and equipments (CMV, heat pump) fitted to each of the zones but also the effects of solar radiations and occupants. In a practical context, in order to facilitate the determination of a hygrothermal model of a multi-zone building, the proposed nodal method is based on the use of the thermal electric analogy. In perspectives, it would be interesting to implement an efficient control law based on this model.

Fig. 1 .

 1 Fig. 1. The "Eco-Safe" platform.

Fig. 2 .

 2 Fig. 2. Hygrothermal modeling of "Eco-Safe" platform.

Fig. 4 .

 4 Fig. 4. Measured (solid line) and estimated (dashed line) humidity rates in zone 1 .

Fig. 5 .

 5 Fig. 5. Outside temperature (solid line) and humidity rate (dashed line) of outdoor air (3/5/2018 -3/9/2018).

TABLE III THEORETICAL

 III VALUES OF HYGROTHERMAL PARAMETERS[START_REF]Hygrothermal Properties. Tabulated Design Values[END_REF] Fig.3. Measured (solid line) and estimated (dashed line) temperatures for zoom1 . An estimation of the operative temperature is represented in dotted line.

			ρ k	δ pk (e -12 )	C pk	λ k
	Concrete	2 400	1.4	1 000	2
	Air		1.23	185	1 006 0.025
	Polystyrene	40	3.1	1 450	0.04
	Glass		2 580	0	720	1.2
	26
	25
	24
	23
	22
	21
	20
	19
	18
	Mar 05	Mar 06	Mar 07	Mar 08	Mar 09
						2018
	48
	46
	44
	42
	40
	38
	36
	34
	32
	30
	28
	Mar 05	Mar 06	Mar 07	Mar 08	Mar 09
						2018


Farah Gabsi
email: farah.gabsi@univ-lorraine.fr

Frédéric Hamelin
email: frederic.hamelin@univ-lorraine.fr

Nathalie Sauer
email: nathalie.sauer@univ-lorraine.fr

Hygrothermal modeling

Keywords: buildings, hygrothermal modeling, modelpredictive control, energy efficiency, thermal comfort

published or not. The documents may come

I. INTRODUCTION

Since the entry into force of the Kyoto Protocol on 2005, the building sector has been encouraged to reduce its consumption of energy and materials. These efforts have resulted in thermal regulations, environmental quality or energy performance labels. More recently, in France, the Grenelle Environment has confirmed the potential of the building sector in terms of energy saving.

The envelope of a building is the interface between the inner ambience and the external environment. It is the seat of transfers of heat, humidity, and air, which determine the internal climate of the building. Most often, it must be supplemented by a regulation system (heating, humidification, refreshment, air renewal...). The introduction of the Model Predictive Control (MPC) in the sector of building turns out very promising to maintain the comfort of the occupants and to generate energy savings. This type of control is favored because it has the advantage of taking into account the disturbance predictions that greatly influence the interior comfort of buildings (outside temperature, sunshine, occupations,...) while anticipating their effects [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF], [2], [START_REF] Salazar | Predictive control of a renewable energy microgrid with operational cost optimization[END_REF], [START_REF] Walker | Application of distributed model predictive approaches to temperature and co2 concentration control in buildings[END_REF]. Recent studies have highlighted the risks associated with humidity in buildings: deterioration of conditions of comfort or hygiene, decrease of energy performance. However, very little work considers the hygrothermal control of buildings [START_REF] Ferroukhi | Effect of coupled heat, air and moisture transfers modeling in the wall on the hygrothermal behavior of buildings[END_REF].

This work has financial support from the Contrat de Plan Etat-Région (CPER) Lorraine 2015-2020, project "Matériaux, Energie, Procédés".

In this paper, a complete methodology of modeling/ controlling a multi-zone office building is proposed, the objective is to improve the hygrothermal comfort of the occupants while minimizing the energy consumed. The modeling step, which is essential for the design of a MPC, is based on a physical description of the building [START_REF] Gabsi | Building hygrothermal modeling by nodal method[END_REF]. Knowing that each control action will be associated with a different descriptor state-space model of the system, the predictive control strategy will be of the hybrid type. In view of the control objectives, the optimization criterion associated to the MPC will integrate the energy cost of each actuator in euros, and will be solved by an iterative algorithm in order to avoid a risk of combinatory explosion. The first part of this work is focused on the hygrothermal modeling of a set of zones. Then a MPCbased control strategy is proposed to ensure the hygrothermal comfort of a set of zones with a minimal energy consumption.

II. HYGROTHERMAL COMFORT

The notion of thermal comfort in a building is related not only to the quality of the interior atmospheres, but also to the amount of energy to be supplied by the equipments. In ISO 7730 [START_REF] Nicol | Adaptive thermal comfort standards in the hot-humid tropics[END_REF], a method was presented to predict the thermal sensation and the degree of discomfort (thermal dissatisfaction) of people exposed to moderate thermal environments and to specify thermal conditions acceptable for comfort. Considering a volume of air i (called afterwards zone i ) bounded by n walls, this standard defines a felt comfort temperature, also known as "operative temperature" or "dry resulting temperature" by:

T o i ≈ T mr i + T i 2 (1)
where T i is the ambient temperature of the air, T k i (1 k n) are the temperatures of the n walls with surfaces equal to S k i and T mr i is the mean radiant temperature [START_REF]ANSI/ASHRAE Standard 55, Thermal Environmental Conditions for Human Occupancy[END_REF] defined by:

T mr i = k S k i × T k i k S k i (2)
A first criterion to consider for the design of the control law is to ensure that operational temperature T o i (1) belongs to a comfort temperature range defined by T comf ± 2 K where T comf (k) is a comfort temperature [START_REF] Mccartney | Developing an adaptive control algorithm for europe[END_REF] related to the outside air temperature. Nevertheless, a source of thermal discomfort remains if the difference between the ambient temperature T i and the surface temperature T k i of each wall is too big. So, the MPC objectives have to include |T i -T k i | < 8 K for glazed surfaces and 5 K for opaque walls.

The last appellation of thermal comfort used by the HQE (High Quality Environmental) standard is the "hygrothermal comfort", because the latter depends not only on the temperature, but also on the air humidity [START_REF] Ebert | Green Building Certification Systems: Assessing sustainability -International system comparison -Economic impact of certifications[END_REF]. More precisely, the hygrothermal comfort depends on the air and surface temperatures, humidity and air movements. To measure the moisture content of the air, one defines the relative humidity of the air φ i according to the relation:

φ i (%) = 100 P v i P vs i = 100 ρ v i ρ vs i ( 3
)
where ρ v i and ρ vs i are the densities of water vapor in zone i associated respectively with vapor pressure P v i and saturation vapor pressure P vs i that is such that [START_REF]ANSI/ASHRAE Standard 55, Thermal Environmental Conditions for Human Occupancy[END_REF]:

P vs i = 610, 78 × exp 17.269(T i -273.15) T i -35.85 (4)
According to the air quality standards recommended for France, the MPC will have to ensure a humidity rate φ i (%) of between 40 and 60 %.

III. HYGROTHERMAL MODELING

In this paragraph is presented the set of balance equations governing the dynamic behaviour of temperature and relative

S k i Area of surface (Σ k ) in contact with zone i m 2 V i Volume of zone i V k Volume of material (Σ k ) m 3 δ pk Vapor permeability of material (Σ k ) kg/(m•s•Pa) ( e k δ pk = R k is the moisture resistance) λ k Thermal conductivity of material (Σ k ) W/(m•K) ( e k λ k = R k is the thermal resistance) ρ da i Density of dry air in zone i kg/m 3
humidity within a zone i . The approximations made are justified by taking as an example the city of Nancy (France, ≈ 230 m altitude). The definition of the different parameters is specified either in the text or in the Table I.

A. Energy and Moisture Balances

The energy balance equation in zone i delimited by n walls matches the sum of all incoming or outgoing heat fluxes to the temporal variation of enthalpy H i (t) of moist air:

Ḣ i (t) = (ρ da i (t)C p da + ρ v i (t)C pv ) V i Ṫ i (t) = C p da ρ da i (t) 1 + C pv ρ v i (t) C p da ρ da i (t) V i Ṫ i (t) = C p da P da i (t) R da T i (t) 1 + 1.15 P v i (t) P da i (t) V i Ṫ i (t) = j q j i (t) (5)
knowing that

ρ da i (t) = P da i (t) R da T i (t) and ρ v i (t) = P v i (t) R v T i (t) .
This nonlinear differential equation can be approximated by a linear differential equation by considering local weather conditions. Indeed, according to (4) and inequality 288 K < T i (t) < 308 K, it follows knowing that generally in France 980 hPa < P 0 i (t) < 1040 hPa:

1.11 -1.25e -4 z < P da i (t) R da T i (t) 1+1.15 P v i (t) P da i (t) < 1.26 -1.42e -4 z (6)
taking into account atmospheric pressure variations P z i (t) as a function of altitude (z in metres):

P z i (t) ≈ P 0 i (t) 1 - 0.0065z 288.15 5.255 (7)
Therefore, for the city of Nancy, (5) can be approached by the following linear differential equation with a maximum error of ±9% on the estimation of the term weighting Ṫ i (t):

1.15 × C p da V i Ṫ i (t) = j q j i (t) (8)
As for the moisture balance equation in zone i , it is given by:

V i ρv i (t) = j g j i (t) (9)
If the surfaces (walls/windows/doors/ceiling/floor) surrounding zone i are now considered, the energy balance equation associated with each homogeneous and isotropic material (Σ k ) (1 k n) of thickness e k delimited by the two parallel surfaces (Σ k i ) and (Σ k j ) is given by:

ρ k C pk V k Ṫk (t) = -2S k × l = i , j λ k e k (T k (t)-T k l (t)) - 2S k σh v ρ k C mk × l = i , j δ pk e k (ρ vk (t)-ρ vk l (t)) (10)
with:

ρ vk (t) = ρ k C mk P vk (t) ρ vk l (t) = ρ k C mk P vk l (t) (11)
and where the temperature, the partial vapor pressure and the vapor density characterizing the middle of (Σ k ) and each of the two surfaces (Σ k i ) and (Σ k j ) are given by (T k (t), P vk (t), ρ vk (t)), (T k i (t), P vk i (t), ρ vk i (t)) and (T k j (t), P vk j (t), ρ vk j (t)) respectively.

As for the mass balance equation at (Σ k ), it is given by:

V k ρvk (t) = -2S k ρ k C mk × l = i , j δ pk e k (ρ vk (t)-ρ vk l (t)) (12)
The link between equations ( 8), ( 9) and ( 10), ( 12) is established by considering the convective heat / moisture flux between air of zone i and the n surfaces [START_REF] Webb | Standard nomenclature for mass transfer processes[END_REF]:

           g 1 i (t) = - n k=1 P i (t) n a h m S k i (ρ v i (t)-ρ vk i (t)) q 1 i (t) = - n k=1 h c S k i (T i (t)-T k i (t))+h v g 1 i (t) (13)
B. Main Mass and Heat Exchanges within a Zone 1) Heating/cooling systems: Air treatment systems in buildings do not just cool/heat but also treat, filter and humidify indoor air. In this paper, we assume that the temperature T h , T c (in K) and the vapor density ρ v,h , ρ v,c of the forced-air (in kg/m 3 ) in heating / cooling mode are known. The focus is on reversible heat pumps (HP). These are thermodynamic devices for transferring heat from a "cold source" to a "hot source", although naturally the heat is diffused from the warmest to the coldest place. Thus, the thermal comfort in a building can be ensured by heat pumps that refresh in summer and heat in winter. Essentially, and without taking into account the transitional aspects, the equations translating this energy flow in zone i is in heating/cooling mode:

g 2 i (t) = -Q HP ρ v i (t) -ρ v,h/c (t) q 2 i (t) = -1.15 × Q HP C p da T i (t) -T h/c (14)
with Q HP the air flow of the HP (in m 3 /s).

2) Occupancy and equipments: Table II shows the energy and the water vapour released by one person depending on his activity and the temperature of the area in which he is located. According to [START_REF] Calow | Blackwell's concise encyclopedia of environmental management[END_REF], q 3 i , presented within Table II, represents an estimate of the total heat emitted by a person. It includes body heat (related to body temperature 310 K) and the latent heat by the production of water vapour in breath and perspiration.

q 3 i (t) ≈ -γ(T i (t) -310)+h v g 3 i (t) (15)
Many other equipments produce heat in a room. Being of a diverse nature and of a different power (lighting, computer equipment, furnace,...), they will not be considered in this work.

T i g 3 i q 3 i g 3 i q 3 i g 3 i q 3 i (C) (g/h) (W) (g/h) (W) (g/
3) Solar radiations: The total solar radiation received by a surface, called incident solar irradiation (or even global energy irradiance), is defined as the sum of three components:

• Direct solar radiation. This component cancels if the sun is hidden by clouds or by an obstacle; • Diffused radiation, corresponding to the radiation received from the celestial vault, without direct radiation. This energy diffused by the atmosphere and directed towards the surface of the Earth, can reach 50% of the overall radiation received, when the sun is low on the horizon, and 100% for a fully covered sky; • Reflected radiation, corresponding to the radiation reflected by the external environment, especially the soil, whose reflection coefficient is called "albedo".

To quantify the global radiation, it is necessary to consider the local latitude (Φ), the hour angle (ω) in the local solar time, the solar declination (δ s ) and the solar azimuth (φ s ) [START_REF] Calow | Blackwell's concise encyclopedia of environmental management[END_REF]. The energy related to direct solar radiation is then such that:

q 4 i (t) = 1260 S wk exp -1 2.3 (α s (t)+3) sin(α s (t)) × cos (β(t)) (16)
where S wk is the glazed surface of (Σ k ), α s (t) is the solar elevation angle defined by:

α s (t) = arcsin (sin(Φ) sin (δ s (t)) + cos(Φ) cos (δ s (t)) cos (ω(t))) (17)
and β(t) is the angle of incidence between the suns rays and face of the windows.

β(t) = arccos (cos (α s (t)) sin(i) cos (φ s (t) -or) + sin (α s (t)) cos(i)) (18)
The determination of α s (t) and β(t) requires the calculation of the solar angle ω(t) (rad):

ω(t) = π 12 (12 -h(t)) (19)
of the solar azimuth angle φ s (t):

φ s (t) = arcsin cos (δ s (t)) sin (ω(t)) cos (α s (t)) (20)
and of the current sun declination δ s (t), which is the angle between the sun direction and the plane of the Equator : δ s (t) = 0.18 π (6.92-400 cos (Γ(t))-6.8 cos (2Γ(t))

-2.7 cos (3Γ(t))+70 sin (Γ(t)) +0.9 sin (2Γ(t))+1.5 sin (3Γ(t)))

where:

• h(t) is the true solar time;

• Γ(t) = 2π (d n (t) -1) /365 corresponds to the day angle calculated from day number d n (t) ∈ R taking into account hours and minutes; • latitude Φ is determined according to the position of the building (Φ = 0.85 rad for Nancy, France); • the orientation of the windows is indicated by variable or; • i corresponds to the inclination of the windows according to the horizontal (i = π/2 rad).

As for the energy associated to the diffused radiation energy, it is calculated by the relationship:

q 5 i (t) = S wk 62.5 ×(1 + cos(i)) sin α s (t) 0.4 +108 ×(1 -cos(i)) sin α s (t) 1.22 (22)
with an albedo approximatively equal to 20 %.

Remark 1: When the sun's rays reach an opaque wall, part of the solar radiation is absorbed by the various materials that make up this wall. There is no radiation transmitted since the wall is opaque. The inertia of the wall determines the response time. In order to reflect the solar energy transmitted by these opaque walls and to consider the inertia of the latter, it is possible to use the concept of equivalent outdoor temperature, which integrates the solar radiation and the characteristics of the opaque walls (absorption capacity and inertia).

C. State-Space Modeling

According to the previous developments, the dynamical hygrothermal behavior of zone i bounded by n walls (Σ k ) can be represented by a set of N s descriptor state-space models, each of them being associated with a particular control scenario (HP (ON-OFF),...):

                   E i Ẋ i (t) = A i X i (t) + l = j , e f l T l (t) ρ v l (t) + n k=1 B k i T k (t) ρ vk (t) +f HP T h/c ρ v,h/c + j≥5 q j i (t) g j i (t) T i (t) ρ v i (t) = C i X i (t) (23) with X i (t) =            T i (t) ρv i (t) T1 i (t) . . . Tn i (t) ρv1 i (t) . . . ρvn i (t)            , A i =            * * * • • • * * • • • * 0 * 0 • • • 0 * • • • *       * . . . *             0 . . . 0             0 . . . 0             * . . . *                         * 0 • • • • • • • • • 0 0 * 0 • • • • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 • • • • • • 0 * 0 0 • • • • • • • • • 0 *                              , B k i = 0 (k+1)×2 * * 0 (2n-k)×2 , C i = [ I 2 0 2×2n ], E i = * 0 0 * 0 2×2n 0 2n×2 0 2n×2n
.

The connection with the other zones j of the building is carried out by the state-space model of all the walls (Σ k ) surrounding zone i :

Ẋk (t) = A k X k (t) + B k     T k i (t) T k j (t) ρ vk i (t) ρ vk j (t)     (24) with X k (t) = T k (t) ρ vk (t) , A k = * * 0 * , B k = * * * * 0 0 * * .
In the particular case where wall (Σ k ) is in contact with the outside air ( j = e ) characterized by a temperature T e (t) and a density of water equal to ρ v e (t), (24) becomes:

Ẋk (t) = A k X k (t) + B k     T k i (t) T e (t) ρ vk i (t) ρ v e (t)     (25)
knowing, according to equations ( 8), ( 9), ( 10), ( 12) and ( 13), that

T k e (t) = a k T k (t) + b k T e (t)

IV. MODEL PREDICTIVE CONTROL

The principle of the predictive control [2] is to optimize a cost function allowing to describe the aim of control over a finite time horizon. At any moment, an optimal control sequence is calculated to minimize the cost function on a prediction horizon. Only the first element is applied to the system.

A. Cost Function and Constraints

Since the objective is to ensure the thermal comfort with a minimum of energy consumption, the cost function must reflect these performances in a mathematical formulation. In the light of sections II and III, the cost function is chosen such that:

J = min i1...,i N Np j=1 l α l ,j T o l ,ij (k+j)-T comf (k+j) 2 + β l ,j φ l ,ij (k+j)-50 2 + e ij (k+j) (26)
with :

• N p : the prediction horizon;

• T o l ,ij (k + j): the predicted operative temperature T o l (1) in zone l at time k+j according to the sequence of control scenario i 1 , . . . , i j ; • φ l ,ij (k +j): the predicted relative humidity φ l (%, [START_REF] Salazar | Predictive control of a renewable energy microgrid with operational cost optimization[END_REF] in zone l at time k + j according to the sequence of control scenario i 1 , . . . , i j . It is determined by means of the prediction of ρ v l (k + j), equations (3-4) and with relation ρ vs l (k+j) = P vs l (k+j) R v T l (k+j) ;

• T comf (k+j): the comfort temperature at time k+j;

• e ij (k+j): the energy cost (in e) of the control scenario i j . It depends on the electricity pricing C e (k + j) (in e/kWh), which varies significantly between off-peak and peak hours. Indeed, national electricity consumption varies by time of day depending on the activities of businesses and individuals. The means of production required to respond instantly to this demand are activated according to a "merit order" logic: the higher the power consumption, the more electricity is expensive to produce.

under the constraints:

• of N s hygrothermal models (23-25). At each step of prediction k + j, the model to be adopted depends on the control scenario i j selected; • of comfort (according to section II):

               T comf (k+j)-2 ≤ T o l ,ij (k+j) ≤ T comf (k+j)+2 T l ,ij (k+j)-T k l ,ij (k+j) < 5 (or 8) K 40% < φ l ,ij (k+j) < 60% ρ vk l ,ij (k+j) < ρ vs l ρ vk,ij (k+j) < ρ vsk
These constraints can be easily integrated in the optimization problem (26) because all variables T l , T k l , ρ vk l , ρ v l and ρ vk are components of state vectors X l and X k (23-25).

Remark 2: In the context of minimizing the energy consumption, it could be interesting to specify that the hygrothermal constraints of zone l are not the same in the absence or in the presence of people. This notion can be reflected using time-varying values for weights α l ,j and β l ,j (for example by taking α l ,j = β l ,j = 0 for an unoccupied zone at time k+j).

Remark 3:

The disturbance inputs T e (k +j), ρ v e (k +j), q 4 l (k + j) and q 5 l (k + j) (16-22) are estimated through short-term weather forecasts. q 4 l (k + j) and q 5 l (k + j) require a basic calculation according to the orientation of the windows and predicted solar radiation. q 3 l (k + j) and g 3 l (k + j) (section III-B2) can be predicted based on the occupancy planning in zone l .

B. Iterative Approach to Solving Optimization Problem

The minimization of criterion J (26) for important values of N p and/or N s quickly becomes impossible using Matlab routines due to the combinatorial explosion in computation time. In practical terms, N p and/or N s can not exceed a few units. Yet it can be very interesting to consider a large value for prediction horizon N p . Indeed, an optimal control scenario on a low prediction horizon may become inappropriate on a larger horizon. This happens especially when the heating and/or cooling capabilities of the building are modest or when weights α l ,j and β l ,j are time-dependent (See Remark 2). In order to consider larger value for the prediction horizon N p , we have developed an iterative approach [START_REF] Gabsi | Energy efficiency of a multizone office building: MPC-based control and simscape modelling[END_REF], whose purpose is to estimate the optimal control scenario, solution of (26), with a controlled computational load. The idea consists in conserving at time of prediction k + j only a limited number (l) of scenarios from among all those possible.

An illustrative example is given in Fig. 1 with N p = 5, N s = 4 and l = 10. In this three-dimensional space, (x, y, z) represents j, φ l ,ij (k+j), T o l ,ij (k+j) . The upper view shows all the N Np s paths associated to all the combinations of i j , while the lower view shows some paths judiciously selected according to their energy costs. Each arc extending from one vertex (x = j -1) to another (x = j) is associated with a control scenario i j and its associated energy cost e ij . Figures 2 more accurately represent vertices φ l ,ij (k+j), T o l ,ij (k+j) (for j = 1, ..., 5) belonging to the cyan transverse planes in Fig. 1.

At each time of prediction k+j, blue points of Fig. 2 denote the end-vertices of all the paths generated in a comprehensive manner (related to Fig. 1a). On the other hand, points represented by red circles and red discs in Fig. 2 correspond to the end-vertices of arcs (related to Fig. 1b), which have been generated at the previous time of prediction k+j-1 from a set of l "optimal" points φ l ,ij-1 (k+j -1), To l ,ij-1 (k+j -1) . The latter, represented by red discs, are selected owing to their representativeness of the points circled in red (and so of the blue points) and because of the low cost of paths (φ l (k), T o l (k)) → φ l ,ij-1 (k+j -1), To l ,ij-1 (k+j -1)

according to criterion J.

These l "optimal" points (red discs) are determined at each time of prediction using principal component analysis [START_REF] Jolliffe | Principal Component Analysis, ser. Springer Series in Statistics[END_REF]. The two principal components of the data set containing the coordinates φ l ,ij (k+j), To l ,ij (k+j) of all the points circled in red define two orthogonal directions (V 1 , V 2 ), which appear in thick line in Fig. 2. The projection of all the points represented by red circles on V 1 and V 2 makes it possible to define a rectangle encompassing all these points. This rectangle is then cut into l small rectangles in which l "optimal" points (represented by red discs) φ l ,ij-1 (k+j -1), To l ,ij-1 (k+j -1) are conserved for the next time of prediction.

It has to be mentioned that more the number l is large, more the approximation of the optimal control scenario is good with this iterative approach. On the other hand, the computational load is greatly reduced with proposed iterative approach. Indeed, if the number of arcs is at the maximum equal to l•N s •N p for this algorithm , it can reach (N Np+1 s -N s )/(N s -1) by an exhaustive search for all possible control scenarios.

V. CONCLUSION

In this paper, a complete methodology has been established for modeling a multi-zone office building and controlling its hygrothermal comfort with the objective of minimizing energy costs. Based on a set of descriptor state-space models and on an energy cost function, an original model predictive control has been designed for considering large prediction horizon without any risk of combinatory explosion. Future research will consist in validating this methodology on a laboratory platform.

  and ρ vk e (t) = c k ρ vk (t) + d k ρ v e (t) (coefficients a k , b k , c k and d k being related to h c , h m , λ k and δ pk ).

Fig. 1 .

 1 Fig. 1. Prediction of φ l ,i j (k+j), T o l ,i j (k+j) for j = 1, ..., Np using a) all the combinations i j ∀j, b) the proposed iterative approach.

TABLE I

 I

		:
		NOMENCLATURE
	Par.	Description	SI units
	C mk Specific moisture capacity of material (Σ k )	kg/(kg•Pa)
	Cp da Specific heat capacity of dry air (= 1 006)	J/(kg•K)
	C pk	Specific heat capacity of material (Σ k )	J/(kg•K)
	Cp v	Specific heat capacity of water vapor	J/(kg•K)
		(= 1 860)
	hc	Heat transfer coefficient (≈ 8.3 (interior)	W/(m 2 •K)
		and ≈ 34 (exterior; normal conditions)
		[11])
	hm	Water vapor transfer coefficient	kg/(m 2 s•Pa)
		(≈ 7e -9 × hc)
	hv	Latent heat of vaporization of water	J/kg
		(= 2.5e 6 )
	na	Number of moles of air per m 3 (≈ 41)	mol/m 3
	P da i Dry air pressure of zone i	Pa
	R da	Gas constant for dry air = 287.1	J/(kg•K)
	Rv	Gas constant for water vapor = 461.5	J/(kg•K)

TABLE II

 II

		:
	HEAT/MOISTURE FLUX
	Rest to	Light to sustai-	Heavy
	light activity	ned activity	physical activity


Yuri Goegebeur

Armelle Guillou

Jing Qin

Robust estimation of the Pickands dependence function under random right censoring

Keywords: Pickands dependence function, censoring, Kaplan-Meier integral

We consider robust nonparametric estimation of the Pickands dependence function under random right censoring. The estimator is obtained by applying the minimum density power divergence criterion to properly transformed bivariate observations. The asymptotic properties are investigated by making use of results for Kaplan-Meier integrals. We investigate the finite sample properties of the proposed estimator with a simulation experiment and illustrate its practical applicability on a dataset of insurance indemnity losses.

Introduction

Multivariate extreme value statistics deals with the estimation of the tail of a multivariate distribution function based on a random sample. When studying multivariate extremes, a natural question is how to quantify extreme dependence between two or more random variables. Usually, the copula function is used as a margin-free description of the dependence structure between several random variables. Indeed, according to Sklar's theorem [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 marginal distribution functions F Y p1q and F Y p2q 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 Y p1q py 1 q, F Y p2q py 2 qq .

(

) 1
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 Y ˆlogpy 2 q logpy 1 y 2 q ˙˙,

where A Y : r0, 1s Ñ r1{2, 1s is the Pickands dependence function, which is convex and satisfies maxtt, 1 ´tu ď A Y ptq ď 1, see [START_REF] Pickands | Multivariate extreme value distributions[END_REF]. Throughout the paper we assume that pY p1q , Y p2q q follows a joint distribution with an underlying extreme value copula.

Since a copula function allows to model efficiently the dependence between several random variables, it becomes more and more popular in financial or actuarial applications. To illustrate our methodology, we consider the insurance company loss and expense application by [START_REF] Frees | Understanding relationships using copulas[END_REF]. The dataset, included in the R package copula, comprises 1500 pairs containing information on general liability claims, the first component being the indemnity payment (loss) and the second one an allocated loss adjustment expense (ALAE). The latter is related to the settlement of individual claims, e.g., expenses for lawyers or claim investigation. A crucial question is the possible dependence between the two components, loss and ALAE, which has to be accounted for if we are interested in actuarial applications, such as, e.g., pricing an excessof-loss reinsurance treaty when the reinsurer shares the claims settlement costs. For instance [START_REF] Micocci | Loss-ALAE modelling through a copula dependence structure[END_REF] (see also [START_REF] Cebrián | Analysis of bivariate tail dependence using extreme value copulas: an application to the SOA medical large claims database[END_REF] motivate the use of copula functions in that context with the aim of building a reinsurance strategy in presence of policy limits and insurer's retentions. As outlined in these contributions, a substantial mispricing can result from the usual independence assumption where the joint distribution is assumed to be the product of the marginals, and thus using copulas is the correct way to model dependence and as such to avoid the undervaluation of the reinsurance premium. However, the estimation of the joint distribution for the losses and expenses is complicated due to the presence of censoring. More specifically, for each claim there is a policy limit, and hence the losses cannot exceed this limit.

In the loss-ALAE dataset, 34 observations have censored losses and these censored observations cannot be ignored since for instance the mean loss of censored claims is much higher than the corresponding mean for the uncensored claims (217 491 against 37 110, see Table 4 in [START_REF] Frees | Understanding relationships using copulas[END_REF]. The scatterplot of the data is given in Figure 1, where the uncensored observations are in grey and the censored ones in black. Overall the scatterplot indicates a reasonably strong relationship between the two variables, but the picture is somehow obscured by censoring in the largest observations and also by potential outliers.

We have thus two issues in the dataset under consideration: the presence of censoring and potential outliers.

Concerning the first issue, we will explore in the present paper nonparametric estimation of the Pickands dependence function when there is random right censoring in the marginal distributions. More precisely, we consider the situation where pY p1q , Y p2q q is right censored by pC p1q , C p2q q, also following a bivariate distribution with an extreme value copula, but now with Pickands dependence function A C . Thus, we observe pminpY p1q , C p1q q, minpY p2q , C p2q q, δ p1q , δ p2q q, where δ pjq :" 1l tY pjq ďC pjq u , j " 1, 2, with 1l E the indicator function on the event E, and interest is in estimating A Y .

Random right censoring has been studied to some extent in the univariate extreme value literature, see, e.g. In q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 500000 1000000 1500000 2000000 0e+00 1e+05 2e+05 3e+05 4e+05 5e+05 Loss ALAE q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1e+01 1e+03 1e+05 1e+07 1e+01 1e+03 1e+05 1e+07 Loss ALAE q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 1: LOSS-ALAE data: scatterplot on the original scale (left) and on a logarithmic scale (right).

multivariate extreme value statistics a central topic is the modelling and estimation of extreme dependence between two or more random variables. Similarly to classical statistics, extreme dependence can be summarised by properly chosen dependence coefficients, like the coefficient of tail dependence, see, e.g., [START_REF] Ledford | Modelling dependence within joint tail regions[END_REF] Concerning the second issue, robust methods must be proposed to prevent possible isolated outliers from completely disturbing the estimate of the joint distribution. To reach this goal, we will propose a robust estimator of the Pickands dependence function of pY p1q , Y p2q q based on the density power divergence method introduced by Basu et al. (1998). In particular, the density power divergence between two density functions f and h is defined as follows

∆ α pf, hq :" $ ' ' & ' ' % ż R " h 1`α pyq ´ˆ1 `1 α ˙hα pyqf pyq `1 α f 1`α pyq  dy, α ą 0, ż R log f pyq hpyq f pyqdy, α " 0.
Here the density function h is assumed to depend on a parameter vector θ and if Y 1 , . . . , Y n is a sample of independent and identically distributed (i.i.d.) random variables according to the density function f , then the minimum density power divergence estimator (MDPDE) of θ is the point p θ minimizing the estimated version (up to a constant independent of θ)

p ∆ α pθq :" $ ' ' & ' ' % ż R h 1`α pyqdy ´ˆ1 `1 α ˙żR h α pyqd p F Y pyq, α ą 0, ´żR log hpyqd p F Y pyq, α " 0,
where p F Y is a suitable estimator of the distribution of Y 1 . In the case of no censoring, p F Y is typically the empirical distribution function p F E Y pyq :" p1{nq ř n i"1 1l tY i ďyu , whereas in the censoring framework where Y i is censored by an independent random variable C i , the famous Kaplan-Meier product-limit estimator (see [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF]

defined as p F KM Y pyq :" 1 ´n ź i"1 " 1 ´δri,ns n ´i `1  1l tZ i,n ďyu , (3)
can be used. Here Z i :" minpY i , C i q, i " 1, . . . , n, Z i,n denotes the i´th order statistic of tZ 1 , . . . , Z n u and δ ri,ns is the concomitant order statistic with respect to Z i,n , i.e., δ ri,ns "

δ k if Z i,n " Z k , i " 1, . . . , n.
The MDPD criterion depends on a parameter α which allows to make a trade-off between efficiency and robustness of the resulting estimator. Indeed, we can observe that for α " 0 one recovers the log-likelihood function, up to the sign, which leads to an efficient but not robust estimator. By increasing α we increase the robustness of the estimator, but decrease its efficiency.

The remainder of the paper is organised as follows. In the next section we introduce the nonparametric MDPDE of the Pickands dependence function A Y under random right censoring. The asymptotic properties of this estimator, consistency and finite dimensional weak convergence, are investigated in Section 3, where we use the asymptotic properties of Kaplan-Meier integrals.

We illustrate the finite sample performance of the estimator with a simulation experiment in Section 4 and in Section 5 we apply the method to the dataset of insurance indemnity losses. The proofs of our results are given in Section 6.

Construction of the estimator

Throughout the paper, for any random variable W , we denote by F W its distribution function.

For convenience we reformulate the model as stated in (1) and (2) into standard exponential margins. Assume F Y pjq , j " 1, 2, are continuous. After applying the transformations r Y pjq " ´log F Y pjq pY pjq q, j " 1, 2, we obtain the following bivariate survival function

G Y py 1 , y 2 q :" P ´r Y p1q ą y 1 , r Y p2q ą y 2 ¯" exp ˆ´py 1 `y2 qA Y ˆy2 y 1 `y2 ˙˙,
for all y 1 , y 2 ą 0. A similar assumption is made for the distribution G C of the random vector p r C p1q , r C p2q q. Let t P r0, 1s. Considering the univariate random variable

r Y t :" min ˜r Y p1q 1 ´t , r Y p2q t ¸,
it is clear that P ´r Y t ą z ¯" e ´AY ptqz , @z ą 0.

Consequently, the distribution of r Y t is an exponential distribution with parameter A Y ptq. Similarly, by defining r

C t :" min ´r C p1q 1´t , r C p2q
t ¯, the random variable r C t follows an exponential distribution with parameter A C ptq. Now, remarking that

r Z t :" minp r Y t , r C t q " min ˜minp r Y p1q , r C p1q q 1 ´t , minp r Y p2q , r C p2q q t ¸, r δ t :" 1l t r Ytď r Ctu " δ t r δ p1q `p1 ´δt q r δ p2q ,
where δ t :" 1l

t minp r Y p1q , r C p1q q 1´t ď minp r Y p2q , r
C p2q q t u and r δ pjq :" 1l t r Y pjq ď r C pjq u , j " 1, 2, the pair p r Z t , r δ t q can actually be observed.

Let p r Z t,i , r δ t,i q, i " 1, . . . , n, be independent copies of the random pair p r Z t , r δ t q. We are now in the classical univariate censoring framework and we want to propose a nonparametric robust estimator for A Y ptq by means of the MDPD criterion, adjusted to the censoring, i.e., we minimize for α ą 0 the function p ∆ α,t paq :"

ż 8 0 `ae ´az ˘1`α dz ´ˆ1 `1 α ˙ż 8 0 `ae ´az ˘α d p F KM r Yt pzq.
The MDPDE p A Y,α,n ptq for A Y ptq satisfies the estimating equation

p ∆ p1q α,t p p A Y,α,n ptqq " 0 (4)
where p ∆ pjq α,t p.q denotes the derivative of order j of p ∆ α,t p.q. Our aim in this paper is to show the joint convergence in distribution of

" ? n ´p A Y,α,n pt j q ´AY pt j q ¯, j " 1, . . . , m ı ,
where tt 1 , . . . , t m u is a grid of values in r0, 1s. In order to achieve this goal, a crucial step is the study of statistics of the type

T n pt, ξq :" ż 8 0 φ t,ξ pzqd p F KM r Yt pzq, (5)
with φ t,ξ pzq :" z ξ e ´αA Y ptqz , t P r0, 1s and ξ P N, as p ∆ α,t and its derivatives are essentially linear combinations of such statistics, see, e.g., Section 6.1.

Asymptotic properties

In this section we derive the asymptotic properties of our estimator p A Y,α,n ptq under suitable assumptions. As a first step we need to establish the limiting behavior of (5). Kaplan

T n pt, ξq " 1 n n ÿ i"1 η t,ξ p r Z t,i q `Rn,t,ξ where η t,ξ p r Z t,i q :" φ t,ξ p r Z t,i qγ p0q t p r Z t,i q r δ t,i `γp1q t,ξ p r Z t,i qp1 ´r δ t,i q ´γp2q t,ξ p r Z t,i q,
and R n,t,ξ " o P pn ´1{2 q.

From this representation, we can deduce the convergence in distribution of our key statistic T n .

Corollary 3.1 Under the assumption of Theorem 3.1

? n ˜Tn pt, ξq ´Γpξ `1q pα `1q ξ`1 A ξ Y ptq ¸d ÝÑ N p0, σ 2 pt, ξqq,
where

σ 2 pt, ξq :" V arpη t,ξ p r Z t qq,
and Γ is the gamma function defined as Γprq :" ş 8 0 t r´1 e ´tdt, @r ą 0.

We now derive the limiting distribution of a vector of statistics of the form (5), when properly normalized. Let T n and T be pm ˆ1q vectors defined as

T n :" pT n pt 1 , ξ 1 q, . . . , T n pt m , ξ m qq T , and
T :" ˜Γpξ j `1q pα `1q ξ j `1A ξ j Y pt j q , j " 1, . . . , m ¸T ,
for some positive integer m, where T stands for the transpose matrix. The aim of the next theorem is to provide the finite dimensional convergence result which will allow us to establish the convergence in distribution of our robust estimator of the Pickands dependence function A Y .

Theorem 3.2 Under the assumptions of Theorem 3.1, we have

? n pT n ´T q d ÝÑ N m p0, Σq ,
where N m denotes a m´dimensional normal distribution and Σ the pm ˆmq covariance matrix with elements pσ j,k q 1ďj,kďm :" pCovpη t j ,ξ j p r Z t j q, η t k ,ξ k p r Z t k qqq 1ďj,kďm .

Note that for the result of Theorem 3.2 we need to assume p1 `2αqA Y pt j q ´AC pt j q ą 0 for j " 1, . . . , m, which imposes a constraint on the parameter α. As a worst case scenario we could consider A Y ptq " maxtt, 1 ´tu (corresponding to complete dependence) and A C ptq " 1 (corresponding to independence), and require p1 `2αqA Y ptq ´AC ptq ą 0 for all t P r0, 1s. By some standard calculations one can easily obtain that this will be satisfied if α ą 0.5.

By using the above results we can now prove the existence of a consistent sequence of solutions to the estimating equation (4).

Theorem 3.3 Under the assumptions of Theorem 3.1, we have that with probability tending to one there exists a sequence p p A Y,α,n ptqq ně1 of solutions to (4), such that p A Y,α,n ptq P ÝÑ A Y ptq as n Ñ 8.

We have now all the needed ingredients for proving the finite dimensional weak convergence of the MDPDE for A Y on a grid tt 1 , . . . , t m u of positions in r0, 1s. Let A denote a pm ˆ2mq matrix with elements

A i,j :" $ & % p1 `αqrA Y pt i qs α´1 , if j " 2i ´1, ´p1 `αqrA Y pt i qs α , if j " 2i, 0, otherwise,
B is an pm ˆmq diagonal matrix with entries B i,i :" p1 `αq 2 rA Y pt i qs α´2 p1 `α2 q , and the matrix C is a p2m ˆ2mq matrix with elements C 2i´1,2j´1 :" Covpη t i ,0 p r Z t i q, η t j ,0 p r Z t j qq,

C 2i´1,2j :" Covpη t i ,0 p r Z t i q, η t j ,1 p r Z t j qq, C 2i,2j´1 :" Covpη t i ,1 p r Z t i q, η t j ,0 p r Z t j qq, C 2i,2j :" Covpη t i ,1 p r Z t i q, η t j ,1 p r Z t j qq.
Theorem 3.4 Under the assumptions of Theorem 3.1, we have

? n » - - p A Y,α,n pt 1 q ´AY pt 1 q . . . p A Y,α,n pt m q ´AY pt m q fi ffi fl d ÝÑ N `0, BACA T B T ˘. (6)
In particular, for t P r0, 1s, we have, as n Ñ 8,

? n ´p A Y,α,n ptq ´AY ptq ¯d ÝÑ N `0, q σ 2 ˘,
where q σ 2 :" N {D, and

N :" p1 `αq 2 A 3 Y ptq " p1 `4α `9α 2 `14α 3 `13α 4 `8α 5 `4α 6 qA 2 Y ptq ´2p1 `3α `5α 2 `6α 3 `4α 4 `2α 5 qA Y ptqA C ptq `p1 `2α `3α 2 `2α 3 `α4 qA 2 C ptq ‰ , D :" p1 `α2 q 2 rp1 `2αqA Y ptq ´AC ptqs 3 .

Simulation experiment

In this section we illustrate the finite sample performance of the proposed estimator with a small simulation study. In first instance we consider distributions for pY p1q , Y p2q q and pC p1q , C p2q q, with an extreme value copula and unit exponential margins. The contamination is introduced by the following mixture model F ε py 1 , y 2 q " p1 ´εqF py 1 , y 2 q `εF c py 1 , y 2 q,

where ε P r0, 1s represents the fraction of contamination in the dataset, F is the distribution function of pT p1q , T p2q q :" pminpY p1q , C p1q q, minpY p2q , C p2q qq and F c is the contamination distribution function. For the main distributions of pY p1q , Y p2q q and pC p1q , C p2q q we consider the asymmetric logistic distribution, with survival function G ' py 1 , y 2 q " exp ´´p1 ´ψ1' qy 1 ´p1 ´ψ2' qy 2 ´ppψ 1' y 1 q 1{r' `pψ 2' y 2 q 1{r' q r' ¯, y 1 , y 2 ą 0,

where ' denotes either Y or C, r ' P p0, 1s and ψ 1' , ψ 2' P r0, 1s. In this model independence is obtained for either r ' " 1, or ψ 1' " 0 or ψ 2' " 0, while complete dependence is obtained for ψ 1' " ψ 2' " 1 and r ' Ó 0. The logistic model is a special case of the asymmetric logistic model, and corresponds to ψ 1' " ψ 2' " 1. We consider the following settings

• Setting 1: pr Y , ψ 1Y , ψ 2Y q " p0.25, 1, 1q and pr C , ψ 1C , ψ 2C q " p0.75, 1, 1q,

• Setting 2: pr Y , ψ 1Y , ψ 2Y q " p0.75, 1, 1q and pr C , ψ 1C , ψ 2C q " p0.25, 1, 1q,

• Setting 3: pr Y , ψ 1Y , ψ 2Y q " p0.1, 0.4, 0.6q and pr C , ψ 1C , ψ 2C q " p0.05, 0.6, 0.4q.

These settings for the main distributions are then combined with the following types of contamination:

• First type of contamination: the distribution function F c is given by

F c py 1 , y 2 q " 1 2 1 ´e´y 1 `1 ´e´y 2 ( 1l ty 1 ě0,y 2 ě0u .
This means that the contamination is on the axes according to the unit exponential distribution.

• Second type of contamination: the distribution function F c has completely dependent unit exponential margins.

We simulate N " 100 datasets of size n " 1000, and consider ε " 0, 0.025 and 0.05. We estimate A Y ptq on the grid t0.05, 0.10, . . . , 0.95u for t. In Figures 2 till 7 we show the boxplots of the estimates p A Y,α,n ptq for t P t0.05, 0.10, . . . , 0.95u, together with A Y ptq (blue solid line) and A C ptq (green dashed line). In each of the figures the rows correspond to the levels of contamination, while the columns correspond to α " 0.1, 0.5 and 1, respectively. From these simulations we can draw the following conclusions • In case of no contamination, we can see that increasing α improves the estimation. The robustness of the MDPD method is thus partly used to correct for the censoring of the data.

• If we increase the contamination then estimation becomes more difficult, whatever α. Note that we can handle 5% contamination still reasonably well.

• For a given percentage of contamination we observe that increasing α from 0.1 to 0.5 clearly improves estimation but increasing α further to 1 does not lead to clear further improvements.

• The contamination affects the estimators in the expected direction: axes contamination pulls the estimator up and diagonal contamination down.

• For the three settings considered, contamination on the axes seems to be more difficult than contamination on the diagonal.

• In Setting 3 with asymmetric Pickands dependence functions, we recover A Y quite well, despite the fact that A Y and A C are rather close. q q q q q q q q q q q q q q 0.05 0. 1.0 t Ayt q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0. q q q q q q q q q q q q q q 0.05 0. 1.0 t Ayt q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q 0.05 0. 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0. 1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0. 1.0 t Ayt q q q q q q q q q q q q q q q q q q q q 0.05 0. 1.0 t Ayt q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q 0.05 0. 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt Figure 6: Setting 3, contamination on the axes. First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q 0.05 0. Next, we illustrate the situation where the marginal distributions are not unit exponential. We consider the case of a logistic Pickands dependence function for pY p1q , Y p2q q and pC p1q , C p2q q with r Y " 0.25 and r C " 0.75, respectively. The marginal distributions of Y p1q and Y p2q are Exp [START_REF] Beirlant | Bias reduced tail estimation for censored Pareto type distributions[END_REF], and those of C p1q and C p2q are Fréchet(2) shifted by 0.75, i.e. F C pjq pcq " expp´pc 0.75q ´1{2 q, c ą 0.75, j " 1, 2. This gives about 5% censoring. As before, we combine this setting with the two types of contamination, whereafter the observations are transformed to approximate unit exponential by ´log p F KM Y pjq , j " 1, 2. The results are shown in Figures 8 and9, which have a layout that is the same as before. In case of no contamination, we have that overall we can capture the shape of the Pickands dependence function but the estimate is biased downwards. Using α " 0.5 gives slightly better results than α " 0.1, especially in the centre of the range for t, and increasing α to one does not lead to further improvements. When adding contamination the estimates behave again as expected, in particular they become pulled up under axes contamination and pulled down under diagonal contamination. Using α " 0.5 gives some protection against contamination, in the sense that the results are close to those obtained under the uncontaminated case, and gives slightly less biased results than α " 1.

Data example

In this section we illustrate the nonparametric MDPDE for A Y on the dataset of insurance company indemnity claims introduced in Section 1.

In [START_REF] Frees | Understanding relationships using copulas[END_REF] these data were analysed by fitting parametric copula models to the data using the maximum likelihood method. Based on their findings, the Gumbel-Hougaard copula (corresponding to the logistic model considered here) provided the best fit to the data, which gives an indication that the distribution underlying the data has an extreme value copula. Also in [START_REF] Cebrián | Analysis of bivariate tail dependence using extreme value copulas: an application to the SOA medical large claims database[END_REF], a similar loss-ALAE dataset was modelled with extreme value copulas.

The nonparametric MDPDE of A Y is shown in Figure 10 for α " 0.1 (solid line), α " 0.5 (dotted line) and α " 1 (dashed-dotted line), along with the parametric maximum likelihood estimate obtained by Frees and Valdez (1998) (dashed line). Our nonparametric MDPDE with α " 0.1 is close to the parametric estimate from [START_REF] Frees | Understanding relationships using copulas[END_REF], and this can be seen as a further confirmation of the fit of the parametric Gumbel-Hougaard copula proposed in the latter. Note that both these estimates are not robust with respect to outliers. To overcome this we also applied the MDPD method with α " 0.5 and α " 1, and the obtained estimates differ slightly from the MDPDE with α " 0.1 and the [START_REF] Frees | Understanding relationships using copulas[END_REF] estimate. This might indicate the presence of a few outliers in the loss-ALAE data, as was also suggested by the scatterplots given in Figure 1, though these do not seem to have a disturbing effect on the estimation, which is reasonable given a total sample size of n " 1500. q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q 0.05 0. . First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1. q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q q 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 0.5 0.6 0.7 0.8 0.9

1.0 t Ayt q q q q q q q q q q 0.05 0. . First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1. 6 Appendix: Proofs of the results

Derivatives of the MDPD objective function

The arguments used to establish the consistency and asymptotic normality of the MDPDE depend on the MDPD objective function and its derivates. For convenience, in this section we give p ∆ pjq α,t paq, for j " 1, 2 and 3. Straightforward computations for α ą 0, give 3 with probability tending to one, where b :" d{3. Combining the above we find that with probability tending to 1, S 1 `S2 `S3 ą cr 2 ´p1 `bqr 3 , for a " A Y ptq ˘r. Clearly, since the right-hand side of the above inequality is positive if r ă c{p1 `bq, (9) follows.

A C ptq A C ptq `AY ptq ! 1 ´e´rA C ptq`A Y ptqsz ) , r H 1 t pzq " A Y ptq A C ptq `AY ptq ! 1 ´e´rA C ptq`A Y ptqsz ) , γ p0q
To complete the proof of the existence and consistency we adjust the line of argumentation of Theorem 3.7 in Chapter 6 of [START_REF] Lehmann | Theory of point estimation[END_REF]. For r ą 0, small enough such that A Y ptq ˘r P r1{2 ´δ, 1 `δs we let S n prq :" t p ∆ α,t pA Y ptqq ă p ∆ α,t paq for a " A Y ptq ˘ru.

From the above we have that P A Y ptq pS n prqq Ñ 1 for any such r, and hence there exists a sequence r n Ó 0 such that P A Y ptq pS n pr n qq Ñ 1 as n Ñ 8. By the differentiability of p ∆ α,t we have that v P S n pr n q implies that there exists a value p

  . Alternatively, one can use functions that give a complete characterisation of the extreme dependence like the Pickands dependence function, the stable tail dependence function and the spectral distribution function. Estimation of such dependence functions was considered in Capéraà et al. (1997), Fils-Villetard et al. (2008), Fougères et al. (2015) and Escobar-Bach et al. (2017). To the best of our knowledge censoring in the multivariate extreme value context is unexplored.

  15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0

Figure 2 :

 2 Figure2: Setting 1, contamination on the axes. First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

  15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0

Figure 3 :

 3 Figure3: Setting 1, contamination on the diagonal. First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

  15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 00.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 00.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0

Figure 4 :

 4 Figure4: Setting 2, contamination on the axes. First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

  15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 00.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 00.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0

Figure 5 :

 5 Figure5: Setting 2, contamination on the diagonal. First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

Figure 7 :

 7 Figure 7: Setting 3, contamination on the diagonal. First row: no contamination, second row: 2.5% contamination, third row: 5% contamination. First column: α " 0.1, second column: α " 0.5, third column: α " 1.

Figure 8 :

 8 Figure 8: Setting 1, contamination on the axes, margins transformed to unit exponential using ´log p F KM Y

Figure Setting 1 ,

 1 Figure Setting 1, contamination on the diagonal, margins transformed to unit exponential using ´log p F KM Y

Figure 10 :

 10 Figure 10: LOSS-ALAE data: nonparametric estimate pA Y,α,n for α " 0.1 (solid line), α " 0.5 (dotted line) and α " 1 (dashed-dotted line). The parametric Gumbel-Hougaard estimate of[START_REF] Frees | Understanding relationships using copulas[END_REF] is given by the dashed line.

  , Beirlant et al. (2007), Einmahl et al. (2008), Gomes and Neves (2011), Worms and Worms (2014), Beirlant et al. (2016), among others, where focus was mainly on estimating the extreme value index and extreme quantiles. Recently, extreme value regression problems with censoring were studied by Ndao et al. (2016), Stupfler (2016) and Goegebeur et al. (2018).

  Y ptq is a value between a and A Y ptq. ÝÑ 0, as n Ñ 8, and hence |S 1 | ă r 3 with probability tending to 1. As for S 2 , we have rA Y ptqs α´2 , and thus there exists a c ą 0 such that S 2 ą cr 2 with probability tending to one. For S 3 we use the fact that

	By Taylor's theorem
	p ∆ α,t paq ´p ∆ α,t pA Y ptqq " p ∆	p1q α,t pA Y ptqqpa ´AY ptqq	`1 2	p ∆	p2q α,t pA Y ptqqpa ´AY ptqq 2
						`1 6	p ∆	p3q α,t p q A Y ptqqpa ´AY ptqq 3
					": S 1 `S2 `S3 ,
	where q A Using Section 6.1 and the result of Corollary 3.1, we have p ∆	p1q α,t pA Y ptqq
			p ∆	p2q α,t pA Y ptqq p1 `αq 2 sup P 1 `α2 ÝÑ | p ∆ p3q α,t paq| ď M,
						aPr1{2´δ,1`δs
	where M
						ż 8	ω ξ e ´p1`αqA Y ptqω dω,
						z
	γ	p2q t,ξ pzq "	A Y ptqA C ptq A Y ptq `AC ptq	0 ż 8	ω ξ e ´p1`αqA Y ptqω
	Then, essentially we need to verify the conditions (1.5) and (1.6) in Stute (1995). In our context,
	as soon as p1 `2αqA Y ptq ´AC ptq ą 0, we have
						ż 8	rφ t,ξ pzqγ	p0q t pzqs 2 d r H 1 t pzq ă 8,
						0
			and	ż 8 0	|φ t,ξ pzq|	a	CpzqdF Yt pzq ă 8, r
	where
			Cpzq :"	ż z 0	dF Ct pyq r Zt pyqsr1 ´F r r1 ´F r Ct pyqs	.
	Using the fact that
			ż 8 0	φ t,ξ pωqdF	r Yt pωq "	Γpξ `1q Y ptq pα `1q ξ`1 A ξ	,

t pzq " exptA C ptqzu, γ p1q t,ξ pzq " A Y ptqe rA C ptq`A Y ptqsz ! e rA C ptq`A Y ptqs minpz,ωq ´1) dω. P P ÝÑ d ă 8. Thus |S 3 | ă br

  A Y,n P rA Y ptq ´rn , A Y ptq `rn s for which p ∆ α,t attains a local minimum, and thus p ∆ p1q α,t p p A Y,n q " 0. Now let p A Y,n :" p A Y,n for v P S n pr n q and arbitrary otherwise. ClearlyP A Y ptq p p ∆ α,t p p A Y,n q " 0q ě P A Y ptq pS n pr n qq Ñ 1,as n Ñ 8. Thus with probability tending to 1 there exists a sequence of solutions to the estimating equation (4). Then, for any fixed r ą 0 and n sufficiently largeP A Y ptq p| p A Y,n ´AY ptq|q ă rq ě P A Y ptq p| p A Y,n ´AY ptq|q ă r n q ě P A Y ptq pS n pr n qq Ñ 1,which establishes the consistency of the sequence p p A Y,n q ně1 .6.6 Proof of Theorem 3.4First we prove the result for a specific, single t P r0, 1s, since this is more explicit. Our starting point is the estimating equation (4). By applying a Taylor series expansion around the true value A Y ptq, we get 0 " p ∆ p1q α,t pA Y ptqq `´p A Y,α,n ptq ´AY ptq ¯p ∆ Y ptq is a random value between A Y ptq and p A Y,α,n ptq. A straightforward rearrangement of the above display gives ? n ´p A Y,α,n ptq ´AY ptqRemark that, according to Section 6.1, the numerator of the right-hand side of (10) can be rewritten asAfter some tedious computations we obtainV arpη t,0 p r Z t qq " α 2 A Y ptq p1 `αq 2 rp1 `2αqA Y ptq ´AC ptqs , V arpη t,1 p r Z t qq " p1 `2α `2α 4 qA 2 Y ptq ´2p1 `α ´α2 qA Y ptqA C ptq `A2 C ptq p1 `αq 4 A Y ptqrp1 `2αqA Y ptq ´AC ptqs 3 , Covpη t,0 p r Z t q, η t,1 p r Z t qq " ´αp1 `α ´α2 qA Y ptq ´αA C ptq p1 `αq 3 rp1 `2αqA Y ptq ´AC ptqs 2 . rA Y ptqs 2α´1 p1 `αq 2 rp1 `2αqA Y ptq ´AC ptqs 3 p1 `4α `9α 2 `14α 3 `13α 4 `8α 5 `4α 6 qrA Y ptqs 2 ´2p1 `3α `5α 2 `6α 3 `4α 4 `2α 5 qA Y ptqA C ptq `p1 `2α `3α 2 `2α 3 `α4 qrA C ptqs 2 ( . rA Y ptqs α´2 . A Y ptqq ´p A Y,α,n ptq ´AY ptq ¯P ÝÑ 1 `α2 p1 `αq 2 rA Y ptqs α´2 .Again by Slutsky's theorem, we have that Theorem 3.4 follows for a single t P r0, 1s. Now we focus on the general finite dimensional convergence. Let , 1q ´1 p1`αq 2 A Y pt 1 q ..

	p2q α,t pA Y ptqq ´?n p `1 2 ´p A Y,α,n ptq ´AY ptq ¯2 p ∆ ∆ p1q α,t pA Y ptqq α,t pA Y ptqq `1 2 ∆ p p2q p ∆ p3q α,t p r A Y ptqq ´p A Y,α,n ptq ´AY ptq p3q α,t p r A Y ptqq ¯. α,t pA Y ptqq " p1 `αqrA Y ptqs α´1 A ¯" where r ´?n p ∆ p1q " ? n ˆTn pt, 0q ´1 1 `α ȦY (10) ptq ? n ˆTn pt, 1q ´1 p1 `αq 2 A Y ptq ˙* . Then, by Theorem 3.2 ´?n p ∆ p1q α pA Y ptqq d ÝÑN p0, σ2 q, where σ2 :" Now, concerning the denominator of the right-hand side of (10), according to Section 6.1 and Corollary 3.1, p ∆ p3q α,t paq is bounded in probability and p ∆ p2q α,t pA Y ptqq P ÝÑ 1 `α2 p ∆ p2q α,t pA Y ptqq `1 2 p ∆ p3q α,t p r X n :" » ----------? n ´Tn pt 1 , 0q ´1 1`α ?n ´Tn pt 1 . ? n ´Tn pt m , 0q ´1 1`α ?n ´Tn pt m , 1q ´1 p1`αq 2 A Y ptmq ¯fi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl . Then D n :" ´» ---? n p ∆ ? n p ∆ Then p1 `αq 2 This implies that » p A Y,α,n pt 1 q ´AY pt 1 q ? n --. . .

p1q α,t 1 pA Y pt 1 qq . . . p1q α,tm pA Y pt m qq fi ffi ffi fl " AX n , and we have D n d

ÝÑN m p0, ACA T q. Now let B n be a pm ˆmq diagonal matrix with entries B n,i,i :"

1 p ∆ p2q α,t i pA Y pt i qq `1 2 p ∆ p3q α,t i p r A Y pt i qq ´p A Y,α,n pt i q ´AY pt i q ¯. p A Y,α,n pt m q ´AY pt m q fi ffi fl " B n D n ,

and the result follows.

a direct application of Theorem 1.1 in [START_REF] Stute | The central limit theorem under random censorship[END_REF] achieves the proof of Theorem 3.1.

Proof of Corollary 3.1

Similarly as Corollary 1.2 in [START_REF] Stute | The central limit theorem under random censorship[END_REF], Corollary 3.1 is a direct consequence of the result in our Theorem 3.1, combined with the classical central limit theorem and Slutsky's theorem.

Proof of Theorem 3.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

for all β P R m . A straightforward rearrangement of the terms leads to

According to Theorem 3.1, T 2,n " o P p1q as n Ñ 8. Now, concerning T 1,n , we have T 1,n "

n, where V 1 , . . . , V n are independent and identically distributed centered random variables with variance given by

This implies that V arpT 1,n q " β T Σβ where Σ " pCovpη t j ,ξ j p r Z t j q, η t k ,ξ k p r Z t k qqq 1ďj,kďm . By invoking the central limit theorem, T 1,n d ÝÑ N p0, β T Σβq. Now, by Slutsky's theorem we also have [START_REF] Fils-Villetard | Projection estimators of Pickands dependence functions[END_REF], and Theorem 3.2 follows.

Proof of Theorem 3.3

To prove the existence and consistency of p A Y,α,n ptq we adapt the proof of Theorem 5.1 in Chapter 6 of [START_REF] Lehmann | Theory of point estimation[END_REF], where existence and consistency of solutions of the likelihood equations is established, to the MDPDE framework. First we show that for any r sufficiently small

Note that r should be such that A Y ptq ˘r belongs to the parameter space, here r1{2 ´δ, 1 `δs, where δ ą 0.