Véronique Maume-Deschamps
email: veronique.maume@univ-lyon1.fr

Didier Rullière
email: didier.rulliere@univ-lyon1.fr

Antoine Usseglio-Carleve

D Rulli Ère

Spatial Quantile Predictions for Elliptical Random Fields

Keywords: Elliptical distribution, Quantile regression, Extremal quantile, Spatial prediction, Kriging

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Kriging, introduced by [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF], and formalized by [START_REF] Matheron | Traité de géostatistique appliquée[END_REF], aims at predicting the conditional mean of a random field (Z t ) t∈T given the values Z t1 , ..., Z t N of the field at some points t 1 , ..., t N ∈ T , where typically T ⊂ R d . When using the Kriging techniques, for any x ∈ T , the conditional mean of Z x given Z t1 , ..., Z t N is approximated by a linear combination of Z t1 , ..., Z t N where the weight vector is the solution of a least square minimization problem (see Ligas and Kulczycki, 2010, for example). It seems natural to predict, in the same spirit as Kriging, other functionals by linear combinations. Our starting point is to apply this method to conditional quantiles and then get spatial quantile predictions. In 1978, Koenker and Bassett proposed a conditional quantile estimation as an affine combination of Z t1 , ..., Z t N , called Quantile Regression (cf. [START_REF] Koenker | Regression quantiles[END_REF]. The weight vector is the solution of a minimization problem, with an asymmetric loss function, different from the least squares. Obviously, for Gaussian fields, Z x |(Z t1 , ..., Z t N ) is still Gaussian, and the conditional quantile of Z x given (Z t1 , ..., Z t N ) is easily computed. For non Gaussian fields, explicit formulas for conditional quantiles are more difficult to get. The Quantile Regression approach usually requires time consuming simulations to compute expectations. Moreover, in a non-gaussian setting, we do not have the certainty that the conditional quantile is expressed as a linear combination of the covariates, thus the consistency of the estimation by quantile regression is not guaranteed. In this paper, we focus on elliptical random fields. Elliptical distributions, formalized by [START_REF] Cambanis | On the theory of elliptically contoured distributions[END_REF], have the advantage of being stable under affine transformations. Therefore, explicit formulas for the quantile regression may be obtained for consistent elliptical distributions (cf. Kano (1994)). Nevertheless, the quantile regression is generally not equal to the conditional quantile and the difference may be large, especially for extreme levels of quantile. This is why we propose a spatial quantile prediction that is adapted to extremal quantiles.

The paper is organized as follows. In Section 2, we give some definitions, properties and examples of elliptical distributions satisfying the consistency property. For these models, we give formulas for conditional quantiles in Section 3. The Section 4 is devoted to quantile regression for consistently elliptical random fields: closed formulas are obtained. In Section 5, we propose some extremal predictions and prove asymptotic equivalences when the quantile level is close to 0 or 1. Section 6 provides a numerical study. In particular, we emphasize the fact that quantile regression is generally not consistent, especially for high level quantiles. We illustrate this point on several examples. This work was supported by the LABEX MILYON (ANR-10-LABX-0070) of Universit à c de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR)..

Elliptical distributions

In this section, we recall some useful properties and classical examples of elliptical distributions. Most results may be found, for instance, in [START_REF] Frahm | Generalized Elliptical Distributions: Theory and Applications[END_REF].

Definition 2.1. Let X be a d-dimensional random vector. X is elliptical if and only if there exists a unique µ ∈ R d , a semi-positive definite matrix Σ ∈ R d×d , and a function Φ : R + → R such that the characteristic function of (X -µ) is E [exp(it(X -µ))] = Φ(t Σt) . For such an elliptical random vector, we write X ∼ E d (µ, Σ, Φ)

It seems important to note that Σ is not necessarily the covariance matrix of X. More precisely, Σ is proportional to the covariance matrix K of the random vector X, when it is defined, i.e there exists a positive coefficient τ such that:

(2.1) Σ = τ K For example, for Gaussian distributions, τ = 1, i.e Σ = K. But this is not always the case for all elliptical distributions: for Student distributions with ν degrees of freedom, τ = ν-2 ν . Furthermore, K may not exist (e.g. for Cauchy distributions). In the present paper, we do only consider the case of non-degenerated distributions, i.e. we assume that the matrix Σ is invertible.

The following representation theorem is central in the theory of elliptical distributions. This result may be found in [START_REF] Cambanis | On the theory of elliptically contoured distributions[END_REF].

Theorem 2.1 [START_REF] Cambanis | On the theory of elliptically contoured distributions[END_REF]). The random vector X is elliptical, X ∼ E d (µ, Σ, Φ), if and only if (2.2) X = µ + RΛU (d) ,

where ΛΛ T = Σ, U (d) is a d-dimensional random vector uniformly distributed on S d-1 (the unit disk of dimension d), and R is a non-negative random variable independent of U (d) .

The representation of Theorem 2.1 is not unique (see Cambanis et al., 1981, for details). Given µ and Σ, the elliptic random vector X is characterized by the non-negative random variable R, called the radius of X. We now recall the consistency property of an elliptical distribution. The related definitions and properties may be found in Kano (1994). Kano (1994) established the following relation between Definition 2.2 and the radius R.

Definition 2.2. Let X ∼ E d (µ, Σ, Φ). X is said consistent if Φ is dimension-free, i.e. if Φ does not depend on d.
Proposition 2.2. Let X ∼ E d (µ, Σ, Φ), and let R be the corresponding radius of X. X is consistent if and only if :

(2.3) R d = χ d ,
where χ d is the square root of a χ 2 distributed random variable with d degrees of freedom, is a nonnegative random variable which does not depend on d, and χ d , and U (d) are mutually independent.

Theorem 2.3 (Elliptical density). Let X ∈ R d be an elliptical random vector, X ∼ E d (µ, Σ, Φ), and let R be the corresponding radius of X.

(2.4)

f X (x) = c d | det(Λ)| g d (x -µ) T Σ -1 (x -µ) where c d g d (t) = Γ( d 2 ) 2π d 2 √ t -(d-1) f R ( √ t), and f R (t) is the p.d.f of R.
The coefficient c d is called the normalization constant and the function g d is called the generator of X. Table 1 provides some examples of elliptical distributions, associated with their coefficients and generators. Most of them may easily be found in the literature: Kotz distribution is introduced in [START_REF] Nadarajah | The kotz type distribution with applications[END_REF], Student in [START_REF] Nadarajah | Multivariate T-Distributions and Their Applications[END_REF], and Laplace in [START_REF] Eltoft | On the multivariate Laplace distribution[END_REF] and [START_REF] Kozubowski | Multivariate generalized Laplace distribution and related random fields[END_REF]. We added two other distributions, obtained with Kano's representation : Unimodal Gaussian Mixture (Unimodal GM), and another one which will be called Uniform Gaussian Mixture (Uniform GM).

In Table 1, K m (x) denotes the modified Bessel function of the second kind with order m, and χ 2 m (x) denotes the c.d.f of the χ 2 distribution with m degrees of freedom, evaluated at x. Remark that the

Distribution

Coefficient

c d Generator g d (t) Gaussian 1 (2π) d 2 exp(-t 2 ) Student, ν > 0 Γ( d+ν 2 ) Γ( ν 2 ) 1 (νπ) d 2 1 + t ν -d+ν 2 Logistic Γ( d 2 ) (2π) d 2 +∞ 0 x d 2 -1 . e -x (1+e -x ) 2 dx -1 exp(-1 2 t) (1+exp(-1 2 t)) 2 Kotz, q, r, s > 0 sΓ( d 2 ) π d 2 Γ( 2q+d-2 2s ) r 2q+d-2 2s t q-1 exp(-rt s ) Unimodal GM 1 (2π) d 2 n k=1 π k θ d k exp - θ 2 k 2 t Laplace, λ > 0 2 λ(2π) d 2 K d 2 -1 2 λ t λ 2 t d 2 -1 Uniform GM Γ( d+1 2 ) √ 2π d 2 χ 2 d+1 (t) t d+1 2
Table 1. Some classical d-dimensional elliptical distributions with corresponding normalisation constants and generators

Cauchy distribution corresponds to a Student distribution with ν = 1 degree of freedom.

In this paper, we focus on consistent elliptical distributions. Unfortunately, Logistic and Kotz distributions do not have this property (except the Kotz distribution with s = q = 1, and r = 1 2 , i.e the Gaussian distribution). Therefore, our examples will only refer to the distributions mentioned in Table 2. For these models, the non negative random variable ε is given in Table 2.

Distribution

Gaussian We have seen that an elliptical distribution is characterized by parameters µ, Σ, and by either the characteristic function Φ, the radius R or the generator g d . For this reason, we define the distribution of an elliptical random vector by any of these three possible characterization, using indifferently the notation

1 Student, ν > 0 χν √ ν Unimodal Gaussian Mixture n k=1 π k δ θ k Laplace, λ > 0 1 E( 1 λ ) Uniform Gaussian Mixture U (]0, 1[)
X ∼ E d (µ, Σ, Φ), X ∼ E d (µ, Σ, R) or X ∼ E d (µ, Σ, g d ).
At last, in order to emphasis the role played by the radius and the dimension, we also use the denomination (R, d)-elliptical, as defined hereafter.

Definition 2.3. An elliptical random vector of R d with radius R is called (R, d)-elliptical.
The following proposition, from [START_REF] Hult | Multivariate extremes, aggregation and dependence in elliptical distributions[END_REF], is the basis of our study.

Proposition 2.4 (Affine transformation). Let X a consistent (R, d)-elliptical random vector with parameters µ and Σ. Then for any c ∈ R d , c T X is (R, 1)-elliptical with parameters c T µ and c T Σc. Proposition 2.4 implies that an affine transformation of a (R, d)-elliptical random vector is a (R, 1)elliptical random variable. The proposition below is a direct consequence of this result (see [START_REF] Hult | Multivariate extremes, aggregation and dependence in elliptical distributions[END_REF], for a proof).

Proposition 2.5 (Subvectors distributions). Let X = (X 1 , X 2 ) T be a consistent (R, d)-elliptical random vector with X 1 ∈ R d1 , X 2 ∈ R d2 , d 1 + d 2 = d and parameters µ and Σ. Let us write:

Σ = Σ 11 Σ 12 Σ 21 Σ 22 , µ = µ 1 µ 2 .
Then X 1 and X 2 are respectively (R, d 1 )-and (R, d 2 )-elliptical with parameters µ 1 , Σ 11 and µ 2 , Σ 22 , respectively.

Remark that a p-dimensional subvector of a (R, d)-elliptical random vector with the consistency property is (R, p)-elliptical. As a consequence when p = 1, all marginals are (R, 1)-elliptical. The following proposition gives some indications concerning the conditional distributions of elliptical vectors.

Proposition 2.6 (Conditional distribution). Let X = (X 1 , X 2 ) T be a consistent (R, d)-elliptical random vector with with X 1 ∈ R d1 , X 2 ∈ R d2 , d 1 + d 2 = d and parameters µ and Σ. Let us write:

(2.5) Σ = Σ 11 Σ 12 Σ 21 Σ 22 , µ = µ 1 µ 2 .
The conditional distribution X 2 |(X 1 = x 1 ) has parameters:

(2.6)

µ 2|1 = µ 2 + Σ 21 Σ -1 11 (x 1 -µ 1 ) Σ 2|1 = Σ 22 -Σ 21 Σ -1 11 Σ 12 Furthermore, X 2 |(X 1 = x 1
) is elliptical, with radius R * given by :

(2.7)

R * d = R 1 -β R βU (d) = C -1 11 (x 1 -µ 1 )
where C 11 is the Cholesky root of Σ 11 , and β ∼ Beta( d1 2 , d2 2 ).

At last, the conditional density of X 2 |(X 1 = x 1 ) is given by : q1) , and

