Christophe Dutang
email: christophe.dutang@univ-lemans.fr

Yuri Goegebeur
email: yuri.goegebeur@imada.sdu.dk

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

Robust and bias-corrected estimation of the probability of extreme failure sets

Keywords: failure set, bias-correction, tail dependence, robustness, tail quantile process AMS 2000 Subject classifications. 62G32, 62H12, 62G20

published or not. The documents may come

Robust and bias-corrected estimation of the probability of extreme failure sets

Christophe Dutang, Yuri Goegebeur, Armelle

Introduction

Multivariate extreme value statistics deals with the estimation of the joint tail of a distribution function based on a random sample. Within this area, the problem of estimating the probability of extreme failure sets, i.e. regions of the sample space where one or several random variables are large, is of practical relevance for applications in many scientific disciplines. For instance, financial institutions are often interested in assessing the risk of simultaneous large negative returns on several stocks or other assets [START_REF] Charpentier | Limiting dependence structures for tail events, with applications to credit derivatives[END_REF]. In the design of coastal defence structures, engineers might want to estimate the probability that a seawall fails due to a combination of high waves and high still water levels (de Haan and de Ronde, 1998). The estimation of the probability of extreme failure sets is the topic of the present paper.

The extremal dependence between the components of a continuous random vector (X, Y ) with unit Fréchet margins (note that this can be assumed without loss of generality) can be analyzed with the model of [START_REF] Ledford | Modelling dependence within joint tail regions[END_REF]:

P(X > x, Y > y) = x -d 1 y -d 2 (x, y), x, y > 0,
where d 1 , d 2 > 0 and is a bivariate slowly varying function, i.e.

(tx, ty) (t, t) → ς(x, y) as t → ∞, for all x, y > 0, and the function ς is homogeneous of order zero. The parameter η := (d 1 + d 2 ) -1 is called the coefficient of tail dependence. It satisfies η ∈ (0, 1], and larger values of it indicate a stronger extremal dependence. As we can imagine, several attempts have been made to estimate η from data. Since

P(min(X, Y ) > z) = P(X > z, Y > z) = z -1/η (z, z),
i.e. the transformed variable min(X, Y ) follows a Pareto-type model with index 1/η, one can estimate η with classical estimators for the extreme value index like the Hill [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF], Pickands [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] or moment estimator [START_REF] Dekkers | A moment estimator for the index of an extreme-value distribution[END_REF]. However, this type of estimators typically suffers from bias and also they are not robust with respect to outliers. Moreover, since estimators for a bivariate tail probability are typically based on such estimators for η, one can expect that they will also be affected by bias and outliers. These specific issues will be addressed in the present paper.

In order to obtain a bias-corrected estimator we will, as usual in extreme value statistics, invoke a second order condition. In particular we will work under the following condition from Draisma et al.

(2004), which can be seen as an extension of the above discussed Ledford and Tawn model.

Condition SO: Let (X, Y ) be a random vector with joint distribution function F and continuous marginal distribution functions F X and F Y such that

lim t↓0 q 1 (t) -1 P (1 -F X (X) < tx, 1 -F Y (Y ) < ty) P (1 -F X (X) < t, 1 -F Y (Y ) < t) -c(x, y) =: c 1 (x, y) (1)
exists for all x ≥ 0, y ≥ 0 with x + y > 0, a function q 1 tending to zero as t ↓ 0, and c 1 a function neither constant nor a multiple of c. Moreover, we assume that the convergence is uniform on

{(x, y) ∈ [0, ∞) 2 |x 2 + y 2 = 1}.
Essentially, this condition is a second order multivariate regular variation condition on the function R(x, y) := P (1 -F X (X) < x, 1 -F Y (Y ) < y). It can be shown that R(t, t) is regularly varying at zero with index 1/η, |q 1 | is regularly varying at zero with index τ ≥ 0, and that the function c is homogeneous of order 1/η, that is c(tx, ty) = t 1/η c(x, y). Also, c 1 (x, x) = x 1/η (x τ -1)/τ .

Recently, several papers have addressed bias-corrected estimation of the coefficient of tail dependence η, but not too much attention has been paid to the estimation of the probability of extreme failure sets.

Concerning η, we refer to [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF] and [START_REF] Goegebeur | Asymptotically unbiased estimation of the coefficient of tail dependence[END_REF]. Also, Dutang et al.

(2014) introduced a robust and bias-corrected estimator for the coefficient of tail dependence, based on a submodel of (SO) that is fitted to the data by means of the minimum density power divergence technique. The estimators of probabilities of extreme failure sets already proposed in the literature are quite simple in nature, but their asymptotics is very difficult to handle and requires very technical conditions. This is for instance the case of the estimator proposed by [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF] where only the convergence in probability is established. [START_REF] Peng | Estimation of the coefficient of tail dependence in bivariate extremes[END_REF] also proposed an estimator but without showing convergence nor asymptotic normality. Moreover, bias-correction and robustness were not considered so far for estimators of probabilities of extreme regions. This constitutes the topic of the present paper.

The remainder of this paper is organised as follows. In the next section, we introduce the model and the estimator for probabilities of extreme failure sets. The asymptotic results of the proposed estimator are investigated in Section 3, and the finite sample performance is illustrated through an extensive simulation study in Section 4. Section 5 contains an application of the proposed methodology on a dataset from the actuarial context. The proofs of the main results can be found in the appendix.

Model assumptions and estimators

The robust and bias-corrected estimator for probabilities of extreme failure sets will be obtained from a submodel of the second order condition (SO), that will be fitted to properly transformed bivariate observations by means of the minimum density power divergence (MDPD) technique.

The density power divergence criterion was originally introduced by Basu et al. (1998) for the purpose of developing a robust estimation method. In particular, the density power divergence between arbitrary density functions f and h is given by

∆ α (f, h) :=    R h 1+α (z) -1 + 1 α h α (z)f (z) + 1 α f 1+α (z) dz, α > 0, R log f (z) h(z) f (z)dz, α = 0.
Note that for α = 0 one recovers the Kullback-Leibler divergence, whereas setting α = 1 leads to the L 2 divergence. Assume that the density function h depends on a parameter vector θ, and let f be the true density function of the random variable under consideration. The idea is then to estimate θ by minimizing an empirical version of ∆ α based on a random sample Z

1 , . . . , Z n from f : if α > 0 one considers ∆ α (θ) := R h 1+α (z)dz -1 + 1 α 1 n n i=1 h α (Z i ), whereas for α = 0 ∆ 0 (θ) := - 1 n n i=1 log h(Z i ).
For α = 0, one fits the model h to the data using the maximum likelihood method. The parameter α controls the trade-off between efficiency and robustness of the MDPD estimator: the estimator becomes more efficient but less robust against outliers as α gets closer to zero, whereas for increasing α the robustness increases and the efficiency decreases. Let (X, Y ) be a bivariate random vector with continuous marginal distributions satisfying

P (1 -F X (X) < x, 1 -F Y (Y ) < y) = x d 1 y d 2 g(x, y) 1 + 1 η δ(x, y) , x ≥ 0, y ≥ 0, (2)
where d 1 , d 2 are positive constants, η := (d 1 + d 2 ) -1 ∈ (0, 1) is the tail dependence coefficient, g is a continuous function that is homogeneous of order 0 and δ is a function of constant sign in the neighbourhood of zero, with |δ| being a bivariate regularly varying function, that is, there exists a

function ξ such that lim t↓0 |δ|(tx, ty) |δ|(t, t) = ξ(x, y), (3)
for all x, y ≥ 0. We assume additionally that ξ is continuous, homogeneous of order τ > 0, and that the convergence is uniform on {(x, y) ∈ [0, ∞) 2 |x 2 + y 2 = 1}. Note that we exclude the case η = 1, as was also done in [START_REF] Beirlant | Some comments on the estimation of a dependence index in bivariate extreme value in statistics[END_REF], [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF], and [START_REF] Goegebeur | Asymptotically unbiased estimation of the coefficient of tail dependence[END_REF].

Lemma 1 Model (2) satisfies assumption (SO) with c(x, y) = x d 1 y d 2 g(x, y) g(1, 1) , (4)
q 1 (t) = τ η δ(t, t), (5)
c 1 (x, y) = c(x, y) ξ(x, y) -1 τ . (6)
Many commonly used joint distribution functions satisfy model [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF], as for instance the Farlie Gumbel Morgenstern or the Frank distributions, see Section 4.

For convenience we assume that the marginal distributions are unit Pareto. In this case model [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF] becomes

P (X > x, Y > y) = x -d 1 y -d 2 g * (x, y) 1 + 1 η δ * (x, y) (7)
where g * (x, y) := g(1/x, 1/y) and δ * (x, y) := δ(1/x, 1/y). Note that the choice of unit Pareto margins is not only for convenience in order to use some results from

P(X > x, Y > y) = P X > x, ω 1 -ω Y > x ,
where ω := x/(x + y), ω ∈ (0, 1), can be interpreted as being a radial parameter and it has been introduced in order to estimate probabilities in joint tail regions. This ray parametrization was also used in [START_REF] Ramos | A new class of models for bivariate joint tails[END_REF]. This parameter is of little practical relevance for estimation of the coefficient of tail dependence, but it will play an important role for failure set probability estimation.

Thus, we consider the transformed variable Z ω := min(X, ω 1-ω Y ). For this variable one easily derives the survival function, given by

F Zω (z) := P(Z ω > z) = C ω z -1/η 1 + 1 η δ ω (z) , z > 0,
where Formally, using the common reparametrization ρ = -τ η, for u large

C ω := (ω/(1 -ω)) d 2 g * (1, (1 -ω)/ω),
P(Z ω > uz) P(Z ω > u) ≈ H(z; η, δ ω (u), ρ) :=    1, z ≤ 1, [z(1 + δ ω (u) -δ ω (u)z ρ/η )] -1/η , z > 1, (8)
where H has the density function

h(z; η, δ ω (u), ρ) =    0 z ≤ 1, 1 η z -1/η-1 [1 + δ ω (u)(1 -z ρ/η )] -1/η-1 1 + δ ω (u) 1 -1 + ρ η z ρ/η , z > 1,
(η ∈ (0, 1), ρ < 0, and δ ω (u) > max{- Specifically, for a sample (X 1 , Y 1 ), . . . , (X n , Y n ) of independent random vectors from model (2), one transforms into (approximate) unit Pareto margins by using the empirical distribution functions of the X and Y observations. This gives

Z ω,i := min n + 1 n + 1 -R X i , ω 1 -ω n + 1 n + 1 -R Y i with R X
i and R Y i denoting the rank of X i and Y i , i = 1, . . . , n, in the respective samples. The parameters η and δ ω of the extended Pareto distribution are estimated by fitting the density function h to the relative excesses Z j := Z ω,n-m+j,n / Z ω,n-m,n , j = 1, . . . , m, where 1 ≤ m ≤ n -1, and Z ω,1,n ≤ . . . ≤ Z ω,n,n are the order statistics of Z ω,1 , . . . , Z ω,n , using the minimum density power divergence criterion. Remember that the parameter δ ω reflects in fact the function δ ω (u), where δ ω (u) → 0 as u → ∞, but we do not make this dependence on the threshold explicit in the notation.

The MDPD estimator for η and δ ω satisfies the estimating equations

0 = ∞ 1 h α (z; η, δ ω , ρ) ∂h(z; η, δ ω , ρ) ∂η dz - 1 m m j=1 h α-1 (Z j ; η, δ ω , ρ) ∂h(Z j ; η, δ ω , ρ) ∂η , (9) 0
= ∞ 1 h α (z; η, δ ω , ρ) ∂h(z; η, δ ω , ρ) ∂δ ω dz - 1 m m j=1 h α-1 (Z j ; η, δ ω , ρ) ∂h(Z j ; η, δ ω , ρ) ∂δ ω . ( 10
)
Note that only the parameters η and δ ω are estimated with the MDPD method. The parameter ρ will in this paper be fixed at the correct value, or fixed at some arbitrary value, usually the canonical choice ρ = -1, and thus possibly mis-specified. Several estimators for ρ were proposed in the uni- Therefore, these authors did not use it in their simulation experiment, and also resorted to fixing the second order parameter at some value. We will address the estimation of second order parameters in the bivariate extreme value context in future research.

Our goal in this paper is to estimate the tail probability p n := P(X > z n , Y > y n ), with X and Y being unit Pareto random variables, and where z n → ∞ as n → ∞, and y n = ωz n for some ω > 0, i.e.

we estimate a tail probability along a ray. We have thus that

p n := F Zω (z n ), where ω = (1 + ω) -1 .
To this aim, let m be an intermediate sequence

, i.e. m → ∞ as n → ∞ with m = o(n). Assume that p n satisfies np n m -→ β ∈ [0, 1), n → ∞. (11)
Using the estimators η n and δ ω,n satisfying the estimating equations ( 9)-( 10), we can construct the following estimator for p n based on [START_REF] Dekkers | A moment estimator for the index of an extreme-value distribution[END_REF], setting u = Z ω,n-m,n and replacing F Zω (u) by the empirical proportion m/n:

p n := m n H z n Z ω,n-m,n ; η n , δ ω,n , ρ (12)
where ρ is either the true value of ρ or a mis-specified one.

Note that, at first sight, once the minima Z ω are being considered, our problem of estimation might appear a standard problem of tail estimation in the classical univariate framework. However, this is not the case: indeed, the Z ω,i , i = 1, ..., n, are not independent because of the rank transformation, and hence the theoretical study of their order statistics is very complicated, involving convergence results for stochastic processes. All subsequent results are proved in Section 3 using the stochastic process representation, which requires arguments that are different from the univariate independent and identically distributed (i.i.d.) case.

Asymptotic results

Consider the random variables Z ω,1 , . . . , Z ω,n , with order statistics Z ω,1,n ≤ . . . ≤ Z ω,n,n . In order to obtain the limiting distribution of p n , we need some preliminary results given in Theorems 1 and 2. These theorems are based on Theorem 1 in Dutang et al. (2014) which is for completeness recalled in the appendix. In the following, we denote the tail quantile process as

Q n (t) := Z ω,n-[mt],n , 0 < t < n/m. Set c ω := c(1, ω/(1 -ω)), ξ ω := ξ(1, ω/(1 -ω)), q(t) := P(1 -F X (X) < t, 1 -F Y (Y ) < t)
and let k/n = q ← (m/n), where it is assumed that k → ∞.

The estimating equations ( 9) and ( 10) depend only on the data through the following statistics

           A (1) m,n (s 1 ) := 1 m m j=1 Z ω,n-j+1,n Z ω,n-m,n s 1 = 1 0 Qn(t) Qn(1) s 1 dt A (2) m,n (s 2 ) := 1 m m j=1 Z ω,n-j+1,n Z ω,n-m,n s 2 log Z ω,n-j+1,n Z ω,n-m,n = 1 0 Qn(t) Qn(1) s 2 log Qn(t) Qn(1) dt A (3) m,n (s 3 ) := 1 m m j=1 Z ω,n-j+1,n Z ω,n-m,n s 3 log Z ω,n-j+1,n Z ω,n-m,n 2 = 1 0
Qn(t) Qn(1) s 3 log Qn(t)
Qn( 1)

2 dt
with s 1 , s 2 , s 3 ≤ 0 whose joint asymptotic normality follows from Theorem A (cf. Appendix). However, in order to derive the limiting distribution of p n defined in [START_REF] Drees | On smooth statistical tail functionals[END_REF], we need to add an additional component, Q n (1) = Z ω,n-m,n , correctly rescaled. The aim of Theorem 1 is thus to establish the joint limiting distribution of the main statistics. Here and in the sequel the arrow denotes convergence in distribution.

Theorem 1 Let (X 1 , Y 1 ), . . . , (X n , Y n ) be independent copies of the random vector (X, Y ) which has a joint distribution satisfying (2) such that the function c given in (4) has continuous first order partial

derivatives. For m, k → ∞ as n → ∞ such that √ mq 1 (k/n) → λ ∈ R we have √ m            k n Q n (1) -c η ω A (1) m,n (s 1 ) -1 1-ηs 1 A (1) m,n (s 1 ) -1 1-ηs 1 A (2) m,n (s 2 ) - η (1-ηs 2 ) 2 A (3) m,n (s 3 ) -2η 2 (1-ηs 3 ) 3            N 5 (λµ, Σ)
where

µ :=                      -η 2 c η ω c ρ ω ξω-1 ρ -ξ ω c ρ ω s 1 η 2 (1-ηs 1 )(1-ρ-ηs 1 ) -ξ ω c ρ ω s 1 η 2 (1-ηs 1 )(1-ρ-ηs 1 ) -ξ ω c ρ ω η 2 (1-ρ-η 2 s 2 2 ) (1-ηs 2 ) 2 (1-ρ-ηs 2 ) 2 ξ ω c ρ ω 2η 3 ρ 1 (1-ηs 3 ) 3 - 1-ρ (1-ρ-ηs 3 ) 3                     
, and where Σ has elements given by (1989). Indeed in the univariate case this property can be explained by the fact that for U 1 , . . . , U n i.i.d. U (0, 1) random variables, with order statistics U 1,n ≤ . . . ≤ U n,n , one has that V j,k := U j,n /U k+1,n , j = 1, . . . , k, behave jointly as the order statistics of a random sample of size k from the U (0, 1) distribution, independently from U k+1,n .

σ 11 := η 2 c 2η ω σ 22 := η 2 s 2 1 (1 -s 1 η) 2 (1 -2s 1 η) σ 33 := η 2 s 2 1 (1 -s 1 η) 2 (1 -2s 1 η) σ 44 := η 2 2 (1 -2s 2 η) 3 - 1 (1 -s 2 η) 4 σ 55 := η 4 4! (1 -2s 3 η) 5 - 4 (1 -s 3 η) 6 σ 12 = σ 13 = σ 14 = σ 15 := 0 σ 23 := η 2 s 1 s 1 (1 -s 1 η)(1 -s 1 η)[1 -(s 1 + s 1 )η] σ 24 := η 1 [1 -(s 1 + s 2 )η] 2 - 1 (1 -s 1 η)(1 -s 2 η) 2 σ 25 := 2η 2 1 [1 -(s 1 + s 3 )η] 3 - 1 (1 -s 1 η)(1 -s 3 η) 3 σ 34 := η 1 [1 -(s 1 + s 2 )η] 2 - 1 (1 -s 1 η)(1 -s 2 η) 2 σ 35 := 2η 2 1 [1 -(s 1 + s 3 )η] 3 - 1 (1 -s 1 η)(1 -s 3 η) 3 σ 45 := 2η 3 3 [1 -(s 2 + s 3 )η] 4 - 1 (1 -s 2 η) 2 (1 -s 3 η) 3 .

The proof of this

Note however that the formal establishment of the asymptotic independence property in our situation is different from the classical univariate i.i.d. framework since we work with the order statistics of the (non-independent) variables Z ω,i , i = 1, . . . , n. In our situation, the independence can be intuitively explained by the fact that the asymptotic distribution of (A

(1) m,n (s 1 ), A (1) m,n (s 1 ), A (2) m,n (s 2 ), A (3) m,n (s 3 )) is totally determined by that of Zω,n,n Z ω,n-m,n , ..., Z ω,n-m+1,n Z ω,n-m,n
whereas the asymptotic distribution of the first component in our Theorem 1, Q n (1), is totally determined by that of Z ω,n-m,n .

From now on, we will denote the true value of η and ρ by η 0 and ρ 0 , respectively. Let δ ω,n := δ ω ( Z ω,n-m,n ). In the next theorem, we establish the joint limiting distribution of Q n (1), η n and δ ω,n , when properly normalised. The parameter ρ in the estimating equations ( 9) and ( 10) is denoted in the sequel as ρ and is either the true value ρ 0 or a mis-specified one.

Theorem 2 Under the conditions of Theorem 1 we have that

√ m      k n Q n (1) -c η 0 ω η n -η 0 δ ω,n -δ ω,n           Z Γ ∆      (13)
where

     Z Γ ∆      ∼ N 3 λ[C -1 ( ρ)B( ρ) µ -δ * ], C -1 ( ρ)B( ρ)D( ρ)B ( ρ)C -1 ( ρ) , (14)
with

µ :=               -η 2 0 c η 0 ω c ρ 0 ω ξω-1 ρ 0 ξ ω c ρ 0 ω αη 0 (1+η 0 ) (1+α(1+η 0 ))(1-ρ 0 +α(1+η 0 )) ξ ω c ρ 0 ω η 0 [α(1+η 0 )-ρ] (1-ρ+α(1+η 0 ))(1-ρ 0 -ρ+α(1+η 0 )) -ξ ω c ρ 0 ω η 2 0 [1-ρ 0 -α 2 (1+η 0 ) 2 ] (1+α(1+η 0 )) 2 (1-ρ 0 +α(1+η 0 )) 2               , δ * := [0, 0, -η 2 0 ξ ω c ρ 0 ω /ρ 0 ], B( ρ) :=      1 0 0 0 0 -η -α-1 0 0 η -α-2 0 0 -η -α-1 0 η -α-1 0 (1 -ρ) 0      , the symmetric matrix C( ρ) has elements c 11 ( ρ) := 1 c 12 ( ρ) := 0 c 13 ( ρ) := 0 c 22 ( ρ) := η -α-2 0 1 + α 2 (1 + η 0 ) 2 [1 + α(1 + η 0 )] 3 c 23 ( ρ) := η -α-2 0 ρ(1 -ρ)[1 + α(1 + η 0 ) + α 2 (1 + η 0 ) 2 ] + α 3 ρ(1 + η 0 ) 3 [1 + α(1 + η 0 )] 2 [1 -ρ + α(1 + η 0 )] 2 c 33 ( ρ) := η -α-2 0 (1 -ρ) ρ 2 + α ρ 2 (1 + η 0 )[α(1 + η 0 ) -ρ] [1 + α(1 + η 0 )][1 -ρ + α(1 + η 0 )][1 -2 ρ + α(1 + η 0 )]
and the symmetric matrix D( ρ) has elements

d 11 ( ρ) := η 2 0 c 2η 0 ω d 12 ( ρ) = d 13 ( ρ) = d 14 ( ρ) := 0 d 22 ( ρ) := α 2 (1 + η 0 ) 2 [1 + α(1 + η 0 )] 2 [1 + 2α(1 + η 0 )] d 23 ( ρ) := α(1 + η 0 )[α(1 + η 0 ) -ρ] [1 + α(1 + η 0 )][1 -ρ + α(1 + η 0 )][1 -ρ + 2α(1 + η 0 )] d 33 ( ρ) := [α(1 + η 0 ) -ρ] 2 [1 -ρ + α(1 + η 0 )] 2 [1 -2 ρ + 2α(1 + η 0 )] d 24 ( ρ) := η 0 1 [1 + 2α(1 + η 0 )] 2 - 1 [1 + α(1 + η 0 )] 3 d 34 ( ρ) := η 0 1 [1 -ρ + 2α(1 + η 0 )] 2 - 1 [1 + α(1 + η 0 )] 2 [1 -ρ + α(1 + η 0 )] d 44 ( ρ) := η 2 0 2 [1 + 2α(1 + η 0 )] 3 - 1 [1 + α(1 + η 0 )] 4 .
Note that if ρ = ρ 0 in ( 14), the mean of the limiting normal distribution is (-

λη 2 0 c η 0 ω c ρ 0 ω ξω-1 ρ 0 , 0, 0) .
Thus, the estimator η n is asymptotically unbiased in the sense that the mean of the limiting distribution is zero, whatever the value of λ. Obviously, if ρ is not the true value, this property of unbiasedness is lost. Nevertheless, despite the loss of asymptotic unbiasedness when fixing ρ at a mis-specified value, like e.g. -1, one has in the univariate context often observed that the resulting second order estimators still have some bias-correcting effect and typically outperform estimators that are not corrected for bias, both in terms of bias and mean squared error; see e.g. [START_REF] Gomes | Bias reduction and explicit semi-parametric estimation of the tail index[END_REF], [START_REF] Feuerverger | Estimating a tail exponent by modelling departure from a Pareto distribution[END_REF]. A similar observation was made by [START_REF] Goegebeur | Asymptotically unbiased estimation of the coefficient of tail dependence[END_REF] in the context of biascorrected estimation of η. Keep in mind that in Theorem 2,

Z = η 0 c η 0 ω W (1/c ω ) -λη 2 0 c η 0 ω (c ρ 0 ω ξ ω -1)/ρ 0 ,
where W (1/c ω ) comes from Theorem A, and satisfies W (1/c ω ) ∼ N (0, 1). Now we derive the asymptotic normality of the estimator p n . In first instance we consider the situation where np n /m → β ∈ (0, 1).

Theorem 3 Under the conditions of Theorem 1 and (11) with β ∈ (0, 1), we have

√ m p n p n -1 W 1 c ω - ln β η 0 Γ - 1 -β -ρ η 0 ∆ + λ η 0 ξ ω c ρ 0 ω ρ 0 (β -ρ 0 -β -ρ ).
The result of Theorem 3 indicates that the tail probability estimator is asymptotically unbiased if one uses the correct value for ρ. We now consider the case where np n /m → 0, corresponding to a more

η 0 √ m ln d n p n p n -1 Γ.

Simulation study

In this section, we illustrate the finite sample properties of our robust and bias-corrected estimator p n through a simulation study. In order to make the dependence on the tuning parameter α explicit, we use from now on the notation p n,α .

Note that model (2), under which our estimator has been constructed, is in fact a condition on the copula function. Indeed

P (1 -F X (X) < x, 1 -F Y (Y ) < y) = x + y -1 + C(1 -x, 1 -y).
As already mentioned, this condition is satisfied by many commonly used bivariate distribution functions. In particular, we illustrate the behaviour of our estimator p n,α on the following two well-known copulas:

• the Farlie Gumbel Morgenstern (FGM) copula function given by

C(x, y) = xy [1 + ζ(1 -x)(1 -y)] , (x, y) ∈ [0, 1] 2 , ζ ∈ [-1, 1].
We take two different values of ζ: ζ = 1 which implies η = 1/2, and ζ = -1 which leads to η = 1/3.

For the FGM copula τ = 1.

• the Frank copula function given by

C(x, y) = - 1 2 log 1 - (1 -e -2x )(1 -e -2y ) 1 -e -2 , (x, y) ∈ [0, 1] 2 .
In that case η = 1/2 and τ = 1.

In our simulation, the value of z is chosen such that the (theoretical) survival probability p z, ω = P(X > z, Y > ωz) is equal to 2%, 0.1% and 0.01%, see Table 1.

FGM Frank FGM ζ = 1 ζ = -1 ζ = 1 ω = 0.
= P(X > z, Y > ωz) with ω = (1 + ω) -1 .
The data are generated as follows. First, we simulate n pairs (X i , Y i ), i = 1, ..., n, independently from a FGM copula with unit Fréchet marginals. Similarly, we simulate n pairs from a Frank copula. Each time, we transform the margins into (approximately) unit Pareto using the empirical distribution functions. This gives us Z ω,i , i = 1, ..., n. Finally, we minimize the empirical density power divergence ∆ α .

In each setting (FGM or Frank copulas), we also contaminate the sample as follows:

• we simulate independently Xj and Yj , j = 1, ..., n 0 := nε , from a unit Fréchet distribution;

• all our above methodology is applied to our new n + n 0 pairs (X 1 , Y 1 ), ..., (X n , Y n ), ( X1 + X n,n , Y1 +

Y n,n ), ..., ( Xn 0 +X n,n , Yn 0 +Y n,n ), where X n,n and Y n,n are sample maxima of the X and Y observations, respectively.

The percentage of contamination is set to ε = 0, 2, 5%, while n = 100, thus n 0 = 0, 2, 5. The procedure is repeated 2 10 = 1024 times since we used parallelization on a 8-core computer.

The left panel of Figure 1 represents the mean of our estimator p n,α of p z,1 = 2% and the right panel the mean squared error of the ratio p n,α /p z,1 as a function of m for the FGM copula with ζ = 1 and different values of α: 0 (solid line), 0.1 (dashed line), 0.5 (dotted line), 1 (dash-dotted line). Three levels of contamination have been considered: ε = 0 (first row), 2% (second row), 5% (third row).

The parameter ρ is mis-specified to the value -1. Analogously, Figures 2 and3 are built for the FGM copula with ζ = -1 and the Frank copula, respectively, both in case p z,1 = 2%. Note that in these figures, ω is set to the value 1, i.e. ω to 1/2. Based on our simulations, we can draw the following conclusions:

• on uncontaminated simulations, the choice α = 0 is the best one, whatever the distribution is. This is expected, since in this no contamination framework, a procedure of robustification is not necessary.

However, if we keep α small, we preserve the good behaviour of the estimator, in terms of bias and MSE;

• when contamination occurs, we clearly observe the superiority of our estimator p n,α for α > 0. In terms of bias as well as MSE, α = 0.5 seems to be nearly always the best choice;

• as expected, increasing the percentage of contamination deteriorates the estimation procedure;

• the superiority of our estimator p n,α on p n,0 is stronger in case FGM with ζ = -1 compared to FGM with ζ = 1 in terms of bias and MSE (see the different scale on the y-axis in Figures 1 and2). This can be explained by the fact that one disturbs the sample by outliers that have a different dependence structure in Figure 2 where η = 1/3, whereas for Figure 1 where η = 1/2, the data are only disturbed by outliers with a similar dependence structure as for the uncontaminated data;

• when a proper choice of α is used (all considered values of α in case of no contamination and α = 0.5 or 1 when contamination is present), the estimators are clearly very stable as a function of m, and stay close to the true value of p z,1 for a wide range of values for m, which illustrates the bias-correction of our procedure. We also run our simulation study for smaller values of p z,1 (out of the sample, see the right panel in Table 1) and larger sample size. The results are displayed in Figures 4 to 6. To keep the length of the paper reasonable, we only show the results for the FGM distribution with ζ = 1, which is, as explained above, the more challenging case. From the simulation we deduce:

• when contamination occurs, the superiority of our estimator p n,α for α > 0 is again observed. In terms of bias as well as MSE, α = 0.5 seems to be nearly always the best choice;

• the estimators are again very stable as a function of m illustrating the bias-correction of our procedure;

• the MSEs are larger than in Figures 123due to the fact that such estimation is much more challenging, since they are out of the sample.

As usual in the robust framework, we compute also the breakdown point of our estimator for contamination as follows:

• we simulate 1024 samples of size n (100 and 1000) and take a percentage of contamination ε = 0, 1%, ..., 10%;

• for each value of ε and each sample size n, we compute our estimator p n,α of p z,1 = 2% and the mean squared error of the ratio p n,α /p z,1 as a function of m;

• for a fixed value of m, the breakdown point is defined as the smallest value of ε such that the MSE exceeds 1. This breakdown value is reported in Table 2 for n = 100 and in Table 3 for n = 1000. Each table contains three subtables corresponding to the three distributions under consideration (FGM with ζ = 1, -1, and Frank) and for each distribution we consider four values of α (0, 0.1, 0.5 and 1). For both sample sizes and whatever the distribution, the higher the value of α, the larger is the breakdown point. In particular, for α = 0.5 or 1, the breakdown point is particularly high illustrating the robustness of our estimator p n,α , whereas for α = 0, the breakdown point is low.

In order to examine the accuracy of the asymptotic results in finite samples we calculated the coverage probabilities of 95% confidence intervals for p n = 0.01% based on the result of Theorem 4, for the Farlie Gumbel Morgenstern copula with ζ = 1 in case of a sample of size n = 1000. The results are displayed in Table 4. Note that the coverage probabilities are generally somehow below 95% but can be considered fairly good, especially if one takes into account that the use of the result of Theorem 4 is not so obvious in practice since in fact p n → 0 when n → ∞, together with some additional assumptions on the involved sequences. Furthermore, we consider the estimation of our probability p z, ω , but this time outside the diagonal, that is for ω = 1. Figure 7 illustrates the behaviour of our estimator p n,α in the uncontaminated case for the FGM copula with ζ = 1, the same values of α, and two values of ω: 0.45 (first row) and 0.55 (second row). In both cases, we can observe the almost best behaviour of p n,0 compared to p n,α , α > 0, as expected in this framework. However, the most notable point compared to Figure 1 is that the estimation is much more difficult, with a quite stable bias appearing. This can be partly explained by our Theorems 3 and 4. Indeed, if we look at the bias term of our estimator p n,α , we have to compute the expectations of Γ and ∆ which can be deduced from our Theorem 2. Straightforward computations lead to an expectation for both variables of the form λξ ω c ρ 0 ω f (η 0 , ρ 0 , α) where the function f (η 0 , ρ 0 , α)

α = 0 2 1 1 1 1 1 1 1 1 FGM α 0.1 2 2 2 2 2 2 2 3 3 ζ = -1 α = 0.
α = 0 2 1 1 1 1 1 1 1 1 FGM α = 0.1 2 2 2 2 2 2 2 2 3 ζ = -1 α = 0.5 3
does not depend on ω. As the sense of variation of the function ω → ξ ω c ρ 0 ω depends if ω is smaller or larger than 1/2 for our three bivariate distribution functions, the estimation on the diagonal is simpler than outside.

Finally, to illustrate that mis-specifying the parameter ρ is not crucial, we show in Figure 8 the behaviour of our estimator p n,α for the FGM copula with ζ = 1 in the uncontaminated framework in case ρ is replaced by the true ρ 0 = -1/2. The same values of α have been used and ω is set to 1/2.

This figure is indeed very similar to the first row of Figure 1 where ρ is mis-specified to -1.

Illustration on an actuarial dataset

We consider an actuarial dataset published by the Australian regulator of insurance services which consists of company-wide performance indicators for insurers operating in Australia between 2005 and 2010, see http://www.apra.gov.au/. This dataset contains n = 422 observations of (gross) loss ratios X and expense ratios Y of 173 private insurance companies. Possible time dependency among the data points is ignored in this paper. 5. Since X and Y are ratios (in terms of gross written premium), we may reasonably qualify these points (X i , Y i 1) as outliers. In the sequel, we will denote by X and Y the original random variables, and by X and Y the random variables on the unit Pareto scale, as obtained from the inverse probability integral transform. Now, in order to detect if these three points listed in Table 5 are outliers or not, we consider the Deleting these three points clearly lead to estimators close to each other, for m not too small, whatever the value of α is. This corroborates our feeling that the three points are outliers. Indeed, as observed in Figure 10 (a) the estimator p n,α differs considerably depending on the value of α which is expected in the presence of outliers. As in the simulation study, the estimator p n,0.5 remains stable both in the situation with and without outliers.

Finally, Figure 11 shows the median over m = 50, . . . , 80 of our probability estimates p n,α of P( X > z, Y > z) for increasing values of z and different values of α. Again, two cases have been considered: the whole dataset in panel (a) and the dataset without the three points in panel (b). The same conclusions follow, i.e. the estimators are close to each others if the three points are removed from the initial dataset and p n,0.5 exhibits a more stable pattern than p n,α for other values of α in case of the whole dataset. Both numerical illustrations of this section advocate α = 0.5 as a good value, in line with the results of Section 4. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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 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 0 1000 2000 3000 4000

qq 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 q q q q q q q q q q q q q 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 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 qq 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 qq q q q q q q q q q q q q q q q q q qq 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 qq 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 qq qq 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 qq 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 0 100 200 300 400 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q 0 1 2

Appendix

The key element of our proofs is the following theorem established in [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF].

Theorem A [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF]. Under the assumptions of Theorem 1, we have that there exist suitable versions of Q n , a suitable process W (./c ω ), equal in distribution to a standard Brownian motion, such that for all t 0 , ε > 0

sup 0<t≤t0 t η+ 1 2 +ε √ m k n Q n (t) - t c ω -η -ηt -(η+1) c η ω W t c ω + λη 2 t c ω -η t cω -ρ ξ ω -1 ρ = o P (1).
Proof of Theorem 2 We follow the lines of proof of Theorem 2 in Dierckx et al. (2013). We give here only the adjustements needed. Define

A m,n ( ρ) := √ m          k n Q n (1) -c η 0 ω A (1) m,n -α(1+η 0 ) η 0 - 1 1+α(1+η 0 ) A (1) m,n -α(1+η 0 )-ρ η 0 - 1 1-ρ+α(1+η 0 ) A (2) m,n -α(1+η 0 ) η 0 - η 0 (1+α(1+η 0 )) 2          .
We can easily infer from Theorem 1 that A m,n ( ρ) A( ρ) ∼ N 4 (λ µ, D( ρ)). From the proof of Theorem

√ m      k n Q n (1) -c η 0 ω η n -η 0 δ ω,n      C -1 ( ρ)B( ρ)A( ρ) ∼ N 3 (λC -1 ( ρ)B( ρ) µ, C -1 ( ρ)B( ρ)D( ρ)B ( ρ)C -1 ( ρ)).
Taking into account that √ mδ ω,n P -→ -λη 2 0 ξ ω c ρ 0 ω /ρ 0 , Theorem 2 follows. From Theorem A, we easily infer that

k n Z ω,n-m,n = c η 0 ω + η 0 √ m c η 0 ω W 1 c ω - λ √ m η 2 0 c η 0 ω c ρ 0 ω ξ ω -1 ρ 0 + o P 1 √ m .
Thus, if we denote by Q Zω the quantile function of the variable Z ω and if U k,n denotes the k-th order statistic of a random sample of n uniform (0, 1) random variables, we deduce that

+ o P 1 √ m C η 0 ω q( k n ) -η 0 1 + [ n m U m+1,n -1] -η 0 (1 + δ ω (C η 0 ω U -η 0 m+1,n )(1 + o P (1)))
.

Recall now that

q k n = P 1 -F X (X) < k n , 1 -F Y (Y ) < k n = P X > n k , Y > n k = C 1 2 n k -1/η 0 1 + 1 η 0 δ 1 2 n k = g(1, 1) n k -1/η 0 1 + 1 η 0 δ 1 2 n k (15)
and by Lemma 1

c ω = ω 1 -ω d 2 g(1, ω 1-ω ) g(1, 1) = C ω g(1, 1) . ( 16
)
Thus using the convergence

√ m n m U m+1,n -1 Ω ∼ N (0, 1), (17)
we deduce that ln Z ω,n-m,n Z ω,n-m,n = ln the last equality coming from [START_REF] Alves | A new class of semi-parametric estimators of the second order parameter[END_REF].

1 + η 0 √ m W 1 cω -q 1 ( k n ) η 2 0 c ρ 0 ω ξω-1 ρ 0 + o P 1 √ m (1 -δ 1 2 n k (1 + o(1)))(1 -η 0 √ m Ω + o P ( 1 √ m ))(1 + δ ω (C η 0 ω U -η 0 m+1,n )(1 + o P (1))) = ln 1 + η 0 √ m W 1 cω -q 1 ( k n ) η 2 0 c ρ 0 ω ξω-1 ρ 0 + o P 1 √ m 1 -η 0 √ m Ω -
By an application of Taylor's theorem we finally get that

T 11 = 1 √ m W 1 c ω + 1 √ m Ω + o P 1 √ m .
We have now to study the second term, T 12 : (1 + o P (1)) -q 1 C -η 0 ω q η 0 k n (1 + o P (1))

T 12 =
= o P q 1 k n , by ( 15) and ( 16).

Thus

T 12 = o P 1 √ m .
Now the last term T 13 can be treated as follows

T 13 d = -ln n m F Zω (Q Zω (U n-m,n )) d = - n m U m+1,n -1 (1 + o P (1)) = - Ω √ m + o P 1 √ m
by [START_REF] Goegebeur | Kernel estimators for the second order parameter in extreme value statistics[END_REF]. Consequently

T 1 = 1 √ m W 1 c ω + o P 1 √ m . ( 18
)
Concerning T 2 . Using ( 5), ( 11) and ( 16

  variate context, and we refer here to Fraga Alves et al. (2003), Ciuperca and Mercadier (2010) and Goegebeur et al. (2010) for examples of estimators with a good practical performance. In the bivariate context the estimation of ρ (or τ ) is a challenging problem. Goegebeur and Guillou (2013) introduced a consistent estimator for τ , based on the stochastic process results of Draisma et al. (2004) for the empirical distribution function. However, this estimator does not seem to perform well in practice.

  Theorem 1 is an adaptation of the proof of Theorem 2 in Dutang et al. (2014) combined with Theorem A. Note that the asymptotic independence between the random threshold (first component of the (5×1)-vector in Theorem 1) and the other components has already been observed in the univariate framework, see for instance Lemma 5.3 in Dekkers et al.

Theorem 4

 4 extreme failure set. Let d n := m/(np n ). Under the conditions of Theorem 1, and assuming additionally np n = o(m) and ln np n = o( √ m) we have that

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Farlie Gumbel Morgenstern copula with ζ = 1: mean (left) of p n,α , MSE (right) of p n,α /p z,1 for n = 100, p z,1 = 2 n and different values of α: 0 (solid line), 0.1 (dashed line), 0.5 (dotted line) and 1 (dash-dotted line), based on 1024 simulations; from the top to the bottom: no contamination, contamination with ε = 2%, and ε = 5%. The parameter ρ is mis-specified to the value -1. The horizontal reference line in the left panels of the figure corresponds to the true value of p z,1 .

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Farlie Gumbel Morgenstern copula with ζ = 1: mean (left) of p n,α , MSE (right) of p n,α /p z,1 for n = 100, p z,1 = 1 10 n and different values of α: 0 (solid line), 0.1 (dashed line), 0.5 (dotted line) and 1 (dash-dotted line), based on 1024 simulations; from the top to the bottom: no contamination, contamination with ε = 2%, and ε = 5%. The parameter ρ is mis-specified to the value -1. The horizontal reference line in the left panels of the figure corresponds to the true value of p z,1 .

Table 2 :

 2 Values of the breakdown points, respectively for FGM ζ = 1, -1, and Frank (n = 100).

Table 3 :

 3 Values of the breakdown points, respectively for FGM ζ = 1, -1, and Frank (n = 1000).

Table 4 :

 4 FGM ζ = 1: coverage probabilities of 95% confidence intervals for p n = 0.01% based on Theorem 4 with n = 1000.

Figures 9 (

 9 Figures 9 (a) and (b) show the scatterplot of the original data and the data transformed into unit Pareto margins, respectively. The Pareto quantile plot of the Z 1 2 ,i , i = 1, ..., n, is depicted in Figure 9 (c). The three points labeled with crosses correspond to possible outliers in the dataset. Their values are given in Table5. Since X and Y are ratios (in terms of gross written premium), we may reasonably

Figure 7 :Table 5 :Figure 8 :

 758 Figure 7: Farlie Gumbel Morgenstern copula with ζ = 1: mean (left) of p n,α , MSE (right) of p n,α /p z, ω for n = 100, p z,1 = 2 n and different values of α: 0 (solid line), 0.1 (dashed line), 0.5 (dotted line) and 1 (dash-dotted line), based on 1024 simulations in the uncontaminated case. The first row corresponds to ω = 0.45 and the second row to ω = 0.55, where ω = (1 + ω) -1 . The parameter ρ is mis-specified to the value -1. The horizontal reference line in the left panels of the figure corresponds to the true value of p z, ω .

  Figure 10 shows the estimators p n,α as a function of m for different values of α: 0 (solid line), 0.1 (dashed line), 0.5 (dotted line) and 1 (dash-dotted line). In Figure10(a), we use the whole i ∈ N := {1, . . . , n} to estimate p n,α , whereas in Figure10(b) we use N \ {184, 101, 21}.

Figure 9 :Figure 10 :

 910 Figure 9: Loss/expense ratios dataset: (a) Scatterplot of the original data; (b) Scatterplot of the data transformed into unit Pareto margins; (c) Pareto quantile plot of the Z 1 2 ,i , i = 1, ..., n.

Figure 11 :

 11 Figure 11: Median of the estimator p n,α of P( X > z, Y > z) over m = 50, . . . , 80 for different values of α: 0 (solid line), 0.1 (dashed line), 0.5 (dotted line) and 1 (dash-dotted line). (a) Whole dataset; (b) Dataset without the three points.

2 in

 2 [START_REF] Dierckx | An asymptotically unbiased minimum density power divergence estimator for the Pareto-tail index[END_REF], we deduce that

Proof of Theorem 3 ZConcerning T 1 .

 31 We use the following decomposition ln n-m,n ; η n , δ ω,n , ρF Zω (z n ) = -ln n m F Zω ( Z ω,n-m,n ) + ω,n-m,n ; η 0 , δ ω,n , ρ 0 -ln F Zω zn Z ω,n-m,n Z ω,n-m,n F Zω ( Z ω,n-m,n )    + ln H z n Z ω,n-m,n ; η n , δ ω,n , ρ -ln H z n Z ω,n-m,n ; η 0 , δ ω,n , ρ 0 =: T 1 + T 2 + T 3 .We will study the three terms separately. Let Z ω,1,n ≤ . . . ≤ Z ω,n,n denote the order statistics of a random sample of size n from F Zω . This term can be rewritten asT 1 = -ln F Zω ( Z ω,n-m,n ) F Zω (Z ω,n-m,n ) -ln n m F Zω (Z ω,n-m,n ) = 1 η 0 ln Z ω,n-m,n Z ω,n-m,n -ln 1 + 1 η 0 δ ω,n 1 + 1 η 0 δ ω (Z ω,n-m,n )-ln n m F Zω (Z ω,n-m,n ) =: T 11 + T 12 + T 13 .

  m U m+1,n -η 0 (1 + δ ω (C η 0 ω U -η 0 m+1,n )(1 + o P (1)))

δ 1 2 nk ( 1 +.

 21 o(1)) + δ ω (C η 0 ω U -η 0 m+1,n )(1 + o P (1)) + o P ( 1 √ m )Combining Lemma 1 with the fact that |q 1 | is regularly varying with index τ and (3), we can infer thatδ ω (C η 0 ω U -η 0 m+1,n )) = -ξ ω

- 1 η 0 ( 1 +

 01 (δ ω,n (1 + o P (1)) -δ ω (Z ω,n-m,n )(1 + o P (1o P (1)) -δ ω C η 0 ω U -η 0 m+1,n (1 + o P (1

( 1 +F

 1 ), we have the two convergencesn m F Zω (z n ) = np n m -→ β, n m F Zω ( Z ω,n-m,n ) = n m C ω Z -1/η 0 ω,n-m,n (1 + o P (1)) = o P (1)) = 1 + o P (1), from which we deduce, by using Proposition 2.3 from Beirlant et al. (2009), thatT 2 = H zn Z ω,n-m,n ; η 0 , δ ω,n , ρ 0 -F Zω zn Z ω,n-m,n Z ω,n-m,n Zω ( Z ω,n-m,n ) F Zω (zn) F Zω ( Z ω,n-m,n ) (1 + o P (1)) = o P (δ ω,n ) = o P 1 √ m .

  and |δ ω | is a function of regular variation with index -τ . This second order condition is identical to the one used in the univariate framework in Beirlant et al. (2009) and Dierckx et al. (2013), and therefore, as shown in these papers, one can approximate the distribution of Z ω /u, conditionally on Z ω > u, for large u, by the extended Pareto distribution.

  1, η/ρ}). Using this property, one can estimate η by fitting h to the relative excesses over some large threshold u. Approximation[START_REF] Dekkers | A moment estimator for the index of an extreme-value distribution[END_REF] means that the difference between the two sides is uniformly of order o(|δ ω (u)|) (see Proposition 2.3 in[START_REF] Beirlant | Second-order refined peaks-over-threshold modelling for heavy-tailed distributions[END_REF] and thus this new extended Pareto model can be fitted to a larger portion of the data. This is in line with other approaches to bias-corrected estimation in the extreme value context. Indeed, in general, by taking the second order structure of extreme value models explicitly into account in the estimation stage, one obtains estimators with better bias properties than estimators based on first order models.

Table 1 :

 1 Value of z such that p z, ω

	45 10.486
	ω = 0.5	9.488	4.461	9.828 44.219 140.92
	ω = 0.55 8.579
	p z, ω		2%	0.1%	0.01%

Acknowledgement

This work was supported by a research grant (VKR023480) from VILLUM FONDEN and a PICS project (PICS-6416). The authors thank the editor and two referees for their helpful comments which led to significant improvements of this article.

Concerning T 3 . Remark that

Using the fact that

Combining ( 18)-( 20) with the fact that √ mδ ω,n P -→ -λη 2 0 ξ ω c ρ 0 ω /ρ 0 , our theorem follows.

Proof of Theorem 4 By straightforward calculations one obtains the following decomposition ln

For T 4 , note that

is in fact the well-known Weissman estimator [START_REF] Weissman | Estimation of parameters and larger quantiles based on the k largest observations[END_REF], and hence

From the proof of Theorem 3 we have that T


Piseth Heng

Mihai Bud

Hugues Somja

Mohammed Hjiaj
email: mohammed.hjiaj@insa-rennes.fr

Jean-Marc Battini

Residual Stiffness and Strength of Shear Connectors in Steel-Concrete Composite Beams after Being Submitted to a Pull-out Pre-damaging: An Experimental Investigation

Keywords: Explosion, shear stud, push-out, pullout, residual strength, composite floor, profile steel sheet

Horizontal stability of the medium rise steel frame structures is usually ensured by vertical bracings and diaphragm action of composite floors. Load transfer within the composite floor system is made through shear connectors, e.g. headed studs. In an event of explosion, such connectors must reserve sufficient residual stiffness and strength in order to avoid a sudden or delayed collapse of the building. These remaining capacities have not been experimentally studied yet in the literature. This paper presents large scale horizontal push out tests to determine the residual stiffness of the shear connectors after being initially damaged by explosion. The initial damaging is reproduced by a pull-out test using a quasi-static loading.

Two types of numerical simulation have also been developed using ABAQUS/CAE software to provide a better understanding of the experimental results.

Introduction

The present-day serviceability of steel frame buildings in industry and commerce has exposed the structures with high risk to extreme and accidental loadings, particularly explosion. Many incidents, such as Azote de France Fertilizer factory in Toulouse (AFT),

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Buncefield in UK and CPR Tank Farm in Puerto Rico to name a few, have damaged structures within a large surrounding area, causing serious casualties and substantial property loss. It is thus increasingly important to design structural buildings with such exposure to withstand the effects of blast loading. The Buncefield investigation led to the recommendations on land use planning and the control of societal risk around major hazard sites in addition to Health and Safety Executive (HSE) to take into account not only direct human harm but also property damage and other harms to communities. Coming after the incident of AFT in Toulouse, a regulation was established to allow the design of concrete buildings that can resist blast actions with an intensity up to 12 kPa. However, steel framed buildings are not permitted to tolerate even with low level of explosion except if the spacing between steel frames is less than 5 m. This consideration is based on inertia effects while robustness of steel structures, i.e. high ductility and redundancy, is not considered. Until recently, either large-scale buildings or concrete structures are included in the regulations while low-to-medium rise steel buildings subject to global explosion have not been studied yet. Consequently, a European project called Blast Actions on Structures In Steel (BASIS) has been proposed to investigate and develop a methodology for the safe design of low to medium rise office buildings subjected to large external explosions.

The research work in BASIS project is based on a reference building of a typical business park office building that fits the description of low to medium rise building. This threestorey office structure (See Fig. 1 and2) was designed based on usual loadings in simple construction using the principle of Eurocodes 1993 and 1994. In this type of frame building, the composite floors are supported by steel columns. The wind actions are transmitted from façades to the vertical steel bracings by the diaphragm action of the composite floor and by shear actions of shear connection. From the reference building design in [START_REF] Kyprianou | Reference Building Design. Blast Actions on Structures in Steel (BASIS), RFSR-CT-2013-00020[END_REF], the wind action expected to transmit through the shear connectors was estimated at 10 kN/m. The response of such multi-storey steel frame building subjected to blast loading involves complex interactions between building envelope, composite floor, frame system, and connections. The load transfer between the composite floor and its frames depends strongly on the performance of the shear connectors. Headed steel stud anchors welded to steel beam and encased in IPE 360 [START_REF] Hicks | Strength and Ductility of Headed Stud Connectors Welded in Modern Profiled Steel Sheeting[END_REF] IPE 360 [START_REF] Hicks | Strength and Ductility of Headed Stud Connectors Welded in Modern Profiled Steel Sheeting[END_REF] IPE 300 [START_REF] Hicks | Strength and Ductility of Headed Stud Connectors Welded in Modern Profiled Steel Sheeting[END_REF] IPE 270 [START_REF] Hanswille | Resistance of headed studs subjected to fatigue loading Part II: Analytical study[END_REF] IPE 200 [START_REF]Part 3.1: Structural Use of Steelwork in Building[END_REF] IPE 200 [START_REF]Part 3.1: Structural Use of Steelwork in Building[END_REF] IPE 270 [START_REF] Pallarés | Headed Steel Stud Anchors in Composite Structures, Part II: Tension and Interaction[END_REF] / empirical studies using small-scale standard or full scale push-out tests. After the earlier push-out tests [START_REF] Viest | Investigation of Stud Shear Connectors for Composite Concrete and Steel T -Beams[END_REF][START_REF] Chinn | Pushout Tests on Lightweight Composite Slabs[END_REF], small-scale standard tests were established in CP117 [START_REF]Part 1: Composite Construction in Structural Steel and Concrete[END_REF], BS5950 [START_REF]Part 3.1: Structural Use of Steelwork in Building[END_REF] and Eurocode 4 [START_REF] En | Eurocode 4: design of composite steel and concrete structures-Part 1.1: General rules and rules for buildings[END_REF]. In such tests, the push-out specimen is composed of a profiled steel beam connected to two blocks of concrete by shear connectors. Ollgaard et al in 1971 [START_REF] Ollgaard | Shear Strength of Stud Connectors in Lightweight and Normal Weight Concrete[END_REF] used the small scale push-out test setup to investigate the capacity and behavior of stud shear connectors embedded in lightweight and normal weight concrete. Prakash et al [START_REF] Prakash | Modified Push-out Tests for Determining Shear Strength and Stiffness of HSS Stud Connector-Experimental Study[END_REF] presented modified pushout tests to determine shear strength and stiffness of high strength steel studs, in which the reinforcement, method of casting concrete and test setup were altered from the standard tests with an adaptation to enhance splitting resistance of reinforced concrete. In addition to this, Shim et al [START_REF] Shim | Static Behavior of Large Stud Shear Connectors[END_REF] performed pushout tests on specimens fabricated in accordance with Eurocode 4 [START_REF] En | Eurocode 4: design of composite steel and concrete structures-Part 1.1: General rules and rules for buildings[END_REF] to determine static behavior of large stud shear connectors (diameter 25, 27, 30 mm stud). On the other hand, Lloyd and Wright [START_REF] Lloyd | Shear Connection between Composite Slabs and Steel Beams[END_REF] conducted a series of large scale push-out tests of composite beam specimens (See Fig. 3) with the objective to recommend a standard configuration for through-deck push-out test. In their study, 42 push-out tests have been performed by varying some parameters such as the width of the slab, the detail of mesh reinforcement, and the number of profiled steel sheet pitches. Having similar test setup to the one described in [START_REF] Lloyd | Shear Connection between Composite Slabs and Steel Beams[END_REF], Smith and Couchman [START_REF] Smith | Strength and Ductility of Headed Stud Shear Connectors in Profiled Steel Sheeting[END_REF] carried out 27 push tests on specimens of composite beam using a new push rig to investigate the effect of variables such as meshed position, transverse spacing of shear connectors, number of shear connectors per trough and the slab depth. In 2009, Hicks [START_REF] Hicks | Strength and Ductility of Headed Stud Connectors Welded in Modern Profiled Steel Sheeting[END_REF] also implemented two types of tests on the headed stud connectors welded in profiled steel sheeting and embedded in reinforced concrete slab. One of his setups is a standard push-out test as defined in Eurocode 4 [START_REF] En | Eurocode 4: design of composite steel and concrete structures-Part 1.1: General rules and rules for buildings[END_REF] whereas the other setup is a full-scale composite beam illustrated in Fig. 4. Apart from this, Lam [START_REF] Lam | Capacities of Headed Stud Shear Connectors in Composite Steel Beams with Precast Hollowcore Slabs[END_REF] proposed a new horizontal push-out test setup using full-scale specimens of composite beams with

A C C E P T E D M A N U S C R I P T

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

precast hollowcore slab. He performed 7 push-out tests with a new configuration of headed studs in solid RC slabs to validate the test procedure, and the results agreed well with the ones specified in the available codes [START_REF]Part 3.1: Structural Use of Steelwork in Building[END_REF][START_REF] En | Eurocode 4: design of composite steel and concrete structures-Part 1.1: General rules and rules for buildings[END_REF]. This new full-scale horizontal test arrangement is proposed due to the reasons that the standard push-out test specified in Eurocode 4 does not fit for a specimen of precast hollowcore composite slab. Later in 2011, Qureshi et al [START_REF] Qureshi | Effect of Shear Connector Spacing and Layout on the Shear Connector Capacity in Composite Beams[END_REF] adopted the new horizontal push-out test configuration presented by Lam [START_REF] Lam | Capacities of Headed Stud Shear Connectors in Composite Steel Beams with Precast Hollowcore Slabs[END_REF] to carry the push tests on headed shear studs in composite beams with composite profiled sheeting slab.

Parametric studies were then performed through a simulation of the push-out test using a 3-dimentional nonlinear finite element model in ABAQUS/Explicit [START_REF]ABAQUS documentation[END_REF] to investigate the effect of shear connector spacing and layout on the shear connector capacity.

ACCEPTED MANUSCRIPT

in performing diaphragm action both to transmit lateral loads to the bracing system and to pass on vertical load to frame system. In such a case, the headed shear studs are subjected to combined tension and shear action. Pallarés and Hajjar [START_REF] Pallarés | Headed Steel Stud Anchors in Composite Structures, Part II: Tension and Interaction[END_REF] assessed the limit state formulas for headed stud anchors under shear-tension interaction action against 54 experiments found in literatures [START_REF] Saari | Behavior of Shear Studs in Steel Frames with Reinforced Concrete Infill Walls[END_REF][START_REF] Mcmakin | Headed Steel Anchor under Combined Loading[END_REF] in the context of AISC specifications and compared the results with standard provisions [START_REF]AISC) Specification for Structural Steel Building American Institute for Steel Construction[END_REF][START_REF]American Concrete Institute Committee 318 (ACI) Building Code Requirements for Structural Concrete (ACI 318-08) and Complementary (ACI 318R-08)[END_REF][START_REF]Pre-stressed Concrete Institute (PCI) PCI Design Handbook[END_REF]. However, the structural behavior of headed shear studs embedded in composite slabs with profiled steel sheeting under combined action has not

been assessed yet. Besides, to the best knowledge of the authors, the influence of the initial damage of the composite beam-slab connection caused by uplift blast action onto shear capacity of headed shear stud connectors has not been investigated in the literature yet.

However, the knowledge of this residual strength is important for the design since the load bearing capacities of the composite beam rely on the composite action between the steel beam and the composite slab.

This paper aims at evaluating the influence of explosion on the shear capacity of headed shear stud connectors welded to steel beam and embedded in composite floor with profiled steel sheeting. A large-scale horizontal push-out test of the composite beam is employed to determine both initial and residual shear properties of the headed shear connectors. The initial damage caused by the explosion is experimentally reproduced by a pull-out test of the composite beam specimen. Two types of numerical simulations of the composite beam specimen, namely a simple beam model and a full solid element, have been created in ABAQUS/CAE [START_REF] Abaqus | Abaqus Documentation Version 6.13[END_REF] to validate and analyze the phenomenon in the experimental test.

Experimental Plan

To be able to determine experimentally the residual stiffness of the shear connectors in the composite beam after it has suffered out-of-plan damage caused by high pressure from blast, the following experimental plan was defined. First, three composite beam specimens were submitted to push-out tests at Institut National des Sciences Appliquées de Rennes (INSA de Rennes) in order to find their initial shear characteristics. Next, two of the specimens were sent to Institut National de l'Environnement Industriel et des Risques (INERIS) for

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

the explosion test to generate damage due to blast. The third specimen was prepared to serve as a reference to define the shear capacities of the connectors without initial blast damage. After being exposed to explosion, the two specimen were delivered back to INSA de Rennes to perform the push-out tests in order to determine the residual shear stiffness and the ultimate shear resistance of the specimen. At this stage, the third specimen was also tested for the shear resistance.

Unfortunately, the explosion tests at INERIS had completely destroyed the two composite specimens, which made it impossible to test them for residual stiffness. As a consequence, it was decided to alter the experimental plan. A supplementary pull-out test was performed on the third specimen to quasi-statically reproduce the pre-damaging of the specimen. It has been proven in [START_REF]Final report of Blast Action on Structures in Steel (BASIS)[END_REF] that the response of the specimen under explosion in the experiment was quasi-static. Such behaviour can be reproduced by the pull-out type test. The damaged specimen was then tested for the residual stiffness and the ultimate shear resistance.

Specimen and Material Properties

The composite beam specimen was defined with regard to the aforementioned reference building, with some necessary adaptations following the installation capacities of the labo- 1 where the yield strength, the mean value of yield strength, the ultimate strength and the mean value of the ultimate strength are denoted by f y , f ym , f u and f um , respectively. For all the materials, the obtained mean values passed the nominal yield strength. On the other hand, the concrete material

ratory

Material

Coupon

f y [MPa] f ym [MPa] f u [MPa] f um [MPa]
Comfor of concrete specimens were carried out at the age of 33 days, 99 days and 121 days. The tests at the age of 33 days were to determine the nominal strength of the concrete while the tests at 99 days and 121 days were made before and after the initial stiffness tests respectively in order to know the concrete strength at the approximated time of testing. The results of the cylinder compressive tests are summarized in Table 2 where the concrete strength, the mean value of concrete strengths are denoted by f c and f cm , respectively. Among the cylinder tests, two tests were not taken into consideration due to erroneous results. The concrete class specified in the design of reference building was C30/37; the mean value of the concrete strength from the tests was obtained at f cm = 37.15 MPa, which is quite comparable to the mean value theoretically assumed for the concrete C30/37 (f cm = 38 MPa).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Push-out Test: Experimental Setup and Loading Procedure

The experimental setup of the horizontal push-out test of the composite floor was positioned on a rigid frame system (meshed floor and reacting frame), as shown in Fig. 6. The rigid frame system was designed to have large rigidity so that its displacements in all directions are not allowed to occur. Local buckling of the elements in the rigid system is also avoided by stiffeners. Having one end fixed to the reacting frame by bolts, a hydraulic in which

P Rd,1 = 0.8 × f u × π d 2 4 γ v = 0.8 × 450 × π 19 2 4 γ v = 81.7 kN (2)
P Rd,2 = 0.29 × a × d 2 × √ f ck × E c γ v × 1000 = 0.29 × 1 × 19 2 × √ 35 × 30 × 10 3 1.25 × 1000 = 85.8kN (3)
where f u is nominal yield strength of the stud; γ v , safety factor; d, diameter of the stud; f ck , nominal strength of concrete; E c , Young modulus of concrete; a = 1.

However, the shear resistance of the stud in the sheeting profile, which is the design shear resistance per stud, is determined by multiplying the shear resistance per stud in solid concrete panel P Rd,c with a reduction factor k t as

P Rd,s = k t × P Rd,c = 81.7 × 1.04 = 85 kN (4)
in which

k t = 0.7 b 0 n r h s h stud h s -1 = 0.7 × 105 √ 1 × 46 76 46 -1 = 1.04 (5)
with n r , number of stud per trough; b 0 , average width of trough; h s , height of the sheeting profile.

The estimated shear resistance per stud is then calculated by multiplying the safety

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

factor with the design shear resistance as following

P Rd = γ v × P Rd,s = 1.25 × 85 = 106.2 kN (6)
For 8 studs, the ultimate shear resistance was approximately at 850 kN. The capacity of the force jack of 1500 kN and the design of the bracing system and test setup for the load of 1500 kN were sufficiently strong to ensure proper test conditions.

In the push-out test for initial and residual stiffness, the maximum load was limited at 100 kN to avoid possible cracking in concrete material. Three phases of the servo-controlled loading procedure were exerted:

• The load was increased until 10 kN;

• 10 cycles of loading and unloading procedure between 10 and 20 kN were applied;

• 5 cycles of loading and unloading procedure between 10 and 100 kN were applied.

In push-out test for residual ultimate resistance, the load was increasingly applied until the collapse of the specimen.

The longitudinal slip was measured by LVDT sensors with a maximum displacement capacity of +/-2.5 mm, placed at ten locations as indicated in Fig. 8, for the initial and residual stiffness push-out test. On the other hand, for residual ultimate resistance push-out test, the longitudinal slip was recorded only at four locations where LVDT sensors with a maximum displacement capacity of +/-25 mm were located (see Fig. 9). Fig. 10 gives a general view of specimen and test setup. the maximum applied load that would not completely damage the specimen, the loading procedure follows two phases:

• First, the specimen was loaded up to 70 kN, and the initial stiffness is determined based on force-mean displacement curve.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

• Then, the loading was resumed with low rate so that it could be stopped if the displacement lead to a secant stiffness that was ten time smaller than the computed initial stiffness or if a significant change of the evolution of the specimen could be visually observed.

Pull-out Test: Experimental Results

The pull-out test of the third specimen was carried out on 28/04/2015. Two force actuators applied downward forces on the specimen by imposing displacements, and the force history was recorded by the Deltec force sensors. The relative displacement between the two force jacks was adjusted along the test in order to maintain similar forces in each jack. Two cycles of loading and unloading between a total force of 30 and 70 kN were made to acquire information about the initial stiffness. The force-displacement curve using global displacement is plotted in Fig. 13. The loading was stopped as soon as the force was seen decreased with increasing displacement. The maximum force attained was measured to 162.6 kN, which gives approximately 20 kN per stud. The evolution of the uplift displacements measured by the LVDT sensors is presented in Fig. 15. The uplift displacement grew larger

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

at the extremities. It was considered that the studs had reached their maximum limit load with the evidence of the uplift displacement ranging between 0.6 mm to 2.3 mm. It was concluded that the yielding of rebars happened as a plateau was observed in the forcedisplacement curve. The cracking of concrete after pull-out test is shown in Fig. 14. The concrete zones around the shear connectors were surely damaged because large values of beam-concrete uplifts were measured. Yielding of steel sheet was also noticed. To make a verification of the result with the analytical expression [START_REF] Hanswille | Resistance of headed studs subjected to fatigue loading Part II: Analytical study[END_REF], the force at each stud is calculated using the expression in Eq. ( 7). The force summed for all the studs are presented in Table 3, and the level of magnitude of the force agrees well with the actual

A C C E P T E D M A N U S C R I P T

Table 3: Initial stiffness push-out test: total force of all the studs from slip measurement Specimen 1 2 3

Total force [kN] 106.5 73.5 90 applied load of 100 kN.

P = P u 1 -e -1.22 S 0.59 (7)
where P u = 106.2 kN is the ultimate resistance of the stud (see Eq. ( 6)) and S is the slip of the stud.

Residual stiffness test

After pull-out test, the third specimen was submitted to push-out test for residual stiffness on 07/05/2015. For this test, the contact between the concrete and the support was limited to a width of 60 cm only. The evolution of slips with the force is shown in Fig. 18.

The response behavior in the first cycle of the loading and unloading between 10 and 100 kN differed from the other four cycles. This can be explained by the fact that the concrete zone around the stud must have been damaged by the pull-out test, which required some movement in the first cycle to obtain good contact again. The comparison is also made in

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Residual ultimate resistance test

To acquire the residual resistance of the composite beam specimen after pre-damaging by pull-out force, the push-out test of the specimen up to the ruin of the specimen was performed on 13/05/2015. The evolution of the mean slip against the force is presented in Fig. 21. There was an accidental move of the test setup, which made an unloading at around 400 kN. However, the reloading reached the same envelop curve. The maximum load corresponding to the residual resistance of the specimen was attained at 530 kN. Although the initial resistance was not determined from experimental test, the formula provided in EN 1994.1.1/6.6.3.1 gives an initial resistance of 850 kN. From that, it can be deduced that the resistance of the composite beam is reduced by around 40 percent for a pre-damaging that has led to the yielding of the rebars under pull-out action with a mean uplift of 1.3 mm. In this section, a full model of the composite beam specimen is presented. Created in ABAQUS/Explicit, the concrete, shear stud and IPE-200 steel beam were modelled using a linear 8-node reduced integration hexahedral continuum element (C3D8R) whereas the profiled steel sheet was modelled by a linear 4-node reduced integration quadrilateral shell element (S4R). The meshed rebars were, on the other hand, modelled by a linear 2-node 3D truss element (T3D2). Each of the material parts was partitioned in a way that structured meshing technique could be applied. The concrete floor part was cut out at the zone where the studs were located. The contact interaction with general contact was used to model the A friction coefficient of 0.5 was used for the contacts between the profiled steel sheet and the concrete panel and between the shear studs and the concrete panel, and 0.1 for the other contacts. Tie constraints were used to tie shear studs to the IPE-200 steel beam.

The rebar mesh was embedded in the concrete panel. The HEB-300 steel beams that support the specimen are constrained to a reference point using Rigid-Body Constraint.

The configuration of the modelling used to simulate the pull-out test is shown in Fig. 23.

The boundary conditions for the pull-out test and the push-out test are illustrated in Fig. 24(a) and 24(b), respectively. In the pull-out test, the reference points that constrain the HEB-300 steel beams are fixed in all directions except the rotation around x-axis. In the push-out test, the picked bottom surfaces of the IPE-200 steel beam are fixed in vertical direction (y-axis), the picked edge surface of the concrete panel is fixed in horizontal direction (x-axis), the picked top surface of the concrete panel is fixed in vertical direction (y-axis), and the reference points that constrain the HEB-300 steel beams are fixed in all directions except the rotation around x-axis.

A C C E P T E D M A N U S C R I P T

Material model

The concrete damaged plasticity model was used for concrete material behavior. For the behavior in compression, the stress-strain relationship was used in accordance to EN 1992.1.1:

σ c = f cm kη -η 2 1 + (k -2)η (8)
in which

k = 1.05 E cm ε c1 f cm (9) η = ε c ε c1 (10)
where ε c1 = 0.7 f 0.31 cm is the strain at ultimate compressive strength of concrete; ε c is the compressive strain in concrete; E cm is the secant modulus of elasticity. Following [START_REF] Qureshi | Effect of Shear Connector Spacing and Layout on the Shear Connector Capacity in Composite Beams[END_REF], the tensile behavior of the concrete was modelled by the tensile stress (σ t ) versus the cracking displacement (w) relationship that was defined by:

σ t f t = f (w) - w w c f (w c ) (11)
with

f (w) = 1 + c 1 w w c 3 exp - c 2 w w c (12)
w c = 5.14 G f f t (13)
where f t is the tentile strength, c 1 = 3 and c 2 = 6.93, for normal density concrete, are material constants, and the fracture energy G f is defined by [START_REF] Richard | Guidelines for FE analysis of concrete dams[END_REF]:

G f = 73f 0.18 cm (14)
For steel material model of profiled steel sheet, shear studs and rebar, an elastic perfectly plastic material behaviour is used. The modulus of elasticity of all the steel steel components is taken as E = 200 GPa. The yield stress of the profiled steel sheet, shear studs and profile The push-out test was also simulated in a 3D model. The model was pre-damaged by the simulation of the pull-out test, and then submitted to the push-out load in order to determine the residual stiffness and resistance. The finite element modelling of the postfailure behaviour of this composite beam specimen under pull-out and push-out actions is highly complicated in account of the inelastic brittle material behaviour of concrete, the localized characteristics of shear stud connectors and the simultaneous contact interaction between various components. Various uncertainties such as contact property between the concrete panel and steel sheet also take part in the difficulty of the numerical simulation.

Due to such difficulty to model the coupling effect of the damage caused by pull-out force to the shear resistance, the simulation was not able to reproduce the shear stiffness of the shear connectors in the composite specimen obtained in the experimental experiment.

In spite of that, the curves of the push-out load versus the mean slip are plotted in Fig. 32, illustrating the comparison of the results obtained from the experimental test and the FEM simulation. Considering the mentioned difficulties, it can be observed that the load-slip behaviour in the numerical simulation (with an initial global pull-out displacement 18.75 mm) resembles the one obtained in the experimental test reasonably well.

The applied push-out force is distributed to the shear studs in a way that the stud closer to the applied force is under the most stress (Fig. 33). This can also be illustrated by the equivalent plastic strain in Fig. 34.

The load-slip curves are also illustrated in Fig. 35

Conclusion

This paper presents a series of large-scale push-out tests and pull-out tests carried out with the aim to evaluate the influence of pre-damaging on the shear capacity of headed shear stud connectors welded to steel beam and embedded in composite floor with profiled steel sheeting. Three specimens were submitted to push-out tests to determine the initial shear stiffness of the shear connectors. Next, a pull-out test was carried out on one specimen to reproduce the pre-damaging by explosion of the specimen. The pre-damaged specimen was then submitted to shear load in another push-out test to find the residual shear stiffness.

Finally, the residual resistance of the specimen in shear was determined in a last push-out test that applied the shear load until full collapse of the specimen.

Two models of numerical simulation have also been performed to assess a better understanding of the experimental tests. The results from the numerical simulation provide the following insights. For one thing, experimental results of initial push-out tests are consistent with the information provided by the large historical background about the headed stud behavior under shear loading. Second, the yielding of the specimen attained in the pull-out test was confirmed and quantified.

At the end, the following information can be concluded:

• For studs in a steel-concrete connection of a composite beam that are submitted to a pull-out action that leads to the damage of the concrete around the shear studs, the residual stiffness in shear is reduced by a factor in a range between 10 and 20 compared to their initial stiffness.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

• Their residual resistance in shear is decreased by a factor of about 2 compared to their initial theoretical resistance without pre-damaging.

This preliminary study gives only partial information. However, from a practical point of view, the order of magnitude for the reduction of the initial stiffness due to pre-damaging can be considered accurate enough, bearing in mind the large uncertainties on the level of the blast actions. Regarding the resistance, it can also be inferred that the residual strength of the shear studs is sufficiently large to resist the usual wind actions as this latter induces forces that are small compared to the residual resistance.

As a conclusion, this paper gives a first full practical approach for evaluating the residual capacities of shear stud connection after being submitted to a blast action.

Figure 1 :

 1 Figure 1: Architectural façade view

Figure 2 :

 2 Figure 2: Level 1 and 2, office areas

Figure 3 :Figure 4 :

 34 Figure 3: Lloyd and Wright's push test arrangement series A

  at INERIS. Illustrated in Fig.5, the specimen was composed of a composite floor with steel profile sheet, an IPE-200 steel beam made of S355 steel and a line of 8 shear studs with a nominal yield strength of 450 MPa. The composite floor, with the dimension of 1800mm-by-2000mm and the thickness of 100 mm, was fabricated by casting concrete of a C25/30 type on a corrugated steel sheeting Comflor 46 whose nominal yield limit was 280 MPa (supplied by Tata Steel) with HA-6mm meshed reinforcement. A 20-mm concrete cover was chosen in order to have the same amount of concrete over and under the rebar mesh above the crest of the steel sheet. The spacing of reinforcement was arranged to be 225 mm in both directions. In each flute, an HA-8mm rebar was placed as prescribed in the technical documentation of the Comflor. As the flute at the edges were changed to respect the width of the tunnel used to apply the blast action, a supplementary reinforcement was added to obtain the same steel-concrete ratio. In the other direction, the reinforcement wasA C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPTadapted at the edges because the length of the specimen was also fixed by the dimension of the blast setup. Consequently, for edges perpendicular to crest lines, two HA-6mm rebars with spacing of 75 mm were welded to 6-mm rebar mesh, see Fig. 5. For edges parallel to crest lines, the rebar frame in the edge flutes consisted of four HA-6mm rebars placed along the edges. The lower rebars were fixed to the rebar mesh by stirrup of 6-mm diameter with 75-mm spacing. Similarly, in each flute, two HA-6mm rebars were welded to both the lower 8-mm rebar along the crest line and the two HA-6mm rebars. The rebars of diameter 6 mm and 8 mm had a nominal yield strength of 500 MPa. In order to know the actual yield strength, six tensile coupon tests were made for each material, i.e. Comflor 46, 8-mm rebar, 6-mm rebar and shear stud. The results are summarized in Table

Figure 5 :

 5 Figure 5: Composite beam specimen

  A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPTservo controlled actuator applied a horizontal load on one of the IPE 200 beam's end crosssections. This force jack has a load capacity of 500 kN for the initial and residual stiffness test and 1500 kN for the residual ultimate resistance test. At the other end of the specimen, the composite floor was placed in contact with a rigid steel beam all along its edge's width (See detail A and C in Fig.6). The uplift displacement of the floor at this edge was prevented by an L-shape steel profile (See detail A). The composite beam specimen was placed horizontally on two vertical supporting steel columns called potelet. Each interface between the IPE-200 steel beam and the potelet was equipped with double greased PTFE plates to reduce friction between the surfaces. The poletets are fixed to the rigid meshed floor by bolts. The uplift displacement of the steel beam was also avoided by bolts in slotted holes (only horizontal movement is allowed) as shown in Fig. 7. According to the formula

Figure 6 :Figure 7 :

 67 Figure 6: Push-out test setup

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: LVDT sensors: initial and residual stiffness test

Figure 11 :

 11 Figure 11: Pull-out test set-up

Figure 12 :

 12 Figure 12: (a). Potentiometer sensors (b). LVDT sensors.

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Pull-out test: force-displacement curve using global displacement

Figure 16 :Figure 17 :

 1617 Figure 16: Initial stiffness push-oute test: Slips at the applied load of 100 kN

Fig. 19 betweenFigure 18 :

 1918 Fig.19between the mean slips obtained in the initial stiffness push-out test and the slip observed in the residual stiffness test for the first and the last loading/unloading cycle. The mean value among the slips are summarized in Table4. It can be approximately deduced from the results given in Table4that the residual stiffness of the composite beam equals to the initial stiffness divided by a factor in a range between 10 and 20. The force-mean slip

Figure 19 :Specimen 1 (

 191 Figure 19: Push-out test: comparison of slips at 100 kN

Figure 20 :

 20 Figure 20: Push-out test: comparison of force-slip curves for the last cycle

Figure 21 : 8 . 1 .

 2181 Figure 21: Residual resistance push-out test: force-slip curve

Figure 22 :

 22 Figure 22: FE simulation: simplified model

Figure 23 :Figure 24 :

 2324 Figure 23: FE simulation: Full solid model

Figure 26 :Figure 28 :Figure 30 : 8 . 4 . 3 .

 262830843 Figure 26: Initial stiffness test: Comparison of slips between experimental data and numerical simulations

Figure 32 : 4 (b) Stud 5 to 6 Figure 33 :Figure 34 :

 32463334 Figure 32: Comparison of residual ultimate load between experimental test and FEM simulation

Table 1 :

 1 Results of tensile coupon tests.

Table 2 :

 2 Results of cylinder compressive testsTest date Age[days] f c [MPa] f cm [MPa]

			37.23
	01/09/2014	33	36.26	37.15
			37.96
			42.98
	06/11/2014	99	13.47	43.40
			43.82
			41.87
	28/11/2014	121	42.63	42.25
			16.24

Table 4 :

 4 Residual stiffness push-out test: mean slip of the studs at 100 kN

	Test	Mean slip [mm] Ratio Residual/Initial stiffness
	Initial stiffness test	0.022	1
	Residual stiffness test -last cycle	0.23	10.4
	Residual stiffness test -first cycle	0.38	17.4

Acknowledgements

The authors gratefully acknowledge financial support by the European Commission (Research Fund for Coal and Steel) through the project BASIS under grant agreement RFSR-CT-2013-00020.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

beams is assumed to be 355 MPa while yield stress of the rebars is taken as 535 MPa.

Parametric study

In order to better study the effect of pre-damaging on the resistance and stiffness of the shear stud connectors, the level of the pre-damaging is varied in the numerical simulation by changing the global pull-down displacement, which was 18.75 mm in the experimental test (see Fig. 13). The values considered are 0, 10 mm, 18.75 mm, and 25 mm. The first value will provide information about the initial resistance that is lacking in the experimental test.