(2.8) f X2|X1 (x 2 |x 1 ) = c 2|1 |Σ 2|1 | 1 2 g d q 1 + (x 2 -µ 2|1 ) T Σ -1 2|1 (x 2 -µ 2|1 ) with c 2|1 = c d c d 1 g d 1 (
q 1 = (x 1 -µ 1 ) T Σ -1 11 (x 1 -µ 1 ).
Proof. For Equations (2.6) and (2.7), a proof may be found in [START_REF] Frahm | Generalized Elliptical Distributions: Theory and Applications[END_REF]. Concerning Equation (2.8), the proof is obvious, since

f X2|X1 (x 2 |x 1 ) = f X 1 ,X 2 (x1,x2) f X 1 (x1)
. Then :

(2.9)

f X2|X1 (x 2 |x 1 ) = c d |Σ| 1 2 g d (x -µ) T Σ -1 (x -µ) c d 1 |Σ11| 1 2 g d1 (q 1 ) Since (x -µ) T Σ -1 (x -µ) = q 1 + (x 2 -µ 2|1 ) T Σ -1 2|1 (x 2 -µ 2|1 ) and |Σ| = |Σ 11 ||Σ 2|1 |,
we get the expected result.

We have introduced the main definitions and properties of elliptical distributions. With these tools, we can define elliptical random fields. Indeed, a random field (X(t)) t∈T is R-elliptical if ∀n ∈ N, ∀t 1 , ..., t n ∈ T , the vector (X(t 1 ), ..., X(t n )) is (R, n)-elliptical. Obviously, it implies that all the k-dimensional subvectors of (X(t 1 ), ..., X(t n )) are (R, k)-elliptical. This assumption corresponds exactly to consistent elliptical distributions properties given in Propositions 2.4, 2.5 and in Kano (1994). We thus focus our study on elliptical distributions with the consistency property. In the following section, we focus on conditional quantiles of elliptical distributions, applicated to our problem of spatial prediction.

Conditional quantiles

In this section, we give exact expressions for conditional quantiles of an elliptical random field, and discuss problems that occurs for computing these values. From now on, we consider the following context: (X(t)) t∈T is an R-elliptical random field. We consider N observations at locations t 1 , ..., t N ∈ T , called (X(t 1 ), ..., X(t N )). Our aim is to predict, at a site t ∈ T , the quantile of X(t) given X(t 1 ), ..., X(t N ). Remark that the vector (X(t), X(t 1 ), ..., X(t N )) is (R, N + 1)-elliptical. Thus, we can denote X 2 = X(t) ∈ R and X 1 = (X(t 1 ), ..., X(t N )) ∈ R N and restrict ourselves to the study of the q α (X 2 |X 1 ).

3.1. General case. Let us firstly give an expression of the theoretical conditional quantile. Consider the respective cumulative distribution functions Φ R and Φ * R ,

(3.1) Φ R (x) = P RU (1) ≤ x , Φ R * (x) = P R * U (1) ≤ x ,
where U (1) is 1 or -1 with probability 1 2 . Φ R and Φ R * are respective cumulative distribution functions of the reduced centered (R, 1)-elliptical random variable and (R * , 1)-elliptical random variable. With this notation and the conditional generator given in Equation (2.7), we get the following result.

Proposition 3.1 (Conditional elliptical quantile). Let X = (X 1 , X 2 ) T be a consistent (R, N +1)-elliptical random vector with X 1 ∈ R N , X 2 ∈ R and parameters µ and Σ. Let us write Σ:

Σ = Σ 11 Σ 12 Σ 21 Σ 22 Then the α-quantile of X 2 |(X 1 = x 1
) is given by :

(3.2) q α (X 2 |X 1 = x 1 ) = µ 2|1 + Σ 2|1 Φ -1 R * (α)
where µ 2|1 and Σ 2|1 are given in Equation (2.6), and R * in Equation (2.7).

Proof. We know that 1) . Then, our aim is to get q α such as :

X 2 |(X 1 = x 1 ) ∼ µ 2|1 + R * Σ 2|1 U (
P µ 2|1 + R * Σ 2|1 U (1) ≤ q α = α Since P µ 2|1 + R * Σ 2|1 U (1) ≤ q α = P R * U (1) ≤ qα-µ 2|1 √ Σ 2|1 = Φ R * qα-µ 2|1 √ Σ 2|1
, we easily get the expected result.

A general expression for conditional quantiles of consistent elliptical distributions is thus available. However, Equation (2.7) shows that the conditional radius R * is difficult to exploit and thus, the calculation of the conditional quantile is not tractable in general. Indeed, the term Φ -1 R * (α), is difficult to calculate, since the conditional density in Equation (2.8) leads to the inverse problem in x:

x -∞ c 2|1 g N +1 (q 1 + t 2 )dt = α
where c 2|1 is given in Proposition 2.6. Computing and inversing the function Φ R * from the distribution of R * is not easier, since the latter is hard to obtain. We thus have a general expression for the conditional quantile which is not easy to compute in the general case.

Fortunately, in several cases, we get an explicit formula for conditional quantiles. We propose two classical examples: Gaussian and Student distributions. We thus get explicit formulas for some simple cases, but it would not be possible in some other cases, such as Gaussian Mixture distribution.

3.2. Gaussian case. The Gaussian case is the most famous one and the simplest case of elliptical distribution. Since a conditional Gaussian distribution is still Gaussian, we have:

(3.3) X 2 |(X 1 = x 1 ) ∼ N (µ 2|1 , Σ 2|1 )
with µ 2|1 and Σ 2|1 given in Equation (2.6). Then, the calculation of the conditional α-quantile of X 2 |(X 1 = x 1 ) is immediate, and gives:

(3.4) q α (X 2 |X 1 = x 1 ) = µ 2|1 + Σ 2|1 Φ -1 (α)
where Φ is the standard Gaussian distribution function.

Student case.

In the Student case also, we get explicit formulas for conditional quantiles. For that purpose, we need to introduce some properties of Student distributions. The lemma below and associated proof may be found in [START_REF] Nadarajah | Multivariate T-Distributions and Their Applications[END_REF].

Lemma 3.2. Let X be a d-dimensional Student distribution with ν degree of freedom and parameters

µ ∈ R d and Σ ∈ R d×d . Consider X = (X 1 , X 2 ) T with X 1 ∈ R d1 , X 2 ∈ R d2 and d 1 + d 2 = d.
Then the density function of the conditional random variable X 2 |(X 1 = x 1 ) is given by the following equation.

(3.5) f X 2 |X 1 (x2|x1) = Γ ν+d 2 ((ν + d1)π) d 2 2 Γ ν+d 1 2 |Σ 2|1 | 1 2 1 + 1 ν q 2|1 (x2) 1 + 1 ν q1 -ν+d 2 × ν+d 1 ν 1 + 1 ν q1 d 2 2
, where q 2|1 (x 2 ) and q 1 are the Mahalanobis distances :

q 2|1 (x 2 ) = (x 2 -µ 2|1 ) T Σ -1 2|1 (x 2 -µ 2|1 ) , q 1 = (x 1 -µ 1 ) T Σ -1 11 (x 1 -µ 1 )
. With the conditional density, we can deduce the conditional α-quantile, in the following proposition.

Proposition 3.3 (Conditional Student quantile). Let X a N + 1-dimensional Student distribution with ν degree of freedom and parameters µ ∈ R N +1 and Σ ∈ R (N +1)×(N +1) . Consider X = (X 1 , X 2 ) T with X 1 ∈ R N and X 2 ∈ R 1 . The conditional α-quantile of X 2 |(X 1 = x 1 ) has the following expression (3.6) q α (X 2 |X 1 = x 1 ) = µ 2|1 + Σ 2|1 ν ν + N 1 + 1 ν q 1 Φ -1 ν+N (α) .
Proof. We consider the density (3.5), with

d 1 = N , d 2 = 1. We have q 2|1 (x 2 ) = (x2-µ 2|1 ) 2 Σ 2|1 ∈ R. Let q α be the conditional quantile of X 2 |(X 1 = x 1 ). It satisfies: qα -∞ f X2|X1 (x 2 |x 1 )dx 2 = α .
Hence :

Γ ν+N +1 2 ((ν + N )π) 1 2 Γ ν+N 2 |Σ 2|1 | 1 2 ν+N ν 1 + 1 ν q1 1 2 × qα -∞ 1 + 1 ν + N ν+N ν q 2|1 (x2) 1 + 1 ν q1 -ν+N +1 2 dx2 = α Considering z = √ ν+N ν (x2-µ 2|1 ) √ Σ 2|1 √ 1+ 1 ν q1
we obtain:

Γ ν+N +1 2 ((ν + N )π) 1 2 Γ ν+N 2 × zα -∞ 1 + z 2 ν + N -ν+N +1 2 dz = α, with z α = ν+N ν (q α -µ 2|1 ) Σ 2|1 1 + 1 ν q 1 .
Let Φ ν+N be the one dimensional Student distribution with N + ν degrees of freedom. It is obvious that

Φ ν+N   ν+N ν (q α -µ 2|1 ) Σ 2|1 1 + 1 ν q 1   = α.
As a consequence,

q α = µ 2|1 + ν ν + N Σ 2|1 1 + 1 ν q 1 Φ -1 ν+N (α).
We did not obtain such simple results for other elliptical distributions. In order to predict conditional quantiles for other elliptical distributions, we propose, in what follows, two approaches. In the following section, we apply quantile regression in the case of consistent elliptical distributions. In Section 5, we focus on extreme quantile levels, i.e when the quantile level α is close to 0 or 1.

Quantile regression

In this section, we propose quantile regression predictors for elliptical distributions. We give direct explicit expressions for these predictors as well as their distributions. Quantile regression, introduced by [START_REF] Koenker | Regression quantiles[END_REF], is a classical way to estimate conditional quantiles of a distribution. If

X 1 ∈ R N and X 2 ∈ R, the α-quantile of X 2 |(X 1 = x 1 ) is approximated by (4.1) qα (X 2 |X 1 = x 1 ) = β * T x 1 + β * 0 , where β * and β * 0 are solutions of the following minimization problem (4.2) (β * , β * 0 ) = arg min β∈R N ,β0∈R E φ α (X 2 -β T X 1 -β 0 ) .
and where the scoring function φ α (see [START_REF] Grant | Consistent scoring functions for quantiles[END_REF] is

(4.3) φ α (x) = (α -1)x1 {x<0} + αx1 {x>0} = αx -x1 {x<0} .
Obviously, given Equation (4.2), we need to assume the existence of moments E [X 1 ] and E [X 2 ]. It excludes some elliptical distributions with the consistency property (such as the Cauchy distribution). Most of the time, we need simulations to find the solution, and an appropriate algorithm, for example MM algorithm (see [START_REF] Hunter | Quantile regression via an mm algorithm[END_REF] or stochastic gradient (see [START_REF] Zheng | Gradient descent algorithms for quantile regression with smooth approximation[END_REF]. As an example, the package quantreg of R performs simulations of (X 1 , X 2 ) and solves Equation (4.2) by simplex or interior point algorithm. These simulations are computationally expensive. Furthermore, it is difficult to get the distribution of the estimator. It is a reason why we propose another approach for consistent elliptical distributions.

In the following, we denote by

E 1 c (N + 1) the set of consistent elliptical random vectors X in R N +1 , with X = (X 1 , X 2 ) T , X 1 ∈ R N and X 2 ∈ R
, having an order one moment. We shall get an explicit formula for β * and β * 0 . Let us recall the distribution function introduced in Equation (3.1), Φ R (x) = P RU (1) ≤ x . From Equation (2.8), we get

(4.4) Φ R (x) = x -∞ c 1 g 1 (x 2 )dx.
We will use this notation later in our formulas.

The next lemma is an explicit formula of a truncated moment for a bivariate elliptical distribution with the consistency property.

Lemma 4.1. Let (X, Y ) be a consistent (R, 2)-elliptical vector, which admits an order one moment and

with parameters µ = (µ X , µ Y ) T and Σ = σ 2 X ρσ X σ Y ρσ X σ Y σ 2 Y . Then (4.5) E X1 {Y >0} = µ X Φ R µ Y σ Y + ρσ X +∞ µ Y σ Y c 1 zg 1 (z 2 )dz.
Proof. Let f (X,Y ) be the density of the random vector (X, Y ). Then (4.6)

E X1 {Y >0} = +∞ -∞ +∞ 0 xf X,Y (x, y)dxdy.
We get:

+∞ -∞ xf X,Y (x, y)dx = f Y (y) +∞ -∞ x f X,Y (x, y) f Y (y) dx = f Y (y) +∞ -∞ xf X|Y (x|y)dx .
We recognize the integral expression of the conditional mean of X|Y , which is given in Equation (2.6) :

+∞ -∞ xf X,Y (x, y)dx = f Y (y)E [X|Y = y] = f Y (y) µ X + ρ σ X σ Y (y -µ Y ) .
Integrating with respect to the second variable y,

(4.7) E X1 {Y >0} = +∞ 0 f Y (y) µ X + ρ σ X σ Y (y -µ Y ) dy. Recall that f Y (y) = c1 σ Y g 1 (y-µ Y ) 2 σ 2 Y , this leads to E X1 {Y >0} = µ X +∞ 0 f Y (y)dy + ρ σ X σ Y c 1 +∞ 0 y -µ Y σ Y g 1 (y -µ Y ) 2 σ 2 Y dy.
Using the change of variable z = y-µ Y σ Y , we get

E X1 {Y >0} = µ X Φ R µ Y σ Y + ρσ X +∞ - µ Y σ Y c 1 zg 1 z 2 dz.
Hence the result.

Previous lemma is a keystone to find β * and β * 0 . Let ρ j be the correlation coefficient between X 1j and the random variable

X 2 -β * T X 1 -β * 0 . Indeed, X 1j is (R, 1)-elliptical with parameters µ 1j and σ 2 1j . Furthermore, X 2 -β * T X 1 -β * 0 is (R, 1)-elliptical too, with parameters µ 2 -β * T µ 1 -β * 0 and (-β * , 1) T Σ(-β * T , 1). Hence the (R, 2)-elliptical vector X 1j , X 2 -β * T X 1 -β *
0 admits as second parameter the matrix :

σ 2 1j ρ j σ 1j (-β * , 1) T Σ(-β * T , 1) ρ j σ 1j (-β * , 1) T Σ(-β * T , 1) (-β * , 1) T Σ(-β * T , 1) Proposition 4.2 (Explicit form of β * ). Let X = (X 1 , X 2 ) T ∈ E 1 c (N + 1).
The optimal β * is given by : (4.8)

β * = Σ -1 11 Σ 12 .
Proof.

(4.9) (β * , β * 0 ) = arg min β∈R N ,β0∈R E φ α X 2 -β T X 1 -β 0 Equation (4.9) is equivalent to : (4.10) E -X 1 φ α (X 2 -β * T X 1 -β * 0 ) = 0 E -φ α (X 2 -β * T X 1 -β * 0 ) = 0
with the scoring function derivative :

(4.11) φ α (x) = (α -1)1 {x<0} + α1 {x>0} = (α -1) + 1 {x>0}
We obtain

(4.12) (1 -α)E [X 1 ] -E X 1 1 {X2-β * T X1-β * 0 >0} = 0 (1 -α) -E 1 {X2-β * T X1-β * 0 >0} = 0 Since X 2 -β * T X 1 -β * 0 is (R, 1)-elliptical, with parameters µ = µ 2 -β * T µ 1 -β * 0 and Σ = (-β * , 1) T Σ(-β * , 1). Then E 1 {X2-β * T X1-β * 0 >0} = P X 2 -β * T X 1 -β * 0 > 0 , i.e, if we denote σ = Σ 1 2 : (4.13) E 1 {X2-β * T X1-β * 0 >0} = 1 -Φ R -µ σ = Φ R µ σ
Next, thanks to Equation (4.5), we have, for j ∈ {1, ..., N } :

(4.14) E X 1j 1 {X2-β * T X1-β * 0 >0} = µ 1j Φ R µ σ + ρ j σ 1j +∞ µ σ c 1 zg 1 (z 2 )dz
Finally, Equation (4.12) may be written as :

(4.15)            µ 1j (1 -α) -Φ R µ σ -ρ j σ 1j +∞ µ Σ c 1 zg 1 (z 2 )dz = 0 (1 -α) -Φ R µ σ = 0
, ∀j ∈ {1, ..., N } Therefore :

(4.16)

ρ j σ 1j +∞ µ Σ c 1 zg 1 (z 2 )dz = 0, ∀j ∈ {1, ..., N }
Then ρ j = 0, ∀j ∈ {1, ..., N }. Now, it is not difficult to see that β * = Σ -1 11 Σ 12 . Indeed, the second parameter of the (R, N + 1)-elliptical vector

X 1 , X 2 -β * T X 1 -β * 0 is : Σ 11 0 R N 0 T R N (-β * , 1) T Σ(-β * T , 1)
Hence its determinant is (-β * , 1) T Σ(-β * T , 1)|Σ 11 |. Using the conditional moments of X 2 -β * T X 1 -β * 0 |X 1 given in Proposition 2.6, we get the following equation :

Σ 2|1 = Σ 22 -Σ 21 Σ -1 11 Σ 12 = (-β * , 1) T Σ(-β * T , 1) Since (-β * , 1) T Σ(-β * T , 1) = Σ 22 + β * T Σ 11 β * -2β * T Σ 12
, the previous equation is equivalent to :

2β * T Σ 12 -β * T Σ 11 β * -Σ 21 Σ -1 11 Σ 12 = 0 The convex function 2β * T Σ 12 -β * T Σ 11 β * -Σ 21 Σ -1 11 Σ 12 reaches its minimum 0 at β * = Σ -1 11 Σ 12 .
This result shows that a quantile regression through a linear model is meaningless. Indeed, the vector Σ -1 11 Σ 12 does not depend on α and corresponds to the vector of weights in the Kriging mean (cf. [START_REF] Cressie | Spatial prediction and ordinary kriging[END_REF]. Using a linear predictor, the result of the quantile regression, whatever the quantile level α ∈ [0, 1], would be the conditional mean µ 2|1 given in Equation (2.6). This is why the addition of the affine term β * 0 is required. Let us now calculate this term β * 0 . As mentioned in the last proof that, obviously, X 2 -β * T X 1 -β * 0 is elliptical with the same radius as (X 1 , X 2 ). With Equation (4.8), we are now able to prove that its second parameter is Σ 2|1 , given in Equation (2.6).

Lemma 4.3. Let X = (X 1 , X 2 ) T ∈ E 1 c (N + 1). X 2 -β * T X 1 is (R, 1)-elliptical with parameters µ 2 - Σ 21 Σ -1 11 µ 1 and Σ 2|1 = Σ 22 -Σ 21 Σ -1 11 Σ 12 .
Proof. The random variable X 2 -β * T X 1 , as an affine transformation of a consistent (R, d)-elliptical random vector, is obviously (R, 1)-elliptical. Furthermore, with the expression of β * given in Proposition 4.2,

E X 2 -β * T X 1 = µ 2 -Σ 21 Σ -1 11 µ 1 . Concerning the second parameter, it is equal to (-β * , 1) T Σ(-β * , 1). Thus, (-β * , 1) T Σ(-β * , 1) = σ 2 2 -2 n i=1 β * i ρ ix σ 1i σ 2 + n j=1 β * j n i=1 β * i ρ ij σ 1i σ 1j .
It may be rewritten in matrix form

(-β * , 1) T Σ(-β * , 1) = σ 2 2 -2β * T Σ 12 + β * T Σ 11 β * . Using that β * = Σ -1
11 Σ 12 , we prove the following equality.

(-β * , 1) T Σ(-β * , 1) = σ 2 2 -2Σ 21 Σ -1 11 Σ 12 + Σ 21 Σ -1 11 Σ 11 Σ -1 11 Σ 12 = σ 2 2 -Σ 21 Σ -1
11 Σ 12 , which concludes the proof.

We have seen that the conditional second parameter of our (R, 1)-elliptical random variable is exactly the conditional second parameter Σ 2|1 of Equation (2.6). We are led to our main result of this section on quantile regression for elliptical distributions.

Theorem 4.4 (Quantile Regression Predictor). Let X = (X 1 , X 2 ) T ∈ E 1 c (N +1). The quantile regression vector (β * , β * 0 ) of X 2 |(X 1 = x 1 ), satisfying Equation (4.
2), is given by (4.17)

β * = Σ -1 11 Σ 12 β * 0 = µ 2 -Σ 21 Σ -1 11 µ 1 + σ 2|1 Φ -1 R (α)
The Quantile Regression Prediction with level α ∈ [0, 1] is given by:

(4.18) qα (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 R (α) Furthermore, the distribution of the associated predictor qα (X 2 |X 1 ) is (4.19) qα (X 2 |X 1 ) ∼ E 1 µ 2 + σ 2|1 Φ -1 R (α), Σ 21 Σ -1 11 Σ 12 , g 1
Proof. β * is given in Equation (4.8). Concerning β * 0 , it satisfies:

β * 0 = arg min β0∈R E φ α X 2 -β * T X 1 -β 0 .
Finally, β * 0 is the α-quantile of the random variable X 2 -β * T X 1 . We have seen in Lemma 4.3, that X 2 -β * T X 1 is (R, 1)-elliptical with parameters µ 2 -Σ 21 Σ -1 11 µ 1 and Σ 2|1 . Then, using the quantile formula of an elliptical distribution, and denoting σ 2|1 = Σ 2|1 , we get:

β * 0 = µ 2 -Σ 21 Σ -1 11 µ 1 + σ 2|1 Φ -1 R (α).
With our optimal parameters β * and β * 0 , we can now express our Quantile Regression Prediction of X 2 given X 1 = x 1 , qα (X 2 |X 1 = x 1 ), using Equation (4.1):

qα (X 2 |X 1 = x 1 ) = β * x 1 + β * 0 = µ 2 + Σ 21 Σ -1 11 (x 1 -µ 1 ) + σ 2|1 Φ -1 R (α).
We recognize, on the left, the expression of µ 2|1 given in Equation (2.6).

Concerning its distribution, we know that

qα (X 2 |X 1 ) = µ 2 + Σ 21 Σ -1 11 (X 1 -µ 1 ) + σ 2|1 Φ -1 R (α), with X 1 ∼ E N (µ 1 , Σ 11 , g N ).
As an affine combination of X 1 , using Proposition 2.4, we have the relation-

ship qα (X 2 |X 1 ) ∼ E 1 µ 2 + σ 2|1 Φ -1 R (α), (Σ 21 Σ -1 11 )Σ 11 (Σ -1 11 Σ 12 ), g 1 . Hence the result.
As an illustration of Theorem 4.4, we now propose several examples, such as Gaussian and Student cases, where the theoretical conditional quantiles are known (see Section 3). We also consider the Unimodal Gaussian Mixture and Laplace cases, where we do not have any explicit formula for theoretical quantiles. Numerical applications are proposed in the last section.

4.1. Gaussian case. We recall the theoretical formula of the conditional quantile of a Gaussian distribution :

(4.20) q α (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 (α).
In that case, the Quantile Regression Prediction qα (X 2 |X 1 = x 1 ) is exactly the same :

(4.21) qα (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 (α).
For Gaussian processes, the Quantile Regression Predictor and the theoretical conditional quantile coincide.

4.2. Student case. The Gaussian case is the most simple one. Unfortunately, for Student distributions, Quantile Regression Predictor and conditional quantile do not coincide. Indeed, the conditional quantile of a multivariate Student distribution is given by Equation (3.6):

q α (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 ν ν + N 1 + 1 ν (x 1 -µ 1 ) T Σ -1 11 (x 1 -µ 1 )Φ -1 ν+N (α). The Quantile Regression Prediction is qα (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 ν (α).
The error q α (X 2 |X 1 = x 1 ) -qα (X 2 |X 1 = x 1 ) that is done when using the Quantile Regression Prediction depends on the Mahalanobis distance (x 1 -µ 1 ) T Σ -1 11 (x 1 -µ 1 ) which may be large (see numerical study in Section 5). 4.3. Unimodal Gaussian Mixture case. We consider here a mixture of 2 centered Gaussian distributions, i.e a radius R such that

R = χ d pθ 1 + (1 -p)θ 2 .
Immediately, we get the following multivariate p.d.f :

(4.22) f X (x) = 1 (2π) d 2 pθ d 1 e -θ 2 1 2 (x-µ) T Σ -1 (x-µ) + (1 -p)θ d 2 e -θ 2 2 2 (x-µ) T Σ -1 (x-µ) .
As an illustration, we propose some examples of univariate densities (4.22), with different values of p, θ 1 and θ 2 : For these models, we do not have any explicit formula for q α (X 2 |X 1 = x 1 ). On the other hand, using the conditional density expressed in Equation (2.8), we may get a numerical approximation of this value, and compare it with the Quantile Regression Prediction. (4.23)

f (x) = 1 |Σ| 1 2 (2π) d 2 2 λ K d 2 -1 2 λ (x -µ) T Σ -1 (x -µ) λ 2 (x -µ) T Σ -1 (x -µ) d 2 -1 , ∀x ∈ R d \{µ}.
where K m (x) denotes the modified Bessel function of the second kind and order m, evaluated at x. In this case, we can easily give an expression for the Quantile Regression Prediction qα (X 2 |X 1 = x 1 ). This is the aim of the following proposition.

Proposition 4.5. In the Laplace case, the Quantile Regression Prediction is given by:

(4.24)    qα (X 2 |X 1 = x 1 ) = µ 2|1 + λ 2 σ 2|1 ln(2α), for α ≤ 1 2 qα (X 2 |X 1 = x 1 ) = µ 2|1 - λ 2 σ 2|1 ln (2(1 -α)) , for α > 1 2 .
Proof. From (4.18), we just have to calculate Φ -1 R (α), where Φ R is the distribution function of the univariate reduced and centered Laplace distribution, with density function (4.23), and with

d = 1. Then, Φ -1 R (α) equals x ∈ R so that x -∞ 1 √ 2λ e - √ 2 λ |t| dt = α.
Since the univariate reduced centered Laplace distribution is symmetric, we have the equivalence α ≤ 1 2 ⇔ x ≤ 0. Then, we separate the cases α ≤ 1 2 and α > 1 2 . If α ≤ 1 2 , we have to solve

x -∞ 1 √ 2λ e √ 2 λ t dt = α ⇒ x = λ 2 ln(2α).
If α > 1 2 , we write :

0 -∞ 1 √ 2λ e √ 2 λ t dt + x 0 1 √ 2λ e - √ 2 λ t dt = α. Since 0 -∞ 1 √ 2λ e √ 2 λ t dt = 1 2 and x 0 1 √ 2λ e - √ 2 λ t dt = 1 2 - 1 2 e - √ 2 λ x , we quickly get x = -ln (2(1 -α)). Finally,    qα (X 2 |X 1 = x 1 ) = µ 2|1 + λ 2 σ 2|1 ln(2α) , α ≤ 1 2 qα (X 2 |X 1 = x 1 ) = µ 2|1 - λ 2 σ 2|1 ln (2(1 -α)) , α > 1 2
which is the announced result.

Extremal quantiles

In this section, we show that, for an elliptical distribution, best affine predictors of conditional quantiles are inadequate for extreme quantile levels. We thus propose some new extremal quantiles, which are shown to be equivalent to the true conditional quantiles. In order to simplify the notations, we denote by q 1 the Mahalanobis distance (x 1 -µ 1 ) T Σ -1 11 (x 1 -µ 1 ). We have previously defined a predictor qα based on quantile regression. We have seen that the best affine prediction of conditional quantiles of an R-elliptical distribution with the consistency property is given by:

(5.1) qα (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 R (α) where Φ R (x) = x -∞ c 1 g 1 (t 2 )dt.
The conditional quantile is:

(5.2)

q α (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 R * (α)
with the conditional radius R * defined in Equation (2.7). In this section, the aim is to establish a relation between Φ -1 R and Φ -1 R * for extremal values of α, i.e for α → 0 or α → 1. Extreme quantiles estimation is an active research topic, and we can find several papers in the literature. We think particularly to El Methni et al. ( 2012), but the random variables X 1 , ..., X n are supposed i.i.d. In our case, the covariates are dependent. More recently, [START_REF] Gong | Estimation of extreme quantiles for functions of dependent random variables[END_REF] proposed an approach for dependant random variables, but it requires simulations. It is why, in this paper, we propose another approach. To this aim, we need a kind of regular variations hypothesis. It seems important to note that all elliptic c.d.f Φ R are not regularly varying.

Assumption 1. Their exist 0 < < +∞ and γ ∈ R such as :

(5.3) lim x→+∞ Φ R * (x) Φ R (x γ ) =
where Φ = 1 -Φ is the survival function associated to Φ.

We also recall that the cumulative distribution functions are Φ R (

x) = x -∞ c 1 g 1 (t 2 )dt and Φ R * (x) = x -∞
c * 1 g N +1 (q 1 + t 2 )dt, with c * 1 given by Equation (2.8).

(5.4)

c * 1 = c N +1 c N g N (q 1 )
,

where c N +1 and c N are the normalization coefficients of the elliptical distributions with radius R in dimensions N + 1 and N . Then, coefficients γ and may also be obtained as

(5.5) lim x→+∞ c * 1 g N +1 (q 1 + x 2 ) c 1 γx γ-1 g 1 (x 2γ ) =
Table 3 gives values of γ and for the examples introduced in Section 2. The calculations will be detailed later on. Unfortunately, the Laplace distribution does not satisfy Assumption 1, this is why we will not consider this example in what follows.

Thanks to Equation (5.3), we have, under Assumption 1, the following equivalence when x goes to ∞ :

(

5.6) Φ R * (x) ∼ x→∞ Φ R (x γ )
Our aim is now to get an equivalence relationship between the quantile function Φ -1 R * and Φ -1 R * . For that purpose, we refer to the paper of Djurcić and Torgasev, which gives some conditions to deduce the equivalence of inverse functions if these functions are equivalent (cf. [START_REF] Djurcić | Strong asymptotic equivalence and inversion of functions in the class kc[END_REF].

Definition 5.1. A function f is a ϕ-function if f : [0, +∞[→ [0, +∞[, f (0) = 0, f is continuous, non decreasing on [0, +∞[, and f → +∞ when x → +∞. Distribution γ Gaussian 1 1 Student, ν > 0 N +ν ν Γ( ν+N +1 2 )Γ( ν 2 ) Γ( ν+N 2 )Γ( ν+1 2 ) 1 + q 1 ν N +ν 2 ν N 2 +1 ν+N Unimodal Gaussian Mixture 1 min(θ 1 ,...,θn) N exp - min(θ 1 ,...,θn) 2 2 q 1 n k=1 π k θ N k exp - θ 2 k 2 q 1 Uniform Gaussian Mixture N + 1 Γ( N +2 2 )q N +1 2 1 √ 2 Γ( N +1 2 )(N+1)χ 2 N +1 (q 1 )
Table 3. Coefficients γ and for classical distributions, where

q 1 = (x 1 - µ 1 ) T Σ -1 11 (x 1 -µ 1 ) is a function of x 1 . (see Proposition 2.6).
Clearly, our two equivalent functions Φ R * (x) and lΦ R (x γ ) are not ϕ-functions for several reasons :

Φ R * (0) = 1 2 , Φ R (0) = 1 2 , lim x→+∞ Φ R * (x) = 1 and lim x→+∞ lΦ R (x γ ) = < +∞.
Then, we have to transform these functions in order to find an equivalence. But before, we need more definitions and properties.

The following is the definition of a general class of functions K c , which contains in particular Regularly Varying functions.

Definition 5.2. K c is the set of all ϕ-functions f with the property :

(5.7) lim

x→+∞ λ→1 f (λx) f (x) = 1
In fact, the K c is the set of ϕ-functions for which the limit lim

x→+∞ λ→1 f (λx) f (x)
exists, because if it exists, then automatically we can swap the limits, and lim

x→+∞ λ→1 f (λx) f (x) = lim x→+∞ lim λ→1 f (λx) f (x)
= 1. We can also notice that the condition in Equation (5.7) is more general than the regular variation functions in the sense of Karamata.

These two last definitions are very important for the following lemma, which establishes the relation between the equivalence of two functions and equivalence of their inverses.

Lemma 5.1. [START_REF] Djurcić | Strong asymptotic equivalence and inversion of functions in the class kc[END_REF]] Suppose that f and g are two strictly increasing ϕ-functions, and that at least one of the functions f -1 , g -1 belongs to the class K c , and f (x) ∼ x→∞ g(x). Then

f -1 (x) ∼ x→∞ g -1 (x)
Now, as we said previously, the first step is to transform our fonctions Φ R * (x) and Φ R (x γ ) into ϕ-functions.

Definition 5.3. Let Ψ and Ψ * be

(5.8) Ψ(x) = 1 Φ R (x γ ) -2 Ψ * (x) = 1 Φ R * (x) -2 Lemma 5.2. If γ > 0, then Ψ and Ψ * are ϕ-functions.
Proof. Since Φ R and Φ R * are clearly strictly decreasing, Ψ and Ψ * are obviously strictly increasing. Furthermore, Φ R (0

) = Φ R * (0) = 1 2 , thus Ψ(0) = 0 Ψ * (0) = 0 . Finally, since lim x→+∞ Φ R (x γ ) = lim x→+∞ Φ R * (x) = 0, functions Ψ and Ψ * are such that lim x→+∞ Ψ(x) = lim x→+∞ Ψ * (x) = +∞. Now, given that Φ R * (x) ∼ x→∞ Φ R (x γ ), we have to check whether Ψ(x) ∼ x→∞ Ψ * (x).
Lemma 5.3. Under Asumption 1, we have

(5.9) Ψ * (x) ∼ x→∞ Ψ(x)
Proof. Let us focus on the limit

lim x→+∞ Ψ * (x) Ψ(x) = lim x→+∞ 1 Φ R * (x) -2 1 Φ R (x γ ) -2 . This limit is equal to lim x→+∞ Φ R (x γ ) Φ R * (x) • 1 -2Φ R * (x) 1 -2Φ R (x γ ) .
Thanks to (5.3), lim

x→+∞ Φ R (x γ ) Φ R * (x) = 1. Furthermore, lim x→+∞ Φ R * (x) = lim x→+∞ Φ R (x γ ) = 0. Then : lim x→+∞ Ψ * (x) Ψ(x) = lim x→+∞ Φ R (x γ ) Φ R * (x) • 1 -2Φ R * (x) 1 -2Φ R (x γ ) = 1.
In other words, Ψ * (x)

∼ x→∞ Ψ(x).
In order to use Lemma 5.1, we shall need an additional assumption on Ψ -1 * . Assumption 2. Ψ -1 belongs to the class K c .

The following proposition is the key step to predict conditional quantiles.

Proposition 5.4. Under Assumptions 1 and 2, we have

(5.10) Φ -1 R * (α) ∼ α→1 Φ -1 R 1 - 1 1-α + 2(1 -) 1 γ
Proof. Since Ψ -1 (x) belongs to the K c class, and Ψ * (x) and Ψ(x) are two strictly increasing ϕ-functions, and Ψ * (x) ∼ x→+∞ Ψ(x), then Lemma 5.1 gives

Ψ -1 * (x) ∼ x→+∞ Ψ -1 (x)
In other words

(5.11) Φ -1 R * 1 - 1 x + 2 ∼ x→+∞ Φ -1 R 1 - 1 x + 2 1 γ .
This may also be rewritten as

Φ -1 R * (α) ∼ α→1 Φ -1 R 1 - 1 1-α + 2(1 -) 1 γ with α = 1 -1 x+2
, hence the result. This result leads us to the following conditional extremal quantile prediction.

Definition 5.4 (Extreme Conditional Quantiles Predictors). Define (5.12)

       qα↑ (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 R 1 - 1 1-α +2(1-) 1 γ qα↓ (X 2 |X 1 = x 1 ) = µ 2|1 -σ 2|1 Φ -1 R 1 - 1 α +2(1-) 1 γ
The following proposition shows that qα↑ (X 2 |X 1 = x 1 ) and qα↓ (X 2 |X 1 = x 1 ) are asymptotically equivalent to the theoretical quantile q α (X 2 |X 1 = x 1 ), respectively for α → 1 and α → 0.

Theorem 5.5 (Equivalences of Extreme Conditional Quantiles Predictions). Let (X 1 , X 2 ) T ∈ E 1 c (N +1). Under Assumptions 1 and 2,

(5.13)    qα↑ (X 2 |X 1 = x 1 ) ∼ α→1 q α (X 2 |X 1 = x 1 ) qα↓ (X 2 |X 1 = x 1 ) ∼ α→0 q α (X 2 |X 1 = x 1 ) Proof. Recall Equation (3.2): q α (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1
R * (α) From (5.10) and (5.12), we immediatly prove the first half of the result:

q α (X 2 |X 1 = x 1 ) ∼ α→1 µ 2|1 + σ 2|1 Φ -1 R 1 - 1 1-α +2(1-) 1 γ = qα↑ (X 2 |X 1 = x 1 )
Now, let us prove the equivalence when α → 0. By the symmetry properties of elliptical distributions, we have Φ

-1 R * (α) = -Φ -1 R * (1 -α), ∀α ∈ [0, 1]. Then, q α (X 2 |X 1 = x 1 ) = µ 2|1 -σ 2|1 Φ -1 R * (1 -α).
Using the equivalence (5.10) and Equation (5.12), we finally get:

q α (X 2 |X 1 = x 1 ) ∼ α→0 µ 2|1 -σ 2|1 Φ -1 R 1 - 1 α +2(1-) 1 γ = qα↓ (X 2 |X 1 = x 1 )
We finish this section by checking Assumption 2 is satisfied for examples that are given in Table 3. 5.1. Gaussian case.

Lemma 5.6. Gaussian distribution satisfies Assumptions 1 and 2, with = 1 and γ = 1.

Proof. Firstly, let us calculate and γ. Take γ = 1, it remains the limit of Assumption 1 :

lim x→+∞ c * 1 g N +1 (q1+x 2 ) c1g1(x 2 )
The values c 1 , g N +1 (q 1 + x 2 ) and g 1 (x 2 ) may be easily deduced from Table 1

. Let c * 1 = c N +1 c N g N (q1) , = lim x→+∞ exp -q1+x 2 2 exp -q1 2 exp -x 2 2 = 1.
It remains to check whether the limit

l = lim x→+∞ λ→1 Ψ -1 (λx) Ψ -1 (x) = lim x→+∞ λ→1 Φ -1 1 -1 λx+2 Φ -1 1 -1 x+2
exists, or equivalently, whether the limit = lim (δ,y)→(0,0)

Φ -1 1 - y 2y+δ+1 Φ -1 1 -y 2y+1
exists. Now, we move to polar coordinates, i.e take δ = r cos(θ) and y = r sin(θ), and calculate the limit when r → 0. If this limit do not depend on θ, then the limit when (δ, y) → (0, 0) exists. Consider

lim r→0 Φ -1 1 - r sin(θ) 2r sin(θ)+r cos(θ)+1 Φ -1 1 -r sin(θ) 2r sin(θ)+1 = f (θ).
Clearly, the numerator and denominator both tend to +∞, then we use the L'Hôpital's rule to get

f (θ) = lim r→0 - sin(θ)(2r sin(θ)+r cos(θ)+1)-r sin(θ)(2 sin(θ)+cos(θ)) (2r sin(θ)+r cos(θ)+1) 2 1 φ ( Φ -1 ( 1- r sin(θ)
2r sin(θ)+r cos(θ)+1 ))

sin(θ)(2r sin(θ)+1)-r sin(θ)2 sin(θ) (2r sin(θ)+1) 2 1 φ ( Φ -1 ( 1-r sin(θ) 2r sin(θ)+1 )) On the left, the bulky ratio clearly tends to 1 when r → 0. On the right, we recall that φ is the density function of the standard normal distribution, i.e φ(x) = 1 √ 2π e -x 2 2 . Then :

= lim r→0 φ(Φ -1 (1-r sin(θ) 2r sin(θ)+1 )) φ(Φ -1 (1- r sin(θ)
2r sin(θ)+r cos(θ)+1 )) Now, the numerator and denominator both tend to 0. Once again, we use the L'Hôpital's rule and then,

f (θ) = lim r→0 f 1 (r, θ) f 2 (r, θ) φ Φ -1 1 - r sin(θ) 2r sin(θ) + r cos(θ) + 1 φ Φ -1 1 - r sin(θ) 2r sin(θ) + 1 φ Φ -1 1 -r sin(θ) 2r sin(θ)+1 φ Φ -1 1 - r sin(θ) 2r sin(θ)+r cos(θ)+1 . where       
f 1 (r, θ) = -sin(θ)(2r sin(θ) + 1) -r sin(θ)2 sin(θ) (2r sin(θ) + 1) 2 , f 2 (r, θ) = -sin(θ)(2r sin(θ) + r cos(θ) + 1) -r sin(θ)(2 sin(θ) + cos(θ)) (2r sin(θ) + r cos(θ) + 1) 2 ,

We have seen that the ratio f1(r,θ) f2(r,θ) tends to 1. Furthermore, we have the relationship φ (x) = -xφ(x). Finally,

f (θ) = lim r→0 Φ -1 1 -r sin(θ) 2r sin(θ)+1 Φ -1 1 - r sin(θ) 2r sin(θ)+r cos(θ)+1 = 1 f (θ) ,
thus, f (θ) = 1, hence the result.

Since Ψ -1 (x) belongs to the K c class and = γ = 1, we have the following formulas for the Extremal Predictions qα↑ and qα↓ :

qα↑

(X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1 R (α) qα↓ (X 2 |X 1 = x 1 ) = µ 2|1 + σ 2|1 Φ -1
R (α) Remark that qα↑ and qα↓ equal to the Quantile Regression Predictor qα , which equals the theoretical quantile q α . 5.2. Student case. We have seen in Section 4 that the Quantile Regression Predictor qα was not good in the Student case. Then, for extremal levels of quantile, we can guess that the predictor qα will be worse and worse if α is close to 0 or 1. It is a reason why, in this case, our extremal predictors qα↑ and qα↓ may potentially be especially interesting, and improve the prediction. Firstly, we prove that Assumptions 1 and 2 are filled.

Lemma 5.7. Student distributions satisfies Assumptions 1 and 2, with γ = N +ν ν and given in Table 3.

Proof. Concerning Assumption 1, we have to calculate the limit :

lim x→+∞ c * 1 g N +1 (q1+x 2 ) c1γx γ-1 g1(x 2γ )
From Table 1 and(5

.4), we deduce c 1 = Γ( 1+ν 2 ) √ νπΓ( ν 2 ) , c * 1 = Γ( N +1+ν 2 ) √ νπΓ( N +ν 2 ) 1 + q1 ν N +ν 2 , g N +1 (q 1 + x 2 ) = 1 + q1+x 2 ν -N +1+ν 2 and g 1 (x 2γ ) = 1 + x 2γ ν - 1+ν
2 . Then the previous limit is the following :

lim x→+∞ Γ( N +1+ν 2 )Γ( ν 2 ) Γ( 1+ν 2 )Γ( N +ν 2 ) 1 + q1 ν N +ν 2 1+ q 1 +x 2 ν -N +1+ν 2 γx γ-1 1+ x 2γ ν -1+ν 2
From now, we consider γ = N +ν ν and we focus on the ratio on the right. The leading monomials are ν N +1+ν 2

x -N -1-ν (numerator) and ν 1+ν 2 x -N -1-ν (denominator). Then the limit is simply ν

N 2 . Finally : = Γ( N +1+ν 2 )Γ( ν 2 ) Γ( 1+ν 2 )Γ( N +ν 2 ) 1 + q1 ν N +ν 2 ν N 2 +1
N +ν

Concerning Assumption 2, the proof is similar to the Gaussian case, but we have γ = 1 and = 1 here. We are dealing with the existence and value of the limit:

lim x→+∞ λ→1 Φ -1 R 1 -1 λx+2 1 γ Φ -1 R 1 -1 x+2 1 γ .
As in the Gaussian case, we use the polar coordinates and the L'Hôpital's rule to get the result.

In the case of Student distributions, we shall apply the predictions, with the coefficients and γ given in Table 3 :

   = Γ( ν+N +1 2 )Γ( ν 2 ) Γ( ν+N 2 )Γ( ν+1 2 ) 1 + q1 ν N +ν 2 ν N 2 +1 ν+N γ = N +ν ν 5.3. Unimodal Gaussian Mixture case.
Lemma 5.8. Unimodal Gaussian Mixture distribution satisfies Assumptions 1 and 2, with γ = 1 and given in Table 3.

Proof. Firstly, to prove that Assumption 1 is satisfied, we take γ = 1 and calculate the limit :

= lim x→+∞ c * 1 g N +1 (q1+x 2 ) c1g1(x 2 )
According to Table 1 andEquation (5.4

), c * 1 = 1 √ 2π n k=1 π k θ N k exp - θ 2 k 2 q1 , g d (t) = n k=1 π k θ d k e -θ 2 k 2 t and c 1 = 1 √ 2π . Hence the limit : = lim x→+∞ n k=1 π k θ N +1 k exp - θ 2 k 2 q 1 exp - θ 2 k 2 x 2 n k=1 π k θ N k exp - θ 2 k 2 q 1 n k=1 π k θ k exp - θ 2 k 2 x 2 .
Let k * such as θ k * = min(θ 1 , ..., θ n ). After simplifications, it remains :

= θ N k * exp - θ 2 k * 2 q1 n k=1 π k θ N k exp - θ 2 k 2 q1
Now we have to check Assumption 2, i.e prove that :

lim x→+∞ λ→1 Φ -1 R (1-1 λx+2 ) Φ -1 R (1-1 x+2 ) = 1
Like the Gaussian case, we do the same changes of variables, and the limit becomes tend to +∞ when r → 0, asymptotically, we only have to consider the terms of the sum with the largest coefficient -

= lim r→0 φ R (Φ -1 R (1-r sin(θ) 2r sin(θ)+1 )) φ R (Φ -1 R (1-
θ 2 k 2 , i.e the smallest θ k . Let k * such as θ k * = min(θ 1 , ..., θ n ). We get : = lim r→0 π * k θ * k φ(θ * k Φ -1 R (1-r sin(θ) 2r sin(θ)+1 )) π * k θ k * φ(θ k * Φ -1 R (1- r sin(θ) 2r sin(θ)+r cos(θ)+1 )) = lim r→0 φ(θ * k Φ -1 R (1-r sin(θ) 2r sin(θ)+1 )) φ(θ k * Φ -1 R (1- r sin(θ) 2r sin(θ)+r cos(θ)+1 ))
Finally, we have already calculated this kind of limit in the Gaussian case, and proved that it was equal to 1.

5.4. Uniform Mixture case. In this case, the radius

R d = χ d U (]0,1[) . The density obtained is then (5.14) f (x) = 1 |Σ| 1 2 Γ d+1 2 √ 2π d 2 χ 2 d+1 (x -µ) T Σ -1 (x -µ) [(x -µ) T Σ -1 (x -µ)] d+1 2 , ∀x ∈ R d \{µ}
where χ 2 m (x) denotes the c.d.f of the χ 2 distribution, with m degrees of freedom, evaluated at x.

Lemma 5.9. Uniform Mixture distribution satisfies Assumptions 1 and 2, with γ = N + 1 and given in Table 3.

Proof. Concerning Assumption 1, we consider γ = N + 1, and calculate the limit :

= lim x→+∞ c * 1 g N +1 (q1+x 2 ) c1γx γ-1 g1(x 2 )
According to Table 1 andEquation

(5.4), g d (t) = χ 2 d+1 (t) t d+1 2 , c 1 = 1 √ 2π and c * 1 = Γ( N +2 2 )q N +1 2 1 Γ( N +1 2 ) √ πχ 2 N +1 (q1) . Hence the limit = lim x→+∞ √ 2π Γ N +2 2 q N +1 2 1 Γ N +1 2 √ πχ 2 N +1 (q 1 ) χ 2 N +2 q 1 + x 2 χ 2 2 (x 2N +2 ) x 2N +2 (q 1 + x 2 ) N +2 2 1 (N + 1)x N , which leads to = Γ N +2 2 q N +1 2 1 √ 2 Γ N +1 2 (N + 1)χ 2 N +1 (q 1 )
.

To prove that Assumption 2 is filled, we have, as in the other cases, to calculate the limit :

lim x→+∞ λ→1 Φ -1 R 1 -1 λx+2 Φ -1 R 1 -1 x+2 = lim r→0 φ R Φ -1 R 1 -r sin(θ) 2r sin(θ)+1 φ R Φ -1 R 1 - r sin(θ) 2r sin(θ)+r cos(θ)+1 = f (θ).
We get once more that f (θ) = 1.

In this section, we introduced asymptotic predictions based solely on the function Φ R , that is assumed to be known. These predictions are asymptotically equivalent to the theoretical conditional quantiles, for levels of quantiles close to 0 or 1. Thus, they are expected to improve the quantile regression. In the following section, we perform some numerical applications to study the prediction quality.

Numerical study

In order to give a visual overview of the predictors we have defined, we have plotted in Figure 2, the conditional quantiles of an elliptical process observed at N = 5 points. We call X 1 ∈ R 5 the covariates vector. For x ∈ R, X 2 denotes the process at x and the aim is to predict the quantile of X 2 |X 1 = x 1 . For simplicity, we assume that the process is centered, and stationary (matrices Σ and Σ 11 are obtained through an exponential kernel). But our results would be applicable without these assumptions. Parameters for the Gaussian Mixture example are θ 1 = 1, θ 2 = 3 and p = 0.1.

Of course, for the Gaussian process, the curves coincide. For the other examples, Quantile Regression Predictors seem very far from the theoretical curves, especially in the Uniform Mixture case. On the other hand, Extremal predictors seem significantly better; they look closer to the target conditional quantiles here. We propose to use the following RM SE in order to quantify the error. (6.1) RM SE(q α ) = 1 n n i=1 q α (X

(i) 2 |X 1 = x 1 ) -qα (X (i) 2 |X 1 = x 1 ) 2 .
In the same way, we define the RM SE for Extremal Predictor :

(6.2) RM SE( qα ) = 1 n n i=1 q α (X (i) 2 |X 1 = x 1 ) -qα (X (i) 2 |X 1 = x 1 ) 2 .
The RMSE measures the average error in the prediction of the conditional quantiles. Table 4. r( qα ) = RM SE( qα ) and r(q α ) = RM SE(q α ), for different levels of α, and different consistent elliptical distributions. For non Gaussian distributions, the RM SE( qα ) decreases when α is close to 1, while the RM SE(q α ) increases, which is the expected behavior. We also remark that for the median (α = 0.5), the RMSE are always equal to 0. Indeed, if we replace α by 1 2 in Equation (5.12) and Equation (4.18), our two predictions are equal to the conditional mean µ 2|1 , defined in Proposition 2.6. Since median and mean are coincident for elliptical distributions, we exactly predict the conditional α-quantile for α = 1 2 . In order to have a visual overview of these phenomena, we propose, in Figure 3, some Q-Q plots, where we can observe the convergence of our extreme predictor for extreme values. On the other hand, the Quantile Regression Predictor seems better for reasonable values of α, whereas it is clearly biased. In this case, we can also use simulations, using the conditional density given in Equation (2.8). Of course, this approach is not possible for extreme values of α, and our extremal predictor takes its interest.

Conclusion

In this paper, we focused on conditional quantiles prediction, for elliptical random fields with the consistency property. We proposed two different methods. The first one is to use quantile regression, i.e to express the conditional quantile as an affine transformation of the observed values. This approach is widely used in the literature but it often requires a large number of simulations, especially for extreme levels of quantile (when α → 0 or α → 1). We have seen, in a first time, that we can obtain some explicit formulas in our case of consistent elliptical random fields (4.1). Furthermore, we have given the distribution of the quantile regression (Theorem 4.4). We have seen that regression quantile is not adapted for non Gaussian distributions. A second predictor is given in order to cope with regression quantile problems for extremal quantile levels. We have shown that the proposed extremal quantile predictor is equivalent to the true conditional quantile for extreme quantile levels. We have also illustrated on several numerical examples the better performance of this predictor for extreme levels. As a perspective, these prediction methods require the knowledge of the distribution of the covariates vector X 1 . We have not explored the prediction procedure when the X 1 's distribution is estimated (parametrically e.g.). This is an interesting perspective, which is let for future work. Finally, we would like to emphasize that we have given examples in dimension d = 1, but all the results may be used in higher dimensions (see Figure 4).

Figure 1 .

 1 Figure 1. Univariate mixture distributions, with µ = 0 and Σ = 1

  k φ(θ k t), where φ is the normal p.d.f, i.e φ(t) = e -θ)+r cos(θ)+1

Figure 2 .

 2 Figure 2. Theoretical quantiles, Quantile Regression and Extremal Predictors for Gaussian, Student, Gaussian Mixture and Uniform Mixture processes, with levels of quantile α = 0.995 and α = 0.005. In solid lines, the theoretical quantiles (approximated most of the time numerically), dotted, the Quantile Regression Predictor, and dashed, the Extremal Predictor.

Figure 3 .

 3 Figure 3. Q-Q plots for Student and Unimodal GM examples

Figure 4 .

 4 Figure 4. Quantile regression for a Student random field, with α = 0.0005, in dimension d = 2

Table 2 .

 2 Some classical consistent d-dimensional elliptical distributions with corresponding random variable

  Table 4 is a summary of the RMSE for all treated examples, and different levels of α. Obviously, we only consider the cases α ≥ 1 2 because elliptical distributions are symmetric.

	α	Gaussian	Student	Unimodal GM	Uniform GM
		r(qα) r( qα)	r(qα)	r( qα) r(qα)	r( qα)	r(qα)	r( qα)
	0.5	0	0	0	0	0	0	0	0
	0.6	0	0	0.027	0.372 0.001	0.019	0.067	0.334
	0.7	0	0	0.059	0.376 0.003	0.036	0.151	0.301
	0.8	0	0	0.104	0.358 0.006	0.051	0.297	0.249
	0.9	0	0	0.195	0.322 0.012	0.061	0.770	0.177
	0.95	0	0	0.314	0.289 0.024	0.061	1.880	0.123
	0.9995	0	0	2.880	0.148 0.163	0.000	250.172	0.020
	0.999995	0	0	16.546 0.081 0.109 7.130 • 10 -06 25 178.530 0.008


Benoit Anthony Leroy

Benoit Pauwels

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

After trimming the book binding of incoming theses, each volume is duplexscanned twice. Quality control is then performed with the original in hand. The paper volume is finally boxed for storage using a custom cloud application.

institutional repository

  our library initiated a mass-digitization project of all the Ph.D. theses produced at the university since its creation in 1834.! The goal of the project was twofold: improving the visibility and accessibility of our scientific research and freeing space for the development of our learning centers. Over 10000 volumes and 3 million pages needed to be digitized. We developed a custom-designed digitizing workflow Digitized theses are preserved in SAFE PLN Not to Outsource? That Is the Question. -extra human resources needed! -longer time to start! -taking full responsibility -total control of the production! -less expensive! -once in place, the workflow can be reused We decided to insource this mass digitization project. automaCc filenaming with QR code generated from our catalog duplex cross-checks ensure low error risk The submission representation of the thesis (TIFF files) is automatically transformed into an archival representation (JPEG 2000 files), from which various dissemination representations are generated (PDF, ALTO, txt). ! These representations are referenced in our institutional repository using a custom-designed METS profile. Quality Assurance is supported by custom software server Quality control is always performed with the original in hand. The client-server application allows the operator to report quality issues both in the digitized object and in the original paper object. ! client The SAFE Private LOCKSS Network is an international distributed preservation repository based on the LOCKSS technology operated by seven universities.Each collection of archives is automatically replicated in every box constituting the network, thus ensuring a broad geographical data redundancy. ! Automatic data integrity monitoring and comparison is performed regularly through the network via secured connections. Bit-rot and other forms of data corruption can thus be quickly detected, reported and corrected.


Christopher Badenhop

Jonathan Fuller

Joseph Hall

Benjamin Ramsey

Mason Rice

Evaluating ITU-T G.9959 Based Wireless Systems Used in Critical Infrastructure Assets

Keywords: ITU-T G.9959, Z-Wave, vulnerabilities, wireless sniffing

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The prevalence of wireless connectivity in industrial control and sensor systems is increasing because it extends communications ranges at lower cost than wired alternatives. A 2011 survey of industrial control system operators reported that wireless networks were deployed in 43% of industrial control systems and numerous additional deployments were projected [START_REF] Boyes | All quiet on the wireless front, Control Global, Schaumburg, Illinois[END_REF]. The survey respondents reported the use of IEEE 802.11, WirelessHART, Bluetooth and ZigBee as well as proprietary wireless systems [START_REF] Boyes | All quiet on the wireless front, Control Global, Schaumburg, Illinois[END_REF]. An analysis of vulnerabilities and suggestions for mitigation have been published for these protocols, with the exception of the proprietary systems [START_REF] Reaves | Analysis and mitigation of vulnerabilities in short-range wireless communications for industrial control systems[END_REF].

The growing use of proprietary systems makes it necessary to analyze and discuss their security implications, especially when considering these systems for use in critical infrastructure assets. Of the numerous proprietary systems, wireless systems based on the ITU-T G.9959 recommendation, which specifies a short range narrow-band digital radio communications transceiver operating in the sub-GHz spectrum, have significant potential for growth in the critical infrastructure. The most common commercial instantiation of ITU-T G.9959 is Z-Wave, which is standardized and marketed by the Z-Wave Alliance. The Z-Wave Alliance comprises more than 300 companies and is actively working to increase the adoption of Z-Wave products around the globe. IEEE 802.15.4 networks (e.g., WirelessHART and ZigBee) fulfill low-rate communications roles similar to those of Z-Wave in the critical infrastructure. Recent research efforts (see, e.g., [START_REF] Dubendorfer | ZigBee device verification for securing industrial control and building automation systems[END_REF][START_REF] Ramsey | Defensive rekeying strategies for physical-layermonitored low-rate wireless personal area networks[END_REF][START_REF] Ramsey | Wireless infrastructure protection using low-cost radio frequency fingerprinting receivers[END_REF]) have proposed novel security strategies for these networks. However, Z-Wave product development has been significantly restricted by nondisclosure and confidentiality agreements that stifle open-source security and resilience research. As a result, the security implications of the use of Z-Wave networks in the critical infrastructure are not well understood. In order to address this issue, this chapter introduces open-source techniques and tools for evaluating the security of ITU-T G.9959 wireless networks (including Z-Wave products) used in critical infrastructure assets.

ITU-T G.9959-Based Z-Wave Protocol

The ITU-T G.9959 recommendation specifies the physical (PHY) and media access control (MAC) layers for short-range, narrow-band digital radio communications transceivers. All manufactures adhere to the PHY/MAC specifications to ensure interoperability, but market their devices based on the network and application layers.

ITU-T G.9959-based networks operate in unlicensed frequency bands (e.g., 908.4 MHz in North America, 860.4 MHz in Europe and additional frequencies in other regions). Data rates of 9.6 Kbps (Rate 1), 40 Kbps (Rate 2) and 100 Kbps (Rate 3) are supported, depending on the transceiver type. Networks have two basic types of nodes: (i) control nodes; and (ii) end device nodes. Control nodes initiate commands while end device nodes respond to commands. Leveraging mesh topologies, end device nodes also forward commands to other nodes that are not directly reachable by a control node. The protocol allows a maximum of four hops between nodes and a maximum of 232 nodes in one network.

Each Z-Wave network is identified by a unique 32-bit Home ID, which is programmed by the manufacturer on each control node. The Home ID allows multiple networks to operate in close proximity without overlap. While a network may contain multiple control nodes, only one control node may be designated as the primary controller and its Home ID uniquely identifies the network. During the network inclusion process, the primary control node assigns the new node an 8-bit Node ID, which is only unique in the local network.

The Z-Wave protocol consists of four layers. Figure 1 shows the four-layer ITU-T G.9959/Z-Wave reference model mapped to the five-layer TCP/IP reference model. The PHY layer controls access to the radio frequency medium, the MAC layer handles transmission and reception of frames between adjacent nodes, the routing layer controls the flow of messages throughout the mesh and the application layer executes commands associated with the end device.

PHY Layer

The PHY layer uses carrier sense multiple access with collision avoidance to control access to the wireless medium. As mentioned above, Z-Wave utilizes unlicensed frequency bands, which differ according to the region. The protocol offers three data rates: (i) 9.6 Kbps using frequency-shift keying with Manchester encoding; (ii) 40 Kbps using frequency-shift keying with non-return-tozero encoding; and (iii) 100 Kbps using Gaussian frequency-shift keying with non-return to zero encoding.

As shown in Figure 2, the PHY protocol data unit (PPDU) consists of three main parts. The frame begins with a start header (SHR), which contains a preamble for symbol and bit synchronization, followed by a start of frame delimiter (SFD). The frame payload or PHY service data unit (PSDU) follows. Finally, for 9.6 Kbps data rate transmissions only, the frame concludes with an end header (EHR).

MAC Layer

The MAC layer (sometimes referred to as the transfer layer) is also detailed by ITU-T G.9959. The layer controls the transfer of data between two nodes and is responsible for frame acknowledgment, retransmission, data validation and notifying battery-operated devices to stay awake pending incoming transmissions. The MAC header contains the Home ID, Source ID, frame control field, frame length and Destination ID (or bitmask in the case of a multicast frame). The MAC footer contains either an 8-bit checksum or a 16-bit cyclic redundancy check depending on the data rate used.

Routing Layer

Z-Wave mesh network topologies are managed by the routing layer to ensure that messages are successfully routed among control and end device nodes. The protocol specifies a maximum of 232 nodes and only one primary control node, although multiple secondary control nodes may exist in order to partition a network. Every node, with the exception of battery-operated devices, participates in routing by forwarding frames between nodes outside of direct wireless transmission range. The protocol also specifies a maximum of four hops between the primary control node and any other node.

The routing layer is responsible for scanning the network topology and maintaining a routing table in the primary control node. The routing table is built by the primary control node based on information received, upon inclusion or request, from each end device node about the neighbors of each node. The Z-Wave protocol stack supports automatic topology discovery and healing to optimize routing tables when the location of a node has changed or a node has been removed from the network (exclusion).

Application Layer

The majority of the application layer is implementation-specific depending on the Z-Wave developer. For brevity, this section discusses only what is applicable to all Z-Wave devices.

The application layer frame consists of the header, command class information and command parameters (see Figure 2). The application layer is responsible for executing the commands passed to it. Commands are broken into two classes: command class and device class.

A command class is related to a specific function or device. An example is the binary switchcommand class. The binary switch uses three commands: (i) SET to turn a device on or off; (ii) GET to request the status of a device; and (iii) REPORT to respond to the request. These three commands are foundational to all Z-Wave devices.

The device class is subdivided into the basic, generic and special device classes. The basic device class distinguishes between controllers, end devices and end devices that are capable of routing. The generic device class defines the function that the device performs as a controller or end device. The special device class allows for more specificity in device functionality.

It is important to note that the Z-Wave protocol supports encryption using the Advanced Encryption Standard (AES) with 128-bit keys. When implemented, the application frame is encrypted and an 8-byte authentication frame header is appended to the end of the MAC service data unit. While data encryption is supported by the protocol, its implementation is left to the manufacturer to decide if the device transmission is sensitive enough to warrant encryption. In wireless sensor networks, memory and power are scarce resources, which discourages developers from implementing encryption "unless required" [START_REF]Freescale Beestack: Application Development Guide[END_REF]. Surveys of similar low-rate networks have demonstrated that the use of encryption or other security measures are far from ubiquitous [START_REF] Ramsey | Watching for weakness in wild WPANs[END_REF].

ITU-T G.9959/Z-Wave Attack Classes

While Z-Wave is a proprietary protocol, vendors may purchase software development kits to produce Z-Wave certified products. Given the nature of hardware and software development, certain vulnerabilities are introduced by developer-specific implementation faults as illustrated in [START_REF] Fouladi | Security evaluation of the Z-Wave wireless protocol[END_REF]. This section focuses on vulnerabilities of the underlying ITU-T G.9959 recommendation, which are common to all devices and attack classes that exploit the vulnerabilities.

Three classes of attacks are considered: (i) reconnaissance; (ii) denial-ofservice; and (iii) packet injection. The three classes of attacks undermine network confidentiality, availability and integrity, respectively. A reconnaissance attack involves the passive collection of traffic or the active probing of a target network to gain information without interfering with normal operations. A denial-of-service attack prevents wireless system access and causes varying degrees of system unavailability. A packet injection attack involves the transmission of specially-crafted packets to manipulate network or device behavior.

Reconnaissance Attacks

Reconnaissance lays the foundation for sophisticated follow-up attacks. The information acquired includes the protocols in use, device types, traffic flow patterns and even encryption keys if they are not handled properly. Information received from reconnaissance can help an attacker obtain accurate mappings of a system, services and/or vulnerabilities, enabling more significant attacks to be conducted in the future. Using a high-gain antenna, observations can be made at long distances, allowing an attacker to remain inconspicuous while gathering information. Apa and Hollman [START_REF] Apa | Compromising industrial facilities from 40 miles away, presented at the Black Hat USA Conference[END_REF] have presented a proven exploitation of a wireless sensor network used in critical infrastructure assets. In particular, they demonstrated the exploitation of three devices from a maximum distance of 64 km.

An attacker armed with a directional antenna can capture ITU-T G.9959 transmissions for further analysis. Information gathered that might be useful to an attacker includes: (i) traffic patterns; (ii) use of encryption during transmission; and (iii) frame header content, which includes the unique Home ID of an ITU-T G.9959 network, Source ID of the device being sniffed and Destination ID.

Two demonstrations of ITU-T G.9959 exploitation have been published to date, both of them relied heavily on the ability to conduct reconnaissance to gain the knowledge required to craft follow-up attacks. Fouladi and Ghanoun [START_REF] Fouladi | Security evaluation of the Z-Wave wireless protocol[END_REF] showed how to obtain a detailed understanding of the manner in which a secure door lock implements encryption and authentication; using this knowledge, they were able to discover a flaw that could be exploited. Picod et al. [START_REF] Picod | Bringing software defined radio to the penetration testing community[END_REF] were able to discern the specific commands and associated bit values that could be used to turn an alarm on and off. Follow-on attacks of these example attacks are discussed later in this chapter.

Several entities have developed sniffers for intercepting and transmitting Z-Wave frames. Sigma Designs [START_REF]Z-Wave Products[END_REF] markets a closed-source development kit for Z-Wave device developers. The kit incorporates several hardware development platforms, technical documentation, software tools and a Z-Wave protocol sniffer. However, the kit comes with a non-disclosure agreement, which restricts the use of the tools.

A second sniffer project is z-force [START_REF] Fouladi | Security evaluation of the Z-Wave wireless protocol[END_REF]. The sniffer includes custom firmware hosted on a CC1110 development board [START_REF]Texas Instruments, Low-Power Sub-1 GHz RF Transceiver, CC1101, SWRS0611[END_REF] and a z-force personal computer application. At the time of this writing, z-force has only been demonstrated on European Z-Wave frequency bands. The closed-source nature of z-force makes it difficult to evaluate and extend.

A third Z-Wave sniffer is Scapy-Radio, part of the open-source hackrf project [START_REF] Picod | Bringing software defined radio to the penetration testing community[END_REF]. Scapy-Radio integrates Scapy, a Python environment for manipulating network traffic and GNU Radio, an open-source signal processing toolbox. The tool includes GNU Radio companion implementations of Z-Wave, Bluetooth and ZigBee transceivers. The Z-Wave implementation receives samples from a software-defined radio, demodulates, synchronizes and decodes NRZ or Manchester encodings. As with z-force, the implementation is tuned to European bands, but can be modified for other bands using GNU Radio tools.

Application layer sniffers are available from OpenZWave.

OpenZWave is an open-source project that provides libraries and drivers to communicate with USB-based Z-Wave controllers (e.g., Z-Stick S2 from Aeon Labs). OpenZWave provides a network querying tool to demonstrate its MinOZW library API. The python-openzwave package, an open-source wrapper for OpenZWave, provides a Python version of MinOZW and a Z-Wave shell interface to interact with a Z-Wave controller.

While Scapy or a netcat listener could be used to collect the encapsulated UDP frames sent over a local loopback device, future work will support the capture and analysis of Z-Wave frames using Wireshark. Wireshark is capable of intercepting UDP datagrams, but is unable to decode the encapsulation header, MAC header or Z-Wave payloads.

As part of this research, a packet dissector was developed for the encapsulation header and Z-Wave MAC header. Figure 3 shows a Z-Wave frame dissection. The dissector also decodes the first byte of the PHY service data unit payload to identify the application payload command class. The remaining payload bytes depend on the value of the command class field. The dissector computes a checksum for each frame and denotes the outcome to the user in the info field of the packet list window. Dissectors for each of the 49 command classes, as identified in the OpenZWave source code, are currently being developed.

Denial-of-Service Attacks

Denial-of-service attacks prevent or degrade legitimate access to resources (e.g., radio frequency spectrum). Such attacks can be easily accomplished using off-the-shelf equipment [START_REF] Pelechrinis | Denial-of-service attacks on wireless networks: The case of the jammers[END_REF]. Denial-of-service attacks on wireless networks can deny or degrade access to resources from the physical layer to the network layer. For example, the physical layer may be susceptible to narrow-band jamming, the carrier sense algorithms may be exploited to deny access to the medium [START_REF] Xu | The feasibility of launching and detecting jamming attacks in wireless networks[END_REF] and routers can be consumed with overflowing interface queues.

Constant and deceptive jamming are effective for conducting denial-of-service attacks due to the MAC layer collision avoidance characteristics of ITU-T G.9959. When an attacker is operating a constant or deceptive jammer, any node within range will sense the channel as busy and wait to transmit. Even more effective and efficient is reactive jamming [START_REF] Li | Optimal jamming attack strategies and network defense policies in wireless sensor networks[END_REF], which is difficult to detect [START_REF] Strasser | Detection of reactive jamming in sensor networks[END_REF]. A reactive jammer, that only transmits after the preamble and start of a frame delimiter of an ITU-T G.9959 PHY frame are detected, merely has to corrupt one bit of the PHY service data unit in order to cause an integrity check error and the complete loss of the frame. The non-correcting integrity checks used by Z-Wave are capable of detecting, but not correcting, single-bit errors. Even worse, the corruption of a single bit in the Z-Wave PHY layer, unlike a PHY layer that uses a spreading technique such as direct-sequence spread spectrum (DSSS), is achievable using narrow-band jamming. The use of error correction codes, as in IEEE 802.11a, is more robust to bitwise jamming [START_REF] Noubir | Robust wireless infrastructure against jamming attacks, in Handbook on Securing Cyber-Physical Critical Infrastructure: Foundations and Challenges[END_REF]. Table 1 presents an estimate of jammer efficiency in terms of bits jammed per bit transmitted against ITU-T G.9959 (i.e., ratio of communications effort to jammer effort) based on the results in [START_REF] Noubir | Robust wireless infrastructure against jamming attacks, in Handbook on Securing Cyber-Physical Critical Infrastructure: Foundations and Challenges[END_REF].

Depending on the objective, an attacker may use any of the methods described above to impact the availability of one or more nodes in a wireless network. For example, a Z-Wave network containing a thermostat (sensor) and water valve (actuator) could be subjected to a denial-of-service attack that prevents the thermostat from reporting the current temperature or obstructs a command to activate the water valve. A similar scenario was successfully demonstrated in [START_REF] Reaves | Analysis and mitigation of vulnerabilities in short-range wireless communications for industrial control systems[END_REF] with a gas pipeline remote terminal unit that included a wireless pressure sensor, pump and relief valve.

Packet Injection Attacks

Due to the broadcast nature of wireless networks, an attacker armed with information gained via reconnaissance may be able to inject forged packets into a network. The ability to conduct packet injection enables the attacker to masquerade as a legitimate network device while transmitting messages to manipulate system operation. Badenhop and Mullins [START_REF] Badenhop | A black hole attack model using topology approximation for reactive ad-hoc routing protocols[END_REF] have investigated network degradation attacks against wireless routing protocols similar to those used in Z-Wave.

Using publicly-available hardware and software, researchers have reported the ability to conduct packet injection attacks to manipulate ITU-T G.9959 devices. Fouladi and Ghanoun [START_REF] Fouladi | Security evaluation of the Z-Wave wireless protocol[END_REF] have developed a packet inspection/injection tool using a Texas Instruments radio transceiver and custom software (now publicly-available, albeit not open source). Using the tool, Fouladi and Ghanoun were able to inject traffic to exploit a vulnerability in a device-specific encryption implementation at the application layer that enabled them to send encrypted commands to perform unauthorized actions.

Picod et al. [START_REF] Picod | Bringing software defined radio to the penetration testing community[END_REF] have also demonstrated a packet injection attack against ITU-T G.9959 devices. Their packet inspection/injection tool uses a softwaredefined radio and open-source software packages. They demonstrated the ability to sense legitimate SWITCH BINARY ON commands and automatically inject subsequent SWITCH BINARY OFF commands, which effectively nullified legitimate messages.

These preliminary case studies demonstrate the attack possibilities when a malicious entity has the ability to inject forged ITU-T G.9959 traffic. In addition to these examples, an attacker may be able to: (i) flood a network with traffic causing a denial of service at the routing or application layers; (ii) send false status messages; or (iii) provide a control node with false routing information to poison the network routing table.

If the target ITU-T G.9959-based network is used in a critical infrastructure asset, the ability to inject commands and report false state information could prove disastrous. For example, an attacker may choose to send a close command to a water valve in a cooling system immediately following every open command sent by the master terminal unit while falsely reporting to the human-machine interface that the valve is open.

Passive Reconnaissance Techniques

Several high-level attack classes for Z-Wave systems have been presented above. A basic, invariant system threat from which other threats originate is passive reconnaissance. In the remainder of this chapter, the focus is on identifying several low-complexity passive reconnaissance techniques to further the exploration of Z-Wave security.

To capture arbitrary Z-Wave frames, a receiver chain was constructed using a USRP N210, GNU Radio, Scapy-Radio and a Z-Wave Wireshark packet dissector. Figure 4 shows the sniffer architecture. The USRP was tuned to a center frequency of 908.40MHz to collect frequency-shift-keyed symbols at ±20 KHz deviation from the center frequency. An Ettus VERT900 dipole antenna was connected directly to the USRP and sampling was performed at 800 Kbps with 20 MHz bandwidth. Filtering, demodulation and symbol synchronization were performed by GNU Radio. Scapy-Radio, a subsystem within GNU Radio, provided preamble detection, byte synchronization, NRZ decoding and frame extraction operations. The extracted frames were transmitted over the localhost interface encapsulated as UDP datagrams to port 52002. Wireshark captured traffic over the localhost interface to intercept the datagrams. A custom Z-Wave packet dissector was used to decode the encapsulated Z-Wave frames, which were preserved in the PCAP format for later analysis. The experimental setup in Figure 4 permitted the passive collection of Z-Wave frames. Unlike application level sniffers, Z-Wave frames were captured regardless of the frame Home-ID. Moreover, the PHY and MAC layers were retained rather than stripped in order to provide insights into system behavior.

Controller-Device Pairing

Observations of the responses of the included device to the controller revealed significant details about the hardware, software and current state. At the end of an inclusion process, the controller interrogates the included device to learn about it. This information is used by the controller to present accurate control options and device information to the user.

Figure 5 shows a frame capture of an Aeon Labs Z-Stick controller using a Jasco JS-45603 smart dimmer. The data in the application payloads, which was reversed by referencing the Open-ZWave command class source code, is summarized to the right of Figure 5. The figure shows a sequence of GET commands from the controller to the dimmer. The dimmer responds to each request with a REPORT message that reveals information. Device-level information, such as the vendor, device type and ID are provided in a manufacturer-specific REPORT message. The controller requests software version information, to which the dimmer device replies with a Z-Wave library, protocol and application version. After the device and software level information is gathered, the controller queries the device on the current state of each configuration and control item. This state information can be used by an observer to determine the set of messages that the device services, which is especially useful when the device-specific data refers to an unknown product type or manufacturer ID. The location state is of particular interest. While location information was not initialized in the capture shown in Figure 5, it is still possible for a device to reveal some information about its configured location. The location information is in the form of a grouping index that corresponds to a user-defined label assigned at the controller interface. A user may group devices by room, floor or functionality, so additional information is required to understand the semantics of this value.

If the pairing is not observable, passive reconnaissance may still capture the interactions between the controller and end devices. As with the latter portion of the pairing process, the command classes of observed messages provide information about the device function. The drawback is that the observation period is significantly longer than if observed during the pairing process. This is especially true for low-power event triggered devices. These devices operate at low duty cycles to extend battery life and only emit frames when a specific monitored event occurs (e.g., a door opening or an excessive sensor reading).

Controller-Specific Behavior

A controller may exhibit discernible differences in how it interacts with devices. To illustrate this observation, the dimmer (multilevel switch class device) was paired with two different controllers to show that the frames are noticeably different.

In the first experiment, the Jasco Smart Dimmer was paired with an Aeon Labs Z-Stick S2 controller (Node 1). The multilevel switch SET command was executed through an OpenZWave control panel while the Z-Wave frames were captured by the Scapy-Radio USRP sniffer and processed by the Wireshark Z-Wave dissector. Note that the Z-Stick controller has an ID of 1 and the dimmer has a node ID of 5. The Z-Stick controller sends a multilevel dimmer command to set the dimmer to a value, shown as frame 52, which is acknowledged by the device. After receiving the acknowledgment (ACK), the controller requests an update of the current dimmer state (output voltage) by sending frame number 56 in Figure 6. The dimmer complies by sending its current value to the controller in frame 60 in Figure 6. This results in six frames being exchanged per transaction. The same pattern occurs when the Aeon Z-Stick controller sends commands to other devices. Essentially, the controller issues a SET command, followed by a GET command on the same value to confirm that the device complied with the command. The targeted device responds to the command with a REPORT command. In the second experiment, the same dimmer device was paired with a Ve-raLite controller. Using the VeraLite controller web interface, the multilevel switch SET command was invoked while capturing Z-Wave frames. Figure 7 shows the results of the capture. Coincidentally, the node IDs are once again one and five. The VeraLite controller sends two versions of the same multi-level switch SET command. In frame 63, the ACK-required bit in the frame control field is not set; however, it is set in frame 65. After the dimmer device receives frame 65, it acknowledges receipt to the controller, which is shown in Figure 7 by the matching sequence numbers. The controller initiates the same command with and without the ACK-required bit set is very likely because the controller is compatible with devices that are not able to reply. Devices that are unable to transmit may ignore frames that require an ACK. As with the Z-Stick controller, this pattern is observable when the VeraLite controller issues SET commands to devices. The SET command is sent both with and without the ACK and does not necessarily expect a reply.

The results show that there is more than one way to execute a switch SET transaction. While it cannot be concluded that the observed patterns are unique to the controllers, it can be deduced that a given command pattern reduces the possible identities of the target controller. For example, when an unknown controller exhibits a pattern consistent with that of a Z-Stick controller, this would imply that the unknown device is not a VeraLite controller.

Device-Specific ACK Times

When the application layer of Z-Wave frames are encrypted, passive reconnaissance of a device using the previously-identified methods may prove ineffective. Application layer encryption obscures the command class fields, making it more difficult to identify controller-specific command patterns. An outside observer must turn to non-encrypted observations to deduce properties of the devices. Implementation differences between hardware, software and the physical environments of devices may result in observable differences outside of the encrypted traffic. Existing techniques such as traffic analysis and side channel analysis have been shown to be effective at thwarting confidentiality mechanisms like encryption [START_REF] Wright | Uncovering spoken phrases in encrypted voice-over-IP conversations[END_REF].

A fingerprinting technique known as preamble manipulation may be applicable to Z-Wave devices [START_REF] Ramsey | Wireless infrastructure protection using low-cost radio frequency fingerprinting receivers[END_REF]. Unfortunately, it was difficult to find a Z-Wave device that uses a transceiver other than the ZM3102N to verify this technique. Since the experimental hardware was homogeneous, firmware implementations could be examined behaviorally to identify implementation and functionality differences. As a matter of fact, the 32 KB of flash memory available in the ZM3102N provides ample opportunities for vendors to customize their implementations.

One ubiquitous observable behavior below the application layer is the time taken by a node to send an ACK upon receiving a request. Time-of-flight fingerprinting metrics are explored in [START_REF] Schauer | Potentials and limitations of WiFi-positioning using time-of-flight[END_REF], where the authors report that measurement variance is due to hardware and environmental factors. In their experiments, laptop computers were used to perform packet injection, ACK response and time measurements. Each had an operating system that managed independent operations, but also competed for CPU clock cycles, resulting in variations in the time-of-flight measurements. While this is a valid issue, Z-Wave devices are far more specialized than laptop computers, having application-specific I/O and fewer, if any, concurrent tasks. This suggests that an ACK response of a Z-Wave device is more deterministic than a general purpose laptop. If this is true, ACK-based fingerprinting would be a valid approach for identifying Z-Wave devices.

The experiments described in the remainder of this section demonstrate that ACK-based fingerprinting of Z-Wave devices shows promise. While the Z-Wave devices used in the experiments were all embedded systems, the sniffer ran on a Dell Latitude E6520 laptop. The laptop incurred a degree of the dynamics identified in [START_REF] Schauer | Potentials and limitations of WiFi-positioning using time-of-flight[END_REF]. The ACK delay time was measured using the packet arrival times observed by the Wireshark Z-Wave dissector. The time values were determined by the Wireshark process when an encapsulated Z-Wave frame was observed on the localhost interface. Encapsulated frames were sent over the localhost interface by the GNU Radio process. The GNU Radio process received samples from the USRP FPGA via a gigabit Ethernet interface. In the experiments, each process, interface and device driver contended with other resources on the sniffer laptop and were subject to contention dynamics. Regardless of the dynamics, adequate sampling of ACK delays eliminated the confounding effects from the measurements. Seventy ACK response times were collected for each device under test to mitigate sampling bias and to facilitate observational comparisons.

ACK Times for Different Vendor Devices. In the first experiment, a VeraLite controller was used to repeatedly issue commands to a FortrezZ water valve. Each command generated several acknowledged frames from the controller and valve. The ACK response to a particular singlecast frame could be identified by the matching four-bit sequence numbers. After carefully pairing the singlecast frames with their associated ACKs, the response time was measured as the difference in time of between the observation of the singlecast frame being sent and the observation of the ACK being sent. The distances between the water valve, VeraLite controller and the Z-Wave sniffer were all within one meter of one another, implying that signal propagation delays were negligible.

Figure 8 shows the 99% confidence intervals estimating the true mean (n = 70) of the ACK response times for the VeraLite controller and FortrezZ water valve. The mean response times are different with a two-sample t-test p-value of 0.00; clearly, the response time of the water valve is greater than that of the VeraLite controller. Moreover, the difference in the response times is an order of magnitude more than variations due to thermal effects [START_REF] Zensys | N Z-Wave Module Datasheet[END_REF]. Since both devices use the same hardware, the differences may be due to firmware implementation differences.

ACK Times for Same Vendor Devices. The same experiment was repeated, except that the selected controller and device originated from the same vendor. The Aeon Labs Z-Stick S2 controller was used as the controller and an Aeon Labs Appliance Switch was used as the device. A total of 70 ACKs were recorded for both systems by repeatedly issuing commands from the controller. Figure 9 shows the observed mean ACK response times for the two devices (n = 70 and a confidence interval of 99%). In this case, the mean ACK response times of the two Aeon Labs devices are not statistically different (two-sample t-test p-value = 0.64). This suggests that both devices have similar ZM3102N implementations.

ACK Times for Polling Commands. In this experiment, instead of manually creating Z-Wave traffic, the VeraLite controller was configured to automatically poll two devices every 60 seconds. The two devices polled were the FortrezZ WV01LFUS075 water valve and an Everspring AN145 lamp socket switch. Figure 10 shows the mean response delays for each device with 99% confidence intervals (n = 70). The figure reveals several interesting points. First, there is a statistically significant difference between the response times of the VeraLite controller and the other two devices (ANOVA p-value = 0.00). The difference of means between the water valve and lamp switch is not statistically significant, but the variance is higher for the Everspring switch. Second, the mean ACK response time of the water valve is different from that seen in Figure 8. The variances of the devices in Figure 10 are larger than those reported in Figures 8 and9. This experiment had 60 seconds between each transaction, whereas the two previous experiments had repetitions at intervals of one or two seconds. During the 60-second period of inactivity, it is hypothesized that the devices went into a lower power state, which increased their ACK response times.

Differences in Mean ACK Response Times. The experimental data suggests that the differences in means are due to implementation differences in how ACKs are handled by the devices. Upon receiving a frame, a checksum is performed. If the checksum passes and the frame has either requested an ACK or the protocol requires it, the recipient generates an ACK using the sequence number of the received frame and sends the response to the originator. The differences arise depending on when these steps are taken during the transceiver chains. The ACK response may be initiated immediately upon receiving the frame or after the frame is queued so that the receiver may quickly return to search for the next frame. This option, for example, is provided by the CC2420 transceiver [19].

After the ACK is generated, another implementation decision involves the transmission of the ACK frame. If the ACK is given priority, it is either placed at the front of the send queue or an interrupt is generated to force the transceiver to immediately transmit the ACK. Other design choices include the sizes of the receive and send queues, buffer exception handling and the queue service rates. In a delay-prone scenario, ACKs are handled by an attached microprocessor via the SPI bus.

Another source of differences may be the ZM3102N configuration settings. The ITU-T G.9959 specification lists parameters relevant to ACK response times such as MacMinAckWaitDuration, TurnaroundTimeRXTX and MacMin CCARetryDuration. The three parameters specify the minimum frame spacing between receiving a packet and transmitting its ACK, the time penalty for switching from the receive mode to the transmit mode and the minimum time to wait between clear channel assessments, respectively. It is not clear which, if any, of these parameters are configurable on the ZM3102N; however, each may impact the ACK response time.

Conclusions

This chapter has examined the security implications of using ITU-T G.9959 wireless networks with Z-Wave devices in critical infrastructure assets. Several techniques for passively discriminating between Z-Wave devices based on functionality and vendor were investigated. Experiments involving passive observations of ACK response times demonstrate that it is possible to identify implementation differences in ZM3102N firmware. In particular, the experimental results reveal that a VeraLite controller and FortrezZ water valve have different mean ACK response times, while an Z-Stick controller and appliance switch from the same vendor (Aeon Labs) have equivalent response times. These results suggest intra-vendor similarities.

Future research will focus on developing new passive techniques and refining existing techniques for use with other types of devices. Additionally, research will attempt to understand the reasons for the differences in the mean ACK response times. Efforts will also explore active fingerprinting techniques. One such technique will involve sending messages corresponding to all the command classes to a target device to discern its capabilities. The ultimate goal is to develop a tool that leverages passive and active techniques in device fingerprinting, with functionality similar to the popular nmap operating system fingerprinting tool.

Figure 1 .

 1 Figure 1. ITU-T G.9959/Z-Wave model mapped to the TCP/IP reference model.

Figure 2 .

 2 Figure 2. ITU-T G.9959 and Z-Wave frame formats.

Figure 3 .

 3 Figure 3. Wireshark Z-Wave dissection.

Figure 4 .

 4 Figure 4. Sniffer architecture.

Figure 5 .

 5 Figure 5. Pairing between a Jasco smart dimmer and an Aeon Z-Stick controller.

Figure 6 .

 6 Figure 6. Z-Stick controller SET command capture.

Figure 7 .

 7 Figure 7. VeraLite controller SET command capture.

Figure 8 .

 8 Figure 8. Mean ACK response times for devices from two different vendors.

Figure 9 .

 9 Figure 9. Mean ACK response times for devices from the same vendor.

Figure 10 .

 10 Figure 10. Mean ACK response times for a 60-second polling interval.

Table 1 .

 1 Jamming efficiency against ITU-T G.9959 frames.

	Integrity	Check*	Max	Jammer	Efficiency
	Data Rate		Payload Bits to Jam
	9.6 Kbps (Rate 1)	8-bit checksum	512 bits	1	512
	40 Kbps (Rate 2)	8-bit checksum	512 bits	1	512
	100 Kbps (Rate 3) 16-bit CRC	1,360 bits	1	1,360
	*ITU-T G.9959 uses non-correcting integrity checks.

Note that the views expressed in this chapter are those of the authors and do not reflect the official policy or position of the U.S. Air Force, U.S. Army, U.S. Department of Defense or U.S. Government.