Laurent Gardes
email: gardes@unistra.fr

Gilles Stupfler

Estimating extreme quantiles under random truncation

Keywords: Asymptotic normality, consistency, extreme quantile, heavy-tailed distribution, tail index. AMS Subject Classifications: 62G05, 62G20, 62G30, 62G32

The goal of this paper is to provide estimators of the tail index and extreme quantiles of a heavy-tailed random variable when it is righttruncated. The weak consistency and asymptotic normality of the estimators are established. The finite sample performance of our estimators is illustrated on a simulation study and we showcase our estimators on a real set of failure data.

Introduction

Studying extreme events is relevant in numerous fields of statistical applications. One can think about hydrology, where one may want to estimate the maximum level reached by seawater along a coast over a given period, or to study extreme rainfall at a given location; in actuarial science, a pivotal problem for an insurance firm is to estimate the probability that a claim so large that it represents a threat to its solvency is filed. In this type of problem, the focus is not in the estimation of "central" parameters of the random variable of interest, such as its mean or median, but rather in the understanding of its behavior in its right tail. A particular relevant case is when the random variable of interest, Y , is heavytailed, namely, when its survival function F can be written F (y) = y -1/γ L(y) for all y > 0; here, γ > 0 shall be referred to as the tail index and L is a slowly varying function at infinity, meaning that it satisfies L(λy)/L(y) → 1 as y → ∞ for all λ > 0. In this case, γ clearly drives the tail behavior of F and its knowledge is necessary if, for instance, we are interested in the estimation of extreme quantiles of Y . The estimation of the tail index is thus one of the central topics in extreme value theory, which is why this problem has been extensively studied in the literature. Recent overviews on univariate tail index estimation can be found in the monographs of [START_REF] Beirlant | Statistics of Extremes[END_REF] and de [START_REF] De Haan | Extreme value theory: An introduction[END_REF]. A further challenge arises when facing incomplete data. An example of such a situation is the estimation of (extreme) survival times based on a follow-up study of patients suffering from a given illness. If at the time the data are collected a patient is still alive, then his/her survival time is not available to the researcher, although it is known that the patient survived until the end of the study. This case is the archetypal example of right-censoring. Estimating the tail index in this situation is much more difficult than when having complete data, since information about the right tail of the variable of interest is missing. In this setting, the estimation of the tail index and extreme quantiles has been considered by [START_REF] Beirlant | Pareto index estimation under moderate right censoring[END_REF], [START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF], [START_REF] Beirlant | Peaks-over-threshold modeling under random censoring[END_REF], [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF], [START_REF] Gomes | Estimation of the extreme value index for randomly censored data[END_REF] and [START_REF] Worms | New estimators of the extreme value index under random right censoring, for heavy-tailed distributions[END_REF]. In this paper, we consider the case when the data are right-truncated. In this framework, one observes the variable of interest if and only if it is less than or equal to a truncation variable T . This situation is different from right-censoring since nothing is known about Y in the case Y > T , which adds a further difficulty to the analysis of the right tail of Y . Truncated data may be collected in various cases, for instance when estimating incubation times for a given disease, see Kalbfleisch andLawless (1989, 1991) and [START_REF] Lagakos | Nonparametric analysis of truncated survival data, with application to AIDS[END_REF]; when studying the luminosity of astronomical objects such as quasars, see [START_REF] Jackson | The analysis of quasar samples[END_REF] and [START_REF]A method of allowing for known observational selection in small samples applied to 3CR quasars[END_REF]; when accounting for reporting lags in insurance data, also referred to as the incurred but not yet reported problem, see [START_REF] Herbst | An application of randomly truncated data models in reserving IBNR claims[END_REF], [START_REF] Kaminsky | Prediction of IBNR claim counts by modelling the distribution of report lags[END_REF] and [START_REF] Lawless | Adjustments for reporting delays and the prediction of occurred but not reported events[END_REF]; or when considering failure or warranty data, see Hu andLawless (1996a, 1996b) and the monographs by [START_REF] Meeker | Statistical methods for reliability data[END_REF] and [START_REF] Lawless | Statistical models and methods for lifetime data[END_REF]. To the best of our knowledge, the estimation of the tail index and extreme quantiles in this context is, up to now, still an open question. The outline of this paper is as follows. In Section 2, we give a precise definition of our model and define our estimators of the tail index and of the extreme quantiles of a truncated random variable. Some asymptotic properties of our estimators are stated in Section 3. The finite sample performance of the extreme quantile estimator is studied in Section 4. A real set of brake failure data is analyzed in Section 5. We offer some concluding remarks in Section 6. Proofs of the main results are deferred to Section 7, while the preliminary results are deferred to the Appendix. Proofs of the preliminary results can be found in a supplementary material.

Framework

Let (Y 1 , T 1 ), . . . , (Y n , T n ) be n independent copies of a random pair (Y, T ) ∈ [y 0 , ∞) × [t 0 , ∞), where Y and T are independent and y 0 , t 0 ≥ 0 are the left endpoints of Y and T . The right endpoints of Y and T are supposed to be infinite. The joint cumulative distribution function (cdf) of the random pair (Y, T ) is then given for all (y, t) ∈ R 2 by H(y, t) := P(Y ≤ y, T ≤ t) = F (y)G(t), where F and G are the cdfs of Y and T . The focus of this paper is on extreme quantiles of Y and, as a first step, on the estimation of the cdf F . Of course, because we only record the Y i and T i such that Y i ≤ T i , the classical nonparametric estimator of F cannot be used. However, the conditional cdfs of Y and T given Y ≤ T , respectively denoted by F * and G * , may be estimated in a nonparametric way. Let N be the total (random) number of observed pairs (Y i , T i ) such that Y i ≤ T i and notice that N is a binomial random variable with parameters n and p := P(Y ≤ T ). Such pairs shall be denoted in the sequel as

(Y * i , T * i ), 1 ≤ i ≤ N . It can be shown (see Lemma 1) that the conditional distribution of {(Y * i , T * i ), i = 1, . . . , N } given N = m is
equal to the distribution of m independent copies of a random vector (Y * , T * ) with cdf H * given by H * (y, t) = P(Y ≤ y, T ≤ t|Y ≤ T ). The standard estimators of the conditional cdfs of Y and T are then

F * N (y) = 1 N N i=1 I {Y * i ≤y} and G * N (t) = 1 N N i=1 I {T * i ≤t} .
Note now that, in order to estimate F , it is sufficient to estimate the function Λ F :=log F . The following result, whose proof can be found in Woodroofe (1985, p.166), shows that this quantity is in fact linked to F * and G * :

Proposition 1 Let C * := F * -G * . Then C * (y) = p -1 F (y)(1 -G(y))
> 0 for all y > y 0 , and

Λ F (y) = ∞ y dF (z) F (z) = ∞ y dF * (z) C * (z) .
This result can be used to build an estimator of the function Λ F and consequently of the cdf F : if y > y 0 , we may estimate Λ F (y) by

Λ F N (y) = 1 N N i=1 I {Y * i >y} C * N (Y * i )
when N > 0 and 0 otherwise. The survival function F := 1 -F and its associated quantile function α → q(α) := inf{y ≥ y 0 | F (y) ≤ α}, which is the right-continuous inverse of F , are then estimated by F N (y) = exp(-Λ F N (y)) and q N (α) = inf{y ≥ y 0 | F N (y) ≤ α}, where we let F N = 1 -F N . The first aim of this paper is to study the asymptotic behavior of the estimator q N (α n ) where α n → 0 as n → ∞. We shall tackle this problem in a framework of regular variation: we write that a function Ψ ∈ RV 1 (a), a ∈ R, if Ψ is nonnegative on (0, ∞) and for all λ > 0, we have Ψ(λy)/Ψ(y) → λ a as y → ∞. We thus consider the following model:

(M) We have F ∈ RV 1 (-1/γ F ) and G ∈ RV 1 (-1/γ G ), with 0 < γ F ≤ γ G .
Model (M) is a standard extreme-value model adapted to right-truncated data; see also [START_REF] Beirlant | Estimation of the extreme value index and extreme quantiles under random censoring[END_REF] and [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF] for closely related models when there is right-censoring. The cdfs of Y and T are thus heavytailed with respective tail indices γ F and γ G . The condition γ F ≤ γ G ensures that we have at our disposal enough observations pertaining to the right tail of Y . In this context, the quantile estimator q N (α n ) is consistent if α n → 0 slowly enough, see Theorem 2. In order to remove the restriction on the rate of convergence of α n , we note that under model (M), the quantile function q is regularly varying at 0 (see Corollary 1.2.10 p.23 in de Haan and Ferreira 2006), so that if β n < α n are two positive sequences tending to 0 such that

β n /α n → 0 then q(β n ) ≈ q(α n ) (α n /β n )
γ F when n is large. In order to derive an estimator of an extreme quantile q(β n ) from that, we first need to build an estimator of γ F . With this aim in mind, we remark that F * and G * are heavy-tailed with respective tail indices γ F * := γ F γ G /(γ F + γ G ) and γ G (see Lemma 3) and we introduce the Hill-type estimators (see [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF])

γ N,F * (k N ) = 1 k N k N i=1 log Y * N -i+1,N Y * N -k N ,N and γ N,G (k ′ N ) = 1 k ′ N k ′ N i=1 log T * N -i+1,N T * N -k ′ N ,N
.

Here we let, given

N = m, k N = k m and k ′ N = k ′ m ,
where k m and k ′ m are integers which belong to {1, . . . , m -1}, and

Y * 1,N ≤ . . . ≤ Y * N,N , T * 1,N ≤ . . . ≤ T * N,N
are the order statistics deduced from the samples (Y * i ) 1≤i≤N , (T * i ) 1≤i≤N . It is fairly easy to prove (see Lemma 9) that γ N,F * (k N ) and γ N,G (k ′ N ) are consistent estimators of γ F * and γ G under mild conditions. This leads us to introduce the class of estimators

γ N,F (k N , k ′ N ) = γ N,F * (k N ) γ N,G (k ′ N ) γ N,G (k ′ N ) -γ N,F * (k N ) . (1)
Under some conditions on (k m ) and (k ′ m ), the quantity γ N,F (k N , k ′ N ) is then a consistent estimator of γ F , see Theorem 3. This motivates the following Weissman-type estimator (see [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF] for a quantile having arbitrary order β n → 0:

q W N (β n | α n , k N , k ′ N ) = q N (α n ) (α n /β n ) γ N,F (k N ,k ′ N ) (2)
where α n → 0 converges slowly enough. The asymptotic properties of this estimator are discussed in Theorem 4.

Main results

In this section, we examine the asymptotic properties of our estimators. In order to establish the asymptotic normality of F N (y n ), we introduce the following additional condition:

∞ y0 dF (z) G(z) < ∞. ( 3
)
This assumption is classical in the study of the estimator of the cdf of a truncated random variable, see for instance [START_REF] Stute | The central limit theorem under random truncation[END_REF] and [START_REF] Woodroofe | Estimating a distribution function with truncated data[END_REF] for related hypotheses when there is left-truncation. Note that under model (M), it is a consequence of Lemma 2 with ϕ = 1/G and ψ = F that (3) automatically holds if γ F < γ G . Besides, it is easy to check that (3) fails to hold if γ F > γ G .

Theorem 1 Let y n → ∞. Assume that (M) and (3) hold, and that nv 2 (y n ) converges to infinity where

v(y) := F (y) ∞ y dF (z) G(z) -1/2 . Then v(y n ) √ n F N (y n ) F (y n ) -1 = ξ n if γ F < γ G , O P (1) if γ F = γ G ,
where ξ n is a random variable which is asymptotically standard Gaussian distributed.

We now establish the asymptotic normality of q N (α n ).

Theorem 2 Let α n → 0. Assume that F is a differentiable function in a neighborhood of infinity such that yF ′ (y)/F (y) → 1/γ F as y → ∞, that (M) and (3) hold, and that nv 2 (q(α n )) → ∞. Then

v(q(α n )) √ n q N (α n ) q(α n ) -1 = ζ n if γ F < γ G , O P (1) if γ F = γ G ,
where ζ n is a random variable which is asymptotically Gaussian centered with variance γ 2 F .

Theorem 2 is a convergence result for the intermediate quantile estimator q N (α n ), provided nv 2 (q(α n )) → ∞, which ensures that α n → 0 slowly enough. To examine the asymptotic properties of the extreme quantile estimator (2), we start by proving a couple of results on the tail index estimator γ N,F (k N , k ′ N ). Before that, we introduce some notation: we will write Ψ ∈ RV 2 (a, ∆) where a ∈ R and ∆ is a bounded measurable function having ultimately constant sign and converging to 0 at infinity, such that |∆| is an ultimately monotonic and regularly varying function, if there exists a positive constant c such that

Ψ(y) = cy a exp y 1 ∆(z) z dz with ∆(y) = ∆(y)(1 + o(1)) as y → ∞.
The following second-order condition is required:

(C) We have F ∈ RV 2 (-1/γ F , ∆ F ) and G ∈ RV 2 (-1/γ G , ∆ G ) where γ F , γ G > 0 and |∆ F | ∈ RV 1 (ρ F ), |∆ G | ∈ RV 1 (ρ G ) with ρ F , ρ G ≤ 0.
It can be shown (see Lemma 8) that if (C) holds, then provided ρ F = ρ G and ρ G = -1/γ F , an analogue of condition (C) also holds for F * and G * . Finally, let U F * , U G * be the left-continuous inverses of 1/F * and 1/G * . The following result, in which we write s ∨ t and s ∧ t for the maximum and the minimum of two real numbers s and t and ⌊•⌋ for the floor function, examines the asymptotic properties of γ N,F (k N , k ′ N ).

Theorem 3 Let (k n ), (k ′ n ) be such that k n ∧ k ′ n → ∞ and (k n ∨ k ′ n )/n → 0. Assume that (M) holds. Then we have γ N,F (k N , k ′ N ) P -→ γ F . Suppose moreover that (C) holds, that ρ F = ρ G and ρ G = -1/γ F , that k n ∆ 2 F * (U F * (n/k n )) ∨ k ′ n ∆ 2 G * (U G * (n/k ′ n )) → 0 and sup r,s∈In k r ∧ k ′ r k s ∧ k ′ s -1 → 0 where I n = [np(1 -n -1/4 ), np(1 + n -1/4 )]. ( 4
)
Then if either k n /k ′ n → 0 or k ′ n /k n → 0, we have k ⌊np⌋ ∧ k ′ ⌊np⌋ ( γ N,F (k N , k ′ N ) -γ F ) d -→ N 0, σ 2 F , (5)
where

σ 2 F is equal to γ 2 F (1 + γ F /γ G ) 2 if k n /k ′ n → 0 and γ 4 F /γ 2 G if k ′ n /k n → 0. In the case k n /k ′ n → 1, then we have k ⌊np⌋ ( γ N,F (k N , k ′ N ) -γ F ) = O P (1). (6)
A careful examination of the proof reveals that contrary to Theorems 1 and 2, Theorem 3 also holds when γ F > γ G . Before combining Theorems 2 and 3 to obtain the rate of convergence of the estimator q W N (β n | α n , k N , k ′ N ), we state three remarks about Theorem 3.

Remark 1 Note that without assuming condition (4) the rate of convergence in Theorem 3 would be the random quantity k N ∧ k ′ N , see the proof for further details.

Remark 2 Conditions k n ∆ 2 F * (U F * (n/k n )) → 0 and k ′ n ∆ 2 G * (U G * (n/k ′ n ))
→ 0 are analogues of the condition classically used to prove the asymptotic normality of the Hill estimator. They ensure that the bias of the estimator is negligible with respect to its standard deviation.

Remark 3 Since in the case k n /k ′ n → 1 the correlation between γ N,F * (k N ) and γ N,G (k ′ N ) is very difficult to evaluate, Theorem 3 only provides the rate of convergence of γ N,F (k N , k ′ N ) -γ F to zero.
This case is interesting when γ F < γ G , i.e. if the tail of T is larger than the tail of Y , because then the tail of Y would not be too contaminated as a result of truncation. On the contrary, in the case when the tail of Y is the heaviest (which we cannot consider for the estimation of extreme quantiles), there would be higher confidence in the estimation of γ G than in that of γ F * , which would lead us to take

k n /k ′ n → 0.
The final result of this section gives some asymptotic properties of the estimator

q W N (β n | α n , k N , k ′ N ): Theorem 4 Let α n → 0, β n → 0, k n ∧ k ′ n → ∞ and (k n ∨ k ′ n )/n → 0. As- sume that (M), (3) and (C) hold. Assume that ρ F = ρ G and ρ G = -1/γ F , that β n /α n → 0, nv 2 (q(α n )) → ∞, nv 2 (q(α n ))∆ 2 F (q(α n )) → 0, k n ∆ 2 F * (U F * (n/k n ))∨ k ′ n ∆ 2 G * (U G * (n/k ′ n )) → 0, (k ⌊np⌋ ∧ k ′ ⌊np⌋ )/nv 2 (q(α n )) → 1 and sup r,s∈In k r ∧ k ′ r k s ∧ k ′ s -1 → 0. Then, if γ F < γ G and either k n /k ′ n → 0 or k ′ n /k n → 0, we have v(q(α n )) √ n log(α n /β n ) q W N (β n | α n , k N , k ′ N ) q(β n ) -1 d -→ N 0, σ 2 F . (7)
In the case

k n /k ′ n → 1, or if γ F = γ G and either k n /k ′ n → 0 or k ′ n /k n → 0, we have v(q(α n )) √ n log(α n /β n ) q W N (β n | α n , k N , k ′ N ) q(β n ) -1 = O P (1). (8)

Simulation study

To illustrate the behavior of our estimators, we shall use the following model:

∀y, t > 0, F (y) = (1 + y 1/δ ) -δ/γ F and G(t) = (1 + t 1/δ ) -δ/γ G ,
where δ > 0 and 0 < γ F < γ G . Note that in this situation, ρ F = ρ G = -1/δ. Thus, the larger the value of δ, the smaller the values of |ρ F | and |ρ G | and the slower γ N,F (k N , k ′ N ) converges to γ F , see e.g. de Haan and Ferreira (2006, p.77) for the related situation when there is no truncation. An important bias in the estimation of q(β n ) can then be expected to appear if δ is large. Moreover, the truncation probability is given in our setting by 1p with p = γ G /(γ F + γ G ). In this simulation study, we examine the finite sample behavior of several estimators of the extreme quantile q(β n ) for β n varying in (0, 0.15]. We first consider the two estimators introduced in the present paper, namely q N (β n ) and

q W N (β n |α n ) = q N (α n ) (α n /β n ) γ N,F (αn) , where (α n ) is a sequence in (0, 1) and γ N,F (α n ) is the estimator of the tail-index γ F defined in (1) with k N = k ′ N = ⌊N α n ⌋.
In order to evaluate how important it is to take into account the fact that the data are right-truncated, we also consider the naive estimators

q * N (β n ) = inf{y ≥ y 0 | F * N (y) ≤ β n }, and its extrapolated version q W, * N (β n ) = q * N (α * n ) (α * n /β n ) γ N,F * (α * n ) , where (α * n ) is a sequence in (0, 1) and γ N,F * (α * n ) = γ N,F * (⌊N α * n ⌋).
These last two estimators are computed using only the observations {Y * i , i = 1, . . . , N }, ignoring the fact that these observations are in fact truncated. Note finally that the Weissman-type estimators only depend on the choice of the parameters α n or α * n .

In this simulation study, we generate R = 1000 samples of size n = 200 from the distributions F and G with δ ∈ {1/3, 1}, for γ F ∈ {1/4, 1/2, 1} and p ∈ {0.7, 0.8, 0.9, 0.95}. In each case and for given α n , α * n and β n , we obtain the observations ( q (r)

N (β n ), q W,(r) N (β n |α n ), q (r), * N (β n ), q W,(r), * N (β n |α * n )), r = 1, . . . ,
log 2 q (r), * N (β) q W,(r), * N (β|α) dβ.
The idea behind these criteria is that for quantiles which are not too large, the estimators q (r)

N (resp. q (r), * N ) and q W,(r) N (.|α n ) (resp. q W,(r), * N (.|α * n )) should be close if α n (resp. α * n ) is well chosen. Next, we compute the errors E(q (r) ) := 0.15 0 log 2 q(r) (β) q(β) dβ,
where q(r) is either q (r)

N , q (r), * N , q W,(r) N (.|α (r) opt ) or q W,(r), * N (.|α (r), * opt ).
The error E is a measure of the overall performance of a quantile estimator when estimating extreme quantiles. For θ = {0.1, 0.5, 0.9}, let r(θ) (resp. s(θ), r * (θ) and s * (θ)) be the replication corresponding to the quantile of order θ of the set {E( q (r) N ), r = 1, . . . , R} (resp. {E( q W,(r) N

), r = 1, . . . , R}, {E( q (r), * N

), r = 1, . . . , R} and {E( q W,(r), * N

), r = 1, . . . , R}). In each situation, the errors corresponding to these replications are given respectively in Tables 1, 2, 3 and 4. It appears that the Weissman-type estimators perform better than the others, which is not surprising since these estimators are specifically adapted to the estimation of extreme quantiles. Besides, we can see that the smaller is the probability p, the larger is the bias of the estimators. This was also expected since in our setting, a smaller probability p means a greater number of observations missing in the right tail of Y . Note also that the importance to take into account the fact that the data are right-truncated appears clearly for the Weissman-type estimators when p is small. For large values of p, the naive estimators and the estimators proposed in this paper have similar performances.

Real data example

The dataset we work on here deals with the lifetime of automobile brake pads. It was already considered by Lawless (2002, Example 2.4.2) and obtained in the following way: in order to study the brake pad lifetime, a car manufacturer selected a random sample of cars which were sold over the previous year. In this situation, one way to obtain the brake pad lifetime is to conduct frequent assessments of the brake pad until it is found to be so worn out that it requires replacement, but this procedure is too time-consuming and difficult to implement. Instead, this lifetime was estimated by the manufacturer, based on a measure of the brake pad thickness. For each car, the observed mileage M and the estimated lifetime L were collected, and only the cars with a response variable L larger than M , or equivalently cars whose brake pad thickness was so large that the brake pads did not require immediate replacement, were kept. The random variables L and M can reasonably be assumed to be independent. At the end of this process, only N = 98 estimated lifetimes {L * i , i = 1, . . . , N } and mileages {M * i , i = 1, . . . , N } remained in the sample. Note that since the variable of interest is the brake pad lifetime, this is actually a randomly lefttruncated sample. In order to work within the framework of the present paper, the following transformation is used: for i ∈ {1, . . . , N }, we define

Y * i = (L * i -m N -εt) -1 and T * i = (M * i -m N -ε) -1 ,
where m N = min{M * 1 , . . . , M * N } and ε = 0.05. Thus, (Y * 1 , T * 1 ), . . . , (Y * N , T * N ) can be seen as randomly right-truncated observations from a random sample of independent copies of a random pair (Y, T ), whose size n is unknown. We now need to check that this random sample presents some evidence of heavy tails. To this aim, using our Hill-type estimators, the estimations of the tail indices γ F * and γ G (using the same notation as before) are represented on Figure 7 as functions of k N ∈ {1, . . . , 90}. It appears clearly that the estimated values are positive and, for intermediate values of k N , the estimations of γ F * and γ G seem to be fairly stable. We thus are led to believe that there is indeed evidence of heavy tails for this random sample. We now estimate the extreme quantiles of Y using the Weissman-type estimator q N (β n |α n ), where α n is chosen as in the simulation study. The selected value of α n for this dataset is 0.123, corresponding to k N = 12. The tail indices γ F * , γ F and γ G are respectively estimated to be 0.32, 0.60 and 0.69. In particular, the estimate of γ F is less than that of γ G , which seems to indicate that condition (M) is satisfied. On Figure 7, the estimated quantile q N (β n |α n ) for β n ∈ (0, 0.15) and the one associated to the original data, namely q N (1 -

β n |α n ) = ( q N (β n |α n )) -1 + m N -ε are represented.
In particular, we may see on the right panel of Figure 7 that the brake pad lifetime is estimated to be less than 13 500 kilometers for only 1% of the cars.

Conclusion

In this paper, we introduced and studied an extreme quantile estimator for randomly right-truncated data. This estimator is built upon an empirical high quantile estimator and a tail index estimator which are both adapted for random right-truncation. Our extreme quantile estimator has satisfying performances, be them theoretical or practical.

Future work on our estimator includes obtaining the asymptotic normality of the tail index estimator and the extreme quantile estimator when k n = k ′ n . This is not only a stimulating theoretical problem but also an interesting practical one, since we consider this case in our simulation study and data analysis. Another question worthy of research would be to build and study extreme-value index estimators and extreme quantile estimators when the distribution functions of Y and T belong to an arbitrary max-domain of attraction (de [START_REF] De Haan | Extreme value theory: An introduction[END_REF]. This would certainly be useful in practical applications, for instance when trying to handle random samples whose underlying marginal distributions are believed to be short-tailed.

Proofs of the main results

Proof of Theorem 1. As a preliminary step, note that when N > 0 (which is true with arbitrarily large probability as n → ∞, by Lemma 1), one has

Λ F N (y n ) -Λ F (y n ) F (y n ) = S n,1 + S n,2 + S n,3 - Λ F (y n ) F (y n ) I {Y * N,N ≤yn} where S n,1 = I {Y * N,N >yn} 1 n N i=1 I {Y * i >yn} pF (y n )C * (Y * i ) - Λ F (y n ) F (y n ) , S n,2 = p N - 1 n N i=1 I {Y * i >yn} pF (y n )C * (Y * i ) I {Y * N,N >yn} and S n,3 = I {Y * N,N >yn} N F (y n ) N i=1 I {Y * i >yn} 1 C * N (Y * i ) - 1 C * (Y * i )
.

Lemma 4 entails that I {Y * N,N ≤yn} is zero with arbitrarily large probability as n → ∞ and thus:

v(y n ) √ n Λ F N (y n ) -Λ F (y n ) F (y n ) = v(y n ) √ n 3 j=1
S n,j + o P (1).

Let us first focus on the term S n,1 which we can rewrite as

S n,1 = I {Y * N,N >yn} n n i=1
W n,i where W n,i :=

I {Yi>yn} I {Yi≤Ti} pF (y n )C * (Y i ) - Λ F (y n ) F (y n ) . (9)
It is easy to check that the W n,i are independent, identically distributed and centered random variables. From [START_REF] Gomes | Estimation of the extreme value index for randomly censored data[END_REF] we get

E(W 2 n,1 ) = 1 p 2 F 2 (y n ) E I {Y >yn} I {Y ≤T } (C * (Y )) 2 - Λ F (y n ) F (y n ) 2 . ( 10
)
We now use the fact that since F is nondecreasing and

F (y n ) → 1, we have Λ F (y n ) F (y n ) -1 = 1 F (y n ) ∞ yn F (z) F (z) dF (z) = O(F (y n )) (11)
and by Proposition 1,

C * = p -1 F G so that 1 p 2 E I {Y >yn} I {Y ≤T } (C * (Y )) 2 = ∞ yn dF (z) G(z)F 2 (z) = ∞ yn dF (z) G(z) (1 + o(1)). ( 12
)
• In the case γ F < γ G , noting that dF (z) = -dF (z), we get from ( 10), ( 11), ( 12) and (36) (see the Appendix):

E(W 2 n,1 ) = γ G γ G -γ F 1 F (y n )G(y n ) (1 + o(1)). ( 13
)
Furthermore, because γ F < γ G , one can pick δ > 0 such that (1+δ)/γ G -1/γ F < 0. Hölder's inequality then gives

E[|W n,1 | 2+δ ] ≤ 2 1+δ E I {Yi>yn} I {Yi≤Ti} pF (y n )C * (Y i ) 2+δ + 1 + o(1) ,
where [START_REF] Herbst | An application of randomly truncated data models in reserving IBNR claims[END_REF] was used. Besides, since F is nondecreasing and F (y n ) → 1, we obtain by Lemma 2 with ϕ = 1/G 1+δ and ψ = F :

E I {Yi>yn} I {Yi≤Ti} pC * (Y i ) 2+δ = ∞ yn dF (z) G 1+δ (z)F 2+δ (z) = O F (y n ) G 1+δ (y n ) . (14)
It follows from ( 13), ( 14) and (36) that

n -δ/2 E[|W n,1 | 2+δ ] [Var(W n,1 )] 1+δ/2 = O [nF (y n )G(y n )] -δ = O [v(y n ) √ n] -δ → 0.
Since the W n,i are independent, identically distributed and centered random variables, Lyapunov's central limit theorem (see e.g. Billingsley 1979, p.312) entails √ nS n,1 / Var(W n,1 ) d -→ N (0, 1). Using [START_REF] Hu | Estimation from truncated life time data with supplementary information on covariates and censoring times[END_REF] and the convergence

I {Y * N,N >yn} → 1 leads to v(y n ) √ nS n,1 d -→ N (0, 1) when γ F < γ G . (15)
• When γ F = γ G , we note that using condition (3), the second part of Lemma 2 with ϕ = 1/G and ψ = F entails that the function v is regularly varying with index -1/γ F < 0. This yields 1

F 2 (y n ) ∞ yn dF (z) G(z) → ∞. (16)
Consequently, from [START_REF] De Haan | Extreme value theory: An introduction[END_REF], [START_REF] Herbst | An application of randomly truncated data models in reserving IBNR claims[END_REF] and ( 16), we get

E(W 2 n,1 ) = O(1/v 2 (y n )), which entails v(y n ) √ nS n,1 = O P (1) when γ F = γ G . (17)
Let us now focus on the term S n,2 . From the previous results, it is clear that

1 n N i=1 I {Y * i >yn} pF (y n )C * (Y * i ) I {Y * N,N >yn} P -→ 1.
Since np/N = 1 + O P n -1/2 from Lemma 1, one has, using Lemma 4, that

S n,2 = O P n -1/2 . Using the convergence v(y n ) → 0 it is now obvious that v(y n ) √ nS n,2 = o P (1). ( 18
)
Let us now control S n,3 . Note that (see Woodroofe 1985, pp.172-173):

√ n sup z∈R F * N (z) -F * (z) = O P (1) and √ n sup z∈R G * N (z) -G * (z) = O P (1).
Therefore, it comes that

√ n sup 1≤i≤N C * N (Y * i ) -C * (Y * i ) I {Y * i >yn} = O P (1).
By Lemmas 1 and 5, we thus obtain:

v(y n ) √ nS n,3 = O P v(y n ) nF (y n ) N i=1 I {Y * i >yn} (C * (Y * i )) 2 I {Y * N,N >yn} . (19)
Since

N i=1 I {Y * i >yn} (C * (Y * i )) 2 = n i=1 I {Yi>yn} I {Yi≤Ti} (C * (Y i )) 2 ,
it follows from ( 12) that

E v(y n ) nF (y n ) N i=1 I {Y * i >yn} (C * (Y * i )) 2 = p 2 ∞ yn dF (z) G(z) (1 + o(1)) → 0, (20)
because the integral in the right-hand side converges to 0. Using [START_REF] Lagakos | Nonparametric analysis of truncated survival data, with application to AIDS[END_REF] and ( 20) entails

v(y n ) √ nS n,3 = O P v(y n ) nF (y n ) N i=1 I {Y * i >yn} (C * (Y * i )) 2 I {Y * N,N ≥yn} = o P (1). ( 21
)
Use finally ( 15), ( 17), ( 18) and ( 21) together to get

v(y n ) √ n Λ F N (y n ) -Λ F (y n ) F (y n ) = ξ n if γ F < γ G , O P (1) if γ F = γ G ,
where ξ n is a random variable which is asymptotically standard Gaussian distributed. Using the delta-method concludes the proof of Theorem 1.

Proof of Theorem 2. We start by the case γ F < γ G and we define σ n = q(α n )/[v(q(α n ))

√ n]. It is enough to show that Φ n (z) := P(σ -1 n ( q N (α n )q(α n )) ≤ z) → Φ(z), for every z ∈ R, where Φ is the cdf of a N (0, γ 2 F ) distribution. Let us introduce the sequence ϑ n := γ F v(q(α n ))

√ n/α n . It is easy to check that Φ n (z) = P(W n ≤ z n ), where

W n = ϑ n F N ( q n ) -F ( q n ) and z n = ϑ n (α n -F ( q n )),
with q n = q(α n ) + σ n z. Let us first focus on the nonrandom term z n . Since F is a differentiable function, there exists θ n ∈ (0, 1) such that α n -F ( q n ) = σ n zF ′ (q(α n ) + θ n σ n z). Since σ n /q(α n ) → 0 as n → ∞, we may use the convergence yF ′ (y)/F (y) → 1/γ F to get

F ′ (q(α n ) + θ n σ n z) = γ -1 F α n /q(α n )(1 + o(1)
). Hence the following equality holds:

z n = ϑ n σ n z 1 γ F α n q(α n ) (1 + o(1)) = z(1 + o(1)). ( 22
)
We now consider the random term W n . One has

W n = ϑ n F ( q n ) v( q n ) √ n Z n where Z n = v( q n ) √ n F N ( q n ) F ( q n ) -1 .
Note that from model (M), F ( q n ) = α n (1+o(1)). Moreover, it is a consequence of Lemma 2 that the function v is regularly varying, so that v(

q n ) = v(q(α n ))(1+ o(1)). Consequently ϑ n F ( q n ) = γ F v( q n ) √ n(1 + o(1)
). Apply then Theorem 1 with y n = q n to obtain that Z n d -→ N (0, 1) and thus W n d -→ N (0, γ 2 F ), which concludes the proof in the case γ F < γ G . Now, if γ F = γ G , we start by showing that if (ε n ) is an arbitrary nonrandom positive sequence tending to 0 at infinity such that ε n v(q(α n ))

√ n = ε n q(α n )/σ n → ∞, we have ε n σ -1 n | q N (α n ) -q(α n )| P -→ 0. ( 23
)
Pick then an arbitrary z > 0. We shall show that ϕ n (z) := P(ε n σ -1 n | q N (α n )q(α n )| > z) → 0. With ϑ n as above, it is easy to check that ϕ n (z) = P(W n,+ > z n,+ ) + P(W n,-< z n,-), where

W n,± = ϑ n F N ( q n ) -F ( q n ) and z n,± = ϑ n (α n -F ( q n )),
where we redefine q n := q(α n ) ± ε -1 n σ n z. Let us first focus on the nonrandom term z n,± . Mimicking the arguments leading to [START_REF]A method of allowing for known observational selection in small samples applied to 3CR quasars[END_REF] in the proof of the first part of Theorem 2, we get that

z n,± = ±ϑ n ε -1 n σ n z 1 γ F α n q(α n ) (1 + o(1)) = ±ε -1 n z(1 + o(1)). ( 24
)
We now consider the random term W n,± . One has

W n,± = ϑ n F ( q n ) v( q n ) √ n Z n,± where Z n,± = v( q n ) √ n F N ( q n ) F ( q n ) -1 .
Since F and v are regularly varying, we have

F ( q n ) = α n (1 + o(1)) and v( q n ) = v(q(α n ))(1 + o(1)) which leads to ϑ n F ( q n ) = γ F v( q n ) √ n(1 + o(1)
). On the other hand, the second part of Theorem 1 implies that ε n Z n,± = o P (1), so that using [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF] we obtain for n large enough

ϕ n (z) ≤ P(ε n Z n,+ > z/2) + P(ε n Z n,-< -z/2) → 0
and the proof of ( 23) is complete. Note now that if (ε n ) is an arbitrary nonrandom positive sequence, we have

ε n ≤ ε ′ n := ε n ∨ (v(q(α n )) √ n) -1/2 with ε ′ n v(q(α n )) √ n → ∞.
It can thus easily be seen that in fact ( 23) holds for every positive sequence (ε n ); applying Lemma 6 completes the proof of Theorem 2.

Proof of Theorem 3. Use Lemma 1 to get P(N ∈ I n ) → 1, where I n is defined in equation ( 4). The consistency statement is thus an immediate consequence of Lemmas 3 and 9.

To prove ( 5) and ( 6), write

γ n,F (k N , k ′ N ) -γ F = γ N,F * (k N ) γ N,G (k ′ N ) γ N,G (k ′ N ) -γ N,F * (k N ) - γ F * γ G γ G -γ F * . Since γ N,F * (k N ) P -→ γ F * and γ N,G (k ′ N ) P -→ γ G , it is straightforward to obtain the equality γ n,F (k N , k ′ N ) -γ F = 1 + γ F γ G 2 ( γ N,F * (k N ) -γ F * ) - γ 2 F γ 2 G ( γ N,G (k ′ N ) -γ G ) + o P ( γ N,F * (k N ) -γ F * ) + o P ( γ N,G (k ′ N ) -γ G ).
Applying Lemma 8 proves that F * and G * satisfy the conditions of Lemma 9. Using then Lemma 9 twice concludes the proof.

Proof of Theorem 4. From condition (C), note that we may write

∀y > 0, F (y) = y -1/γ F L F (y) with L F (y) = c F exp y 1 ∆ F (v) v dv
where c F is a positive constant and ∆ F is asymptotically equivalent to ∆ F . Further, since q is the (generalized) inverse function of F , it satisfies the equation

∀α ∈ (0, 1), q(α) = α -γ F L γ F F (q(α)). ( 25
)
Note that since | ∆ F | is asymptotically equivalent to the ultimately monotonic function |∆ F |, one has for n large enough

ε n := 1 log(α n /β n ) log L F (q(β n )) L F (q(α n )) ≤ 2 |∆ F (q(α n ))| log(α n /β n ) log q(β n ) q(α n ) .
By [START_REF] Shorack | Empirical processes with applications to statistics[END_REF] we obtain for n large enough

ε n ≤ 2γ F |∆ F (q(α n ))|(1+ε n ), which entails ε n ≤ 2γ F |∆ F (q(α n ))| 1 -2γ F |∆ F (q(α n ))| = O (|∆ F (q(α n ))|) .
Using once again (25), we get

log q W N (β n | α n , k N , k ′ N ) q(β n ) = log q N (α n ) q(α n ) + ( γ N,F (k N , k ′ N ) -γ F ) log α n β n + O (|∆ F (q(α n ))|) . ( 26
)
To prove [START_REF] Einmahl | Statistics of extremes under random censoring[END_REF], remark that since log(α n /β n ) → ∞, applying Theorems 2 and 3 together with Slutsky's lemma yields, if either

k n /k ′ n → 0 or k ′ n /k n → 0 with γ F < γ G : v(q(α n )) √ n log(α n /β n ) log q W N (β n | α n , k N , k ′ N ) q(β n ) d -→ N (0, σ 2 F ).
Using the delta-method ends the proof of (7). To prove (8 26), Theorems 2 and 3 to get

) if k n /k ′ n → 1 or γ F = γ G , use (
v(q(α n )) √ n log(α n /β n ) log q W N (β n | α n , k N , k ′ N ) q(β n ) = O P (1).
Applying the mean-value theorem to the exponential function ends the proof of Theorem 4.

• If α < β then the function ϕ is integrable with respect to ψ on a neighborhood of infinity and

∞ y ϕ(z)dψ(z) = - β β -α ϕ(y)ψ(y)(1 + o(1)) as y → ∞. • If α = β and ∞ A ϕ(z)dψ(z) < ∞ then y → - ∞ y ϕ(z)dψ(z)
is slowly varying at infinity and

- 1 ϕ(y)ψ(y) ∞ y ϕ(z)dψ(z) → ∞ as y → ∞.
Proof of Lemma 2. We start by the case α < β. Let 2δ = βα > 0, take y ≥ A so large that y -α-δ ϕ(y) ≤ 1 and write for Y > y:

Y y ϕ(z)dψ(z) ϕ(y)ψ(y) = Y y z -α-δ ϕ(z) y -α-δ ϕ(y) - z y -δ z α+δ dψ(z) y α+δ ψ(y) + Y y z α dψ(z) y α ψ(y) . ( 27
) Let now S(y) = sup λ≥1 (λy) -α-δ ϕ(λy) y -α-δ ϕ(y) -λ -δ .
A uniform convergence result for the function y → y -α-δ ϕ(y) (see e.

Besides, Theorem 1.6.5 in [START_REF] Bingham | Regular variation[END_REF] entails for all θ < β

∞ y z θ dψ(z) = - β β -θ y θ ψ(y)(1 + o(1)) as y → ∞. (29)
Using ( 29) with θ = α + δ thus entails that the expression on the left-hand side of (28) has a finite limit as Y → ∞ and we may write

∞ y z -α-δ ϕ(z) y -α-δ ϕ(y) - z y -δ z α+δ dψ(z) = o(y α+δ ψ(y)) as y → ∞. (30)
Combining ( 27), ( 29) with θ = α and (30) concludes the proof of the first statement.

We turn to the case α = β. Pick an arbitrary µ > 1; since ψ is nonincreasing on [A, ∞), we have for y ≥ A:

- ∞ y ϕ(z)dψ(z) ϕ(y)ψ(y) ≥ - µy y ϕ(z)dψ(z) ϕ(y)ψ(y) . ( 31
)
To control the expression on the right-hand side of (31), write

Further, an integration by parts yields

- µy y z α dψ(z) y α ψ(y) = 1 - (µy) α ψ(µy) y α ψ(y) + α µy y z α ψ(z) y α ψ(y) -1 dz z + α log µ.
Since y → y α ψ(y) is slowly varying at infinity, Theorem 1.5.2 in Bingham et al.

z α dψ(z) y α ψ(y) → α log µ as y → ∞. (1987) gives - µy y
Combine then (31), (32), ( 33) and (34) to get for y large enough

- ∞ y ϕ(z)dψ(z) ϕ(y)ψ(y) ≥ - µy y ϕ(z)dψ(z) ϕ(y)ψ(y) = α log µ(1 + o(1)) (35)
which, since µ > 1 is arbitrary, proves that the left-hand side of this inequality tends to infinity as y → ∞. We conclude the proof by noting that from (35), for all λ > 0,

∞ λy ϕ(z)dψ(z) ∞ y ϕ(z)dψ(z) -1 = - λy y ϕ(z)dψ(z) ∞ y ϕ(z)dψ(z) = - ϕ(y)ψ(y) ∞ y ϕ(z)dψ(z) α log(λ)(1 + o(1))
converges to 0 as y → ∞, which is what we wanted to prove.

If ∞ y0 dF (z)/G(z) < ∞, using Lemma 2 with ϕ = 1/G and ψ = F entails F (y)G(y) v(y) → γ G /(γ G -γ F ) if γ F < γ G ∞ if γ F = γ G as y → ∞, (36)
where as in Theorem 1,

v(y) := F (y) ∞ y dF (z) G(z) -1/2 . Lemma 3 Assume that F ∈ RV 1 (-1/γ F ), G ∈ RV 1 (-1/γ G ) where γ F , γ G > 0.
Then as y, t → ∞:

F * (y) F (y)G(y) → 1 p γ G γ F + γ G and G * (t) G(t) → 1 p .
Proof of Lemma 3. Pick y, t > 0 and recall from (37) that

F * (y) = 1 p ∞ y G(z)dF (z) and G * (t) = 1 p ∞ t F (z)dG(z). (37)
Noting that F (z) → 1 as z → ∞ and dF (z) = -dF (z), the result is then a straightforward consequence of Lemma 2 with ϕ = G and ψ = F .

Lemma 4 For all y ≥ y 0 ,

P(Y * N,N ≤ y) = 1 -pF * (y) n . Consequently, if (M) holds, y n → ∞ and nv 2 (y n ) → ∞ then P(Y * N,N ≤ y n ) → 0. Proof of Lemma 4. From Lemma 1, given N = m, the random variables Y * 1 , .
. . , Y * m are independent and identically distributed with cdf F * . Therefore for all y ≥ y 0 , conditioning on N gives

P(Y * N,N ≤ y) = n m=0 n m [pF * (y)] m (1 -p) n-m
which yields the first part of the result. Lemma 3 now entails

n log(1 -pF * (y n )) = -pnF * (y n )(1 + o(1)) = - γ G γ F + γ G nF (y n )G(y n )(1 + o(1)).
Use then (36) to obtain that nF (y n )G(y n ) → ∞; consequently, P(Y * N,N ≤ y n ) → 0, which concludes the proof.

Lemma 5 If y n → ∞ then sup 1≤i≤N Y * i >yn C * (Y * i ) C * N (Y * i ) = O P (1).
Proof of Lemma 5. We closely follow the proof of Lemma 1.2 in Stute [START_REF] Stute | Almost sure representations of the product-limit estimator for truncated data[END_REF]. For 0 < s, t ≤ 1 and λ ≥ 1, set

H n (s, t) = 1 n n i=1 I {Ui≤s, Vi≤t}
where (U 1 , V 1 ), . . . , (U n , V n ) is a sample of independent copies of a random pair (U, V ) which is uniformly distributed on the unit square. We shall prove that for every 0 < a, b ≤ 1 and every µ ∈ (0, 1)

P    sup a≤s≤1 b≤t≤1 - H n (s, t) st ≥ -µ    ≤ exp (-nab h(µ) + 1) (38)
where h(z) = 1 + z(log z -1). To this end, we follow the steps of Lemma 1.1 in Stute [START_REF] Stute | Almost sure representations of the product-limit estimator for truncated data[END_REF], whose proof is itself similar to the proof of inequality (5) p.415 in Shorack and Wellner [START_REF] Shorack | Empirical processes with applications to statistics[END_REF]. For an arbitrary t ∈ [0, 1] let F t be the σ-algebra

F t = σ ({V i ≤ t ′ }, 1 ≤ i ≤ n, t ≤ t ′ ≤ 1, U 1 , . . . , U n ) .
Then (see e.g. Proposition 3.6.2 in Shorack and Wellner [START_REF] Shorack | Empirical processes with applications to statistics[END_REF]) for every s ∈ [a, 1] the process

- H n (s, t) st , F t b≤t≤1
is a reverse martingale. Since the image of a martingale by a convex function is a submartingale, it follows that for every r > 0 the process

sup a≤s≤1 exp -r H n (s, t) st , F t b≤t≤1
is a nonnegative reverse submartingale. Doob's maximal inequality for such processes and the fact that the exponential function is increasing then entail

P    sup a≤s≤1 b≤t≤1 - H n (s, t) st ≥ -µ    ≤ exp(rµ)E sup a≤s≤1 exp -r H n (s, b) sb . (39)
We shall now control the expectation on the right-hand side by another martingale argument: for an arbitrary s ∈ [0, 1] let G s be the σ-algebra

G s = σ ({U i ≤ s ′ }, 1 ≤ i ≤ n, s ≤ s ′ ≤ 1, V 1 , . . . , V n ) .
The same ideas used above show that for every q > 1 and r > 0 the process exp -r q

H n (s, b) sb , G s a≤s≤1
is a nonnegative reverse submartingale. Doob's maximal inequality for q-th moments of a submartingale yields

E sup a≤s≤1 exp -r H n (s, b) sb ≤ q q -1 q E exp -r H n (a, b) ab .
Let q → ∞ and recall (39) to get, for every r > 0,

P    sup a≤s≤1 b≤t≤1 - H n (s, t) st ≥ -µ    ≤ exp(rµ + 1)E exp -r H n (a, b) ab .
Finally, since the expectation on the right-hand side above can be seen as the moment generating function of a binomial random variable at -r/nab,

E exp -r H n (a, b) ab = (1 -ab + ab exp(-r/nab)) n ≤ exp (-nab(1 -exp(-r/nab)))
with the classical inequality (1x/n) n ≤ e -x valid for every x > 0. Putting r = -nab log µ gives (38).

We can now prove our result. Note that

C * N (Y * i ) ≥ 1/N ≥ 1/n so that      sup 1≤i≤N Y * i >yn C * (Y * i ) C * N (Y * i ) ≥ λ      ⊂      sup t≥yn t: C * (t)≥λ/n n N C * (t) C * N (t) ≥ λ      ⊂      sup t≥yn t: C * (t)≥λ/n - N n C * N (t) pC * (t) ≥ - 1 λp      . (40)
Remarking that t → G((-t) -) := lim s↑-t s<-t G(s) is the cdf of -T , we find that

N n C * N (t) = 1 n n i=1 I {Yi≤t≤Ti} = 1 n n i=1 I {Yi≤t, -Ti≤-t} (41) = 1 n n i=1 I {Ui≤F (t), Vi≤G(t -)} (42)
where for each n, (U 1 , V 1 ), . . . , (U n , V n ) is a sample of independent copies of a random pair (U, V ) which is uniformly distributed on the unit square. Finally, recall that pC * (t) = F (t)G(t) to obtain

C * (t) ≥ λ/n ⇒ G(t -) ≥ G(t) ≥ λp/n. (43)
Use now (40) together with (41) and (43) to get

     sup 1≤i≤N Y * i >yn C * (Y * i ) C * N (Y * i ) ≥ λ      ⊂        sup t: F (yn)≤F (t)≤1 t: λp/n≤G(t -)≤1 - H n (F (t), G(t -)) F (t)G(t -) ≥ - 1 λp        .
Noting that log(x)/x ≤ e -1 for all x > 1, inequality (38) now yields, for every λ > 0 such that λp > 1:

P    sup 1≤i≤N Y * i >yn C * (Y * i ) C * N (Y * i ) ≥ λ    ≤ exp((-λp + log(λp) + 1)F (y n ) + 1) = O(exp(-cλ))
where c is a positive constant, of which the result follows.

Lemma 6 Let (X n ) be a sequence of positive real-valued random variables such that for every positive nonrandom sequence (δ n ) converging to 0, the random sequence (δ n X n ) converges to 0 in probability. Then X n = O P (1).

Proof of Lemma 6. Assume that the sequence (X n ) is not stochastically bounded, that is:

∃ε > 0, ∀M > 0, ∃m ∈ N, P(|X m | > M ) ≥ ε.
This entails that there exists a sequence of integers (m k ), which we may choose to be increasing, such that ∀k ∈ N\{0}, P(|X m k | > k) ≥ ε. Let then (δ n ) be the sequence defined as

δ n = 1/k if m k ≤ n < m k+1 .
The sequence (δ n ) converges to 0, and we have ∀k ∈ N \ {0}, P(|δ m k X m k | > 1) ≥ ε which shows that (δ n X n ) does not converge in probability to 0: Lemma 6 is proven.

For an arbitrary Borel measurable function ψ such that ψ(y) = 0 for y large enough and such that z → ψ(z)/z is integrable in a neighborhood of infinity, we define

I(ψ, y) = 1 ψ(y) ∞ y ψ(z) z dz.
Lemma 7 Let ψ ∈ RV 2 (-α, ∆) with α > 0 and |∆| ∈ RV 1 (ρ). Then we have

I(ψ, y) = 1 α + 1 α(α -ρ) ∆(y)(1 + o(1)) as y → ∞.
Proof of Lemma 7. Start by using the change of variables z = yu to obtain

I(ψ, y) - 1 α = ∞ 1 ϕ y (u) du u α+1 with ϕ y (u) = exp yu y ∆(v) v dv -1 (44)
where ∆ is asymptotically equivalent to ∆. Since we are only interested in the behavior of ∆ at infinity, we may assume without loss of generality that ∆ has constant sign. Using the change of variables v = yw, a uniform convergence result for the regularly varying function | ∆| (see Theorem 1.5.2 in [START_REF] Bingham | Regular variation[END_REF]) yields for all u ≥ 1:

ϕ y (u) = ∆(y)D ρ (u)(1 + o(1)) as y → ∞ with D ρ (u) = u 1 v ρ-1 dv. (45)
Furthermore, Theorem 1.5.6 in [START_REF] Bingham | Regular variation[END_REF] entails that there exists

y 1 > 0 such that ∀w ≥ 1, ∀y ≥ y 1 , ∆(yw) ∆(y) ≤ 2w ρ+α/4 .
It is then clear that for y ≥ y 1 :

∀u ≥ 1, 1 | ∆(y)| yu y ∆(v) v dv ≤ 2D ρ+α/4 (u).
Choose now y 1 so large that for every y ≥ y 1 , | ∆(y)| < α/2. The inequality |e t -1| ≤ |t|e |t| , valid for every t ∈ R, entails

∀u ≥ 1, u -α-1 sup y≥y1 ϕ y (u) ∆(y) ≤ 2u -1-α/2 D ρ+α/4 (u). ( 46
)
Because the right-hand side of this inequality defines an integrable function on [1, ∞), we may combine (44), ( 45), ( 46) and apply the dominated convergence theorem to obtain

I(ψ, y) = 1 α + ∆(y) ∞ 1 u -α D ρ (u) du u (1 + o(1)) as y → ∞.
Since ∆ is asymptotically equivalent to ∆, the result is then a straightforward consequence of this equality.

Lemma 8 Assume that (C) holds. Let ρ F * = ρ F ∨ ρ G and define for all y, t > 0:

∆ F * (y) = 1 1 -γ F * ρ F * - γ F ρ F 1 -γ F * ρ F ∆ F (y) + ∆ G (y) 1 -γ F * ρ F * and ∆ G * (t) = F (t) γ F + γ G + ∆ G (t). If ρ F = ρ G and ρ G = -1/γ F then, defining ρ G * = (-1/γ F ) ∨ ρ G , we have |∆ F * | ∈ RV 1 (ρ F * ), |∆ G * | ∈ RV 1 (ρ G * ) and F * ∈ RV 2 (-1/γ F * , ∆ F * ), G * ∈ RV 2 (-1/γ G , ∆ G * ). Proof of Lemma 8. Since |∆ F | ∈ RV 1 (ρ F ), |∆ G | ∈ RV 1 (ρ G ) and F ∈ RV 1 (-1/γ F ), it is clear that |∆ F * | ∈ RV 1 (ρ F * ) and |∆ G * | ∈ RV 1 (ρ G *
) under the conditions of the result. To prove that F * and G * satisfy the desired secondorder condition, we first remark that condition (C) entails that F and G are differentiable. Furthermore, we obtain from (37) that

y d dy log F * (y) = - yF ′ (y)G(y) pF * (y) and t d dt log G * (t) = - tF (t)G ′ (t) pG * (t) .
We start by controlling F * (y). Using (C) leads to

pF * (y) yF ′ (y)G(y) = I(F G, y) -γ F ∆ F (y)I( ∆ F F G, y) 1 -γ F ∆ F (y)
where ∆ F is asymptotically equivalent to ∆ F . The condition ρ F = ρ G entails that the function ∆ F + ∆ G has ultimately constant sign and |∆ F + ∆ G | is asymptotically equivalent to an ultimately monotonic regularly varying function with index ρ F * , and thus F G satisfies the second-order condition (C).

Furthermore, note that since ∆ F has asymptotically constant sign, we have

I( ∆ F F G, y) = I(| ∆ F |F G, y
) when y is large enough. Applying Lemma 7 to the function F G, Theorem 0.6.(a) in [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF] to the function | ∆ F |F G and using a Taylor expansion, we obtain after some algebra:

y d dy log F * (y) = - 1 γ F * + ∆ F * (y)(1 + o(1)) as y → ∞.
An integration completes the proof of the statement on F * . We now turn to the statement on G * . Define Ψ : z → zG ′ (z) and use (C) to obtain that Ψ ∈ RV 1 (-1/γ G ). Using again (C) entails

pG * (t) tF (t)G ′ (t) = F (t) F (t) (I(Ψ, t) -I(GΨ, t)) + I(G, t) -γ G ∆ G (t)I( ∆ G G, t) 1 -γ G ∆ G (t)
where ∆ G is asymptotically equivalent to ∆ G . Moreover, we have

I( ∆ G G, t) = I(| ∆ G |G,
t d dt log G * (t) = - 1 γ G + ∆ G * (t)(1 + o(1)) as t → ∞.
An integration completes the proof.

Lemma 9 Let γ > 0 and Z be a random variable whose survival function Ψ belongs to RV 1 (-1/γ). Assume that:

• N := N (n) is a sequence of integer-valued random variables such that there exists a positive sequence (u n ) of integers tending to infinity with

N/u n P -→ ∞; • γ N (k N ) is a random variable such that the distribution of γ N (k N ) given N = m is that of γ m (k m ) = 1 k m km i=1 log Z m-i+1,m Z m-km,m
where Z 1,m ≤ • • • ≤ Z m,m are the order statistics related to a sample of independent and identically distributed copies Z 1 , . . . , Z m of Z.

Then for every sequence

(k n ) such that k n → ∞ and k n /n → 0, we have γ N (k N ) P -→ γ. Assume further that Ψ ∈ RV 2 (-1/γ, ∆); then if k n → ∞, k n /n → 0 and k n ∆ 2 (U (n/k n )) → 0 where U is the left-continuous inverse of 1/Ψ, the random variable √ k N ( γ N (k N ) -γ) is asymptotically Gaussian cen- tered with variance γ 2 . Proof of Lemma 9. Note that if I n = [u n , ∞) then P(N ∈ I n ) → 1.
To prove the consistency statement, we write for all t > 0, 2006), so that the estimator is indeed consistent. To show the asymptotic normality of γ N (k N ), we let ϕ be the cdf of a centered Gaussian distribution with variance γ 2 and we write for all t ∈ R, P(

P(| γ N (k N ) -γ| > t) = r 1,n + r 2,n with r 1,n = P(| γ N (k N ) -γ| > t, N / ∈ I n ) and r 2,n = m∈In P(| γ N (k N ) -γ| > t|N = m)P(N = m). It is clear that r 1,n ≤ P(N / ∈ I n ) → 0. Besides r 2,n ≤ sup m∈In P(| γ m (k m ) -γ| > t) → 0 since γ n (k n ) P -→ γ,
√ k N ( γ N (k N ) -γ) ≤ t) -ϕ(t) = r 3,n + r 4,n with r 3,n = P( k N ( γ N (k N ) -γ) ≤ t, N / ∈ I n ) -ϕ(t)P(N / ∈ I n ) and r 4,n = m∈In P( k N ( γ N (k N ) -γ) ≤ t|N = m) -ϕ(t) P(N = m). Again, clearly |r 3,n | ≤ 2P(N / ∈ I n ) → 0. Moreover |r 4,n | ≤ sup m∈In P( k m ( γ m (k m ) -γ) ≤ t) -ϕ(t)
so that it is enough to show that γ n (k n ) converges in distribution to a centered Gaussian distribution with variance γ 2 . To this end first note that since Ψ is continuous, we have for all y > 1 that (1/Ψ)(U (y)) = y and thus p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.08 0.10 0.16 0.07 0.08 0.12 0.06 0.07 0.10 0.05 0.06 0.08 γ F = 1/2 0.31 0.38 0.60 0.26 0.31 0.45 0.23 0.27 0.36 0.21 0.25 0.32 γ F = 1 1.22 1.53 2.27 1.05 1.28 1.82 0.91 1.08 1.49 0.85 0.99 1.29 p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.11 0.17 0.36 0.09 0.12 0.25 0.07 0.1 0.16 0.07 0.09 0.13 γ F = 1/2 0.34 0.46 0.85 0.28 0.36 0.55 0.24 0.29 0.43 0.22 0.27 0.37 γ F = 1 1. [START_REF]A method of allowing for known observational selection in small samples applied to 3CR quasars[END_REF] for δ = 1/3 (top) and δ = 1 (bottom). p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.004 0.03 0.22 0.003 0.02 0.10 0.002 0.01 0.06 0.002 0.01 0.04 γ F = 1/2 0.01 0.10 0.50 0.007 0.05 0.27 0.004 0.03 0.16 0.004 0.03 0.12 γ F = 1 0.04 0.39 1.71 0.03 0.25 1.15 0.02 0.13 0.61 0.01 0.09 0.39 p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.05 0.22 2.84 0.04 0.17 1.00 0.03 0.12 0.49 0.03 0.10 0.30 γ F = 1/2 0.04 0. [START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF] opt ) for δ = 1/3 (top) and δ = 1 (bottom). p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.10 0.15 0.21 0.07 0.10 0.14 0.05 0.07 0.10 0.05 0.06 0.08 γ F = 1/2 0.36 0.51 0.69 0.26 0.36 0.51 0.21 0.27 0.37 0.20 0.24 0.32 γ F = 1 for δ = 1/3 (top) and δ = 1 (bottom). p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.04 0.08 0.14 0.01 0.04 0.08 0.003 0.01 0.04 0.002 0.01 0.03 γ F = 1/2 0.12 0.28 0.49 0.04 0.14 0.29 0.007 0.04 0.15 0.004 0.03 0.10 γ F = 1 0.5 1.11 1.90 0.12 0.52 1.14 0.03 0.18 0.61 0.01 0.12 0.45 p 0.7 0.8 0.9 0.95 θ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 γ F = 1/4 0.14 0.26 0.44 0.06 0.14 0.25 0.03 0.08 0.16 0.02 0.06 0. ) for δ = 1/3 (top) and δ = 1 (bottom).

∀y > 1, log y = -log c + 1 γ log U (y) - U (

  R of the estimators. For each replication, α n and α * n are then taken as α

z

  g. Theorem 1.5.2 in Bingham et al., 1987) entails that S(y) converges to 0 as y → ∞. Furthermore, since ψ is nonincreasing, α+δ dψ(z) ≤ -S(y) Y y z α+δ dψ(z).

  result for the function ϕ on compact sets (see e.g. Theorem 1.5.2 in[START_REF] Bingham | Regular variation[END_REF]) and the convergence ψ(µy)/ψ(y) → µ -β as y → ∞ entail → 0 as y → ∞.

y) 1 ∆∀y > 1 , 1 δ

 111 (v) v dvwith c > 0 and ∆ asymptotically equivalent to ∆. Differentiating and carrying out some computations, we getU ′ (y) U (y) = y γ c exp y (v) v dvwhere c > 0 and δ = ∆ • U . Since (Z 1 , . . . , Z n ) has the same distribution as (U (X 1 ), . . . , U (X n )) where (X 1 , . . . , X n ) is a sample of independent copies of standard Pareto random variables, we may then write γ n (k n ) d = S n,1 + S n,2 with representation and Paul Lévy's central limit theorem entailk n (S n,1γ) (0, γ 2 ). Moreover, since (n/k n )X n-kn,n P -→ 1 (see Corollary 2.2.2 p.41 in de Haan and Ferreira 2006) and |δ| is asymptotically equivalent to the regularly varying and ultimately nonincreasing function |δ| := γ 2 |∆ • U |, we get with arbitrarily large probabilityk n |S n,2 | ≤ k n |δ(X n-kn,n )|S n,1 (1 + o P (1)) = O P ( k n |δ(n/k n )|) P -→ 0as n → ∞. Using Slutsky's lemma completes the proof of Lemma 9.

Figure 1 :Figure 2 :

 12 Figure 1: Hill estimators of γ G (full line) and γ F * (dashed line) as functions of k N .

  t) when t is large enough, because ∆ G has asymptotically constant sign. Apply then Theorem 0.6.(a) in[START_REF] Resnick | Extreme Values, Regular Variation, and Point Processes[END_REF] to the functions Ψ, F Ψ and | ∆ G |G, Lemma 7 to the function G and use a Taylor expansion to get after some computations:

  see Theorem 3.2.2 p.70 in de Haan and Ferreira (

Table 1 :

 1 1.55 2.46 1.05 1.29 1.9 0.92 1.09 1.46 0.86 1.01 1.33 Errors associated with the estimator q

(r(θ)) N

Table 2 :

 2 2.43 0.03 0.14 0.85 0.02 0.09 0.42 0.02 0.07 0.27 γ F = 1 0.05 0.46 2.65 0.03 0.25 1.42 0.02 0.15 0.66 0.02 0.11 0.53 Errors associated with the estimator q

	W,(s(θ)) N	(.|α	(s(θ))

Table 3 :

 3 1.43 1.97 2.64 1.00 1.42 1.95 0.85 1.06 1.48 0.80 0.98 1.31 F = 1/4 0.19 0.33 0.53 0.10 0.18 0.31 0.07 0.1 0.18 0.06 0.08 0.13 γ F = 1/2 0.46 0.71 1.03 0.29 0.45 0.69 0.23 0.30 0.46 0.21 0.26 0.38 γ F = 1 1.43 2.13 2.91 1.06 1.48 2.09 0.86 1.10 1.53 0.81 0.99 1.33 Errors associated with the estimator q

	p		0.7			0.8			0.9			0.95
	θ	0.1	0.5	0.9	0.1	0.5	0.9	0.1	0.5	0.9	0.1	0.5	0.9
	γ

(r * (θ)), * N

Table 4 :

 4 14 γ F = 1/2 0.22 0.42 0.75 0.07 0.20 0.41 0.02 0.08 0.21 0.02 0.06 0.16 γ F = 1 0.47 1.16 2.04 0.16 0.53 1.15 0.03 0.21 0.60 0.02 0.12 0.43 Errors associated with the estimator q

W,(s * (θ)), * N (.|α (s * (θ)), * opt

Appendix -Preliminary results

The first result gives an equivalent of the random variable N and the conditional distribution of (Y * 1 , T * 1 ), . . . , (Y * N , T * N ) given N .

Lemma 1 We have N/np = 1 + O P (n -1/2 ). Furthermore, the conditional distribution of (Y * 1 , T * 1 ), . . . , (Y * N , T * N ) given N = m > 0 is equal to the distribution of m independent copies of a random vector (Y * , T * ) with cdf H * .

Proof of Lemma 1. The first part of the result is a straightforward consequence of Chebyshev's inequality. Let now (A i ) i≥1 be arbitrary Borel subsets of

which concludes the proof.

Lemma 2 Let ϕ ∈ RV 1 (α) and ψ ∈ RV 1 (-β) with α ∈ R and β > 0. Assume that ψ is right-continuous and nonincreasing on some interval [A, ∞), A ≥ 0.


Rachid Ouaret
email: rachid.ouaret@cstb.fr

Anda Ionescu

Olivier Ramalho

Yves Candau

Viorel Petrehus

Lucille Labat

Topic B6: Prediction and measurement MODELLING THE TIME FLUCTUATION OF INDOOR AIR FORMALDEHYDE CONCENTRATION: VARIABILITY STRUCTURE IDENTIFICATION AND FORECASTING USING NONLINEAR MODELS

Keywords: formaldehyde, clustering, forecasting, SOTA, SETAR

This study explores the possibility to forecast formaldehyde (HCHO) concentration from past observations in an office. A monitoring campaign of HCHO was performed during 96 days, with a short time-step (1 minute). The monitored formaldehyde time series exhibited, in particular, abrupt changes and structural breaks in variance, that cannot be modelled implicitly using simple models. To overcome this problem, hybrid model for forecasting HCHO time series (characterized by nonlinear and nonstationary behaviour) was used. The periodicities of the data were modelled by Fast Fourier Transform (FFT). A Self-Exciting Threshold AutoRegressive (SETAR) model was used to model the FFT component. SETAR models are typically designed to accommodate the nonlinear features and can explain two regimes in a time series. The residuals of the FFT component subtracted from raw data was modelled using a second SETAR model. The output of the two models were summed and compared to raw test data.

INTRODUCTION

The indoor concentration of formaldehyde (HCHO) fluctuates over time and remains difficult to model because of a number of different emission sources (including reactivity) and its sensitivity to temperature, humidity and ventilation changes. The concentration dynamic and its relation to sources and meteorological variables are complex as they involve the dependency of sources and sinks as well as other parameters such as temperature, humidity, air flow but also to the concentration of oxidants in the air (i.e. radicals, peroxides, ozone, oxygen). While the outdoor air research has experience in long monitoring studies with advanced statistical modelling at various time scales, the field of indoor air quality lacks monitoring studies of pollutant concentrations with a short time-step. In particular, the temporal variability of indoor formaldehyde fluctuations remains unknown. Identification of predominant scales of temporal variability of HCHO concentration is a prerequisite to understand the heterogeneity of indoor source emissions (type and intensity). Information on the magnitude of variation and underlying mechanism is important both to improve our understanding of HCHO dynamics and key characteristics of source emission.

The indoor concentration of a pollutant, e.g. formaldehyde, can be modelled by two different manners: physical models or statistical models. The physical models are based on general mass balance integral-differential equations. These models are commonly used in experimental chambers, where all parameters are controlled and do not vary in time. The use of such models in real environments require the knowledge of all contributing sources and their emission parameters, as well as the sink effects, and how all these parameters fluctuates over time. On the opposite, statistical models require less a priori knowledge, i.e. information can be retrieved from a concentration time-series by inverse modelling and provide insights on the contributing sources and possible forecast.

The objective of this study is to explore the statistical methods to forecast indoor formaldehyde concentration using past observations. The variability structure of the past observations is the key parameter to achieve this goal. It is determined through the analysis of abrupt point changes and clustering of similar trends in the data, which oriented us towards the use of a hybrid SETAR model (Self-Exciting Threshold Autoregressive) to predict future concentrations.

METHODS

Measurement of HCHO

The indoor formaldehyde concentration was monitored every minute in an open-plan office occupied by 6-8 persons over 96 days, using an AL4021 analyzer (Aerolaser GmbH). The detection of formaldehyde is based on the liquid phase reaction of formaldehyde with acetyl acetone and ammonia that produces 3,5-diacetyl-1,4-dihydrolutidine which is measured by fluorescence at 510 nm. The device was regularly calibrated with liquid solutions of formaldehyde and the internal permeation source in the range 0 to 55 ppb. The detection limit was 0.36 ppb. The sampling point was placed in the exhaust ventilation duct.

The office is located in a suburban area 30 km east of Paris. The total area of the office is 132 m 2 with a volume of 397 m 3 . A constant air flow rate is maintained through a mechanical exhaust ventilation system with air inlets provided in the windows. The total exhaust ventilation rate is 252 m 3 /h. The indoor materials were carpet on the floor, painted walls and ceiling tiles. The furniture comprises typical L-shaped desks in melaminated particleboards and aluminium closets. A laser multifunction copier was in use in the office space.

Statistical methodology

The first part of this analysis is to explore HCHO variability structure. The main daily profiles were obtained by classification. Several classification techniques were employed: hierarchical classification (CAH) based on the Ward algorithm [START_REF] Saporta | analyse des données et statistique[END_REF], k-means [START_REF] Saporta | analyse des données et statistique[END_REF] and SOTA (Self Organizing Tree Algorithm) [START_REF] Dopazo | Phylogenetic reconstruction using an unsupervised growing neural network that adopts the topology of a phylogenetic tree[END_REF][START_REF] Herrero | A hierarchical unsupervised growing neural network for clustering gene expression patterns[END_REF]. The first methods (CAH and k-means) are classical techniques, widely used for clustering. SOTA algorithm has been developed in biology and gene expression, but it can be used in other fields as well.

With the variability of HCHO concentration, the problem is to find the abrupt changes in the time series, which probably indicates a source of variability. We are interested in abrupt point changes detection from a statistical point of view. That corresponds to a problem of estimating the point at which the statistical properties of a sequence of observations change.

For this purpose, we apply the Pruned Exact Linear Time (PELT) method [START_REF] Jackson | An algorithm for optimal partitioning of data on an interval[END_REF]. Unfortunately, the choice of the appropriate penalty of the algorithm optimization depends on many factors including the size of the changes and the length of segments, which is still an open question.

The presence of abrupt changes revealed by PELT method and present in some daily profiles obtained by classification justified the use of a particular autoregressive model to forecast HCHO from past observations, which is called Self-Exciting Threshold Autoregressive model or SETAR (Tong et al.1980[START_REF] Tong | Threshold models in time series analysis -30 years on (with discussions by[END_REF]. A 2-stage procedure was employed. First, a Fast Fourier Transform (FFT) was applied to raw data. This first component was modelled by a first SETAR. The residuals between raw data and the FFT component were calculated and a second FFT was applied to the residuals. The results were further modelled by a second SETAR. The results of the two forecast procedures were summed. Before performing the forecast procedure, the initial raw database was split in two parts: the first one used for 'training' (model fitting) and a second one, independent, used for test. A comparison between forecast and raw data from the test ensemble was made using classical performance criteria. A general framework of the statistical procedure developed in this study is given in Figure 1. In order to define typical daily profiles, a matrix of 96 rows was designed, one row for each day. The mean daily profile over the 96 days is plotted in Figure 3a and the other daily profiles (median, minimum, maximum, quartiles and percentile 95) in figure 3b. The mean concentration of HCHO is around 10 ppb during night hours and decreased during occupancy hours at 8.8 ppb. This may be explained by higher air exchange rates during the day due to higher frequency of windows and doors openings, and thus maybe more ultraviolet radiation entering the office which increased the photolysis rate of formaldehyde.

A) B)

Figure 3: (A) Daily mean HCHO profiles for 96 days (B) Daily profile of HCHO maximum, minimum, mean, median, quartiles (Q25, Q75) and percentile (P95) for 96 days.

Classification

The 96 daily time series were considered as individuals and different clustering algorithms were applied in order to reveal the main patterns. Several classifications methods (HAC, kmeans, SOTA) were used and the results were similar. SOTA classification results are presented in Figure 4.

The interpretation of the different clusters should be done using some exogenous variables.

During the same campaign, other pollutants (ozone, nitrogen oxides) and climatic parameters were monitored. A complementary analysis that integrates all these parameters will be an important part of our future work.

There are some 'flat' profiles, characterized by a small variation amplitude (cluster 1), but also profiles with abrupt changes (especially clusters 3 and 4).

Abrupt change detection

We detected some abrupt changes in the HCHO using the PELT algorithm. An example is given in Figure 5.

Forecast

In order to take into account this type of variations (abrupt changes), we selected a SETAR model to forecast HCHO concentrations. The database comprises 21 days with no missing data (or very short gaps that could be interpolated), which corresponds to 30,240 observations. A training sample of 27,600 observations was created and also a test sample with the remaining 2,640 observations. The model was fitted on the training dataset and evaluated on the test set.

Figure 6a presents the training set and figure 6b the test set. In the upper plot, raw HCHO concentrations are plotted in black and the concentration filtered using the FFT, in red. The residuals obtained as the difference between raw data and FFT component are given in the lower plot, in black. A second FFT component was calculated for theses residuals and plotted in red. A first SETAR model was fitted to the first FFT component (calculated from the raw data) and a second SETAR model was fitted to the second FFT component (calculated from the residuals). Two forecast sets were calculated using the two SETAR models and then summed; the result is plotted in blue on Figure 6b. The forecast (blue line) can be compared to raw test data (in black). In addition, a FFT of the raw test data was plotted in green.

The forecast model fits quite well the reconstructed FFT component (global trend + residuals) up to 900 min (15 hours) and was satisfactory during the 44 hours of the test set (less than 2 ppb difference). The model is thus useful to predict the general variation of HCHO but unsufficient to predict fast variations change. In order to achieve that, other informative data must be input in the model. That is part of the second step of our work. The criteria to evaluate the forecast performance on the filtered FFT components are given in With a forecasting horizon of 44 hours, the performance of the model is satisfactory to predict the general variation of formaldehyde concentration (FFT test set). When compared to the raw data set, more discrepancies appear as the model do not take into account high frequency variations of formaldehyde. However, the performance of the model is not bad with 66% of the total variance explained.

CONCLUSIONS

The analysis of HCHO time series revealed the presence of some abrupt changes that could be modelled using a SETAR model. The forecast performance of the general trend (derived from FFT filtering) was good for a horizon of 44 hours (almost 2 days ahead). The performance of the model should be increased using exogenous variables if the statistical analysis reveals their importance in the ongoing work.

The classification permitted to obtain some diurnal profiles. These profiles should be further analysed using other pollutant concentrations and climatic parameters, as well as ventilation conditions.

Figure 1 .

 1 Figure 1. General framework for statistical analysis and forecast of HCHO concentrations.

  The concentration of HCHO monitored at 1-minute time step during 96 days(138,240 observations) in the open-plan office are presented in Figure2. The indoor concentration of HCHO remains low with no value observed above 25 ppb. The measurements started on 2013-04-27 and ended on 2013-07-31. There are some gaps, especially during summer holidays. Long periods of missing data lead us to reduce the database used for forecast purposes, in order to have a continuous time series. The final database used for prediction was 21 days long, from 2013-05-24 to 2013-06-14.

Figure 2 :

 2 Figure 2: Time series of 96 days of HCHO sampled every minute

Figure 4 :

 4 Figure 4: SOTA clustering of 96 daily time series of HCHO (0 to 1440 minutes).

Figure 5 :

 5 Figure 5: Example of detecting abrupt change using the PELT algorithm.

Figure 6 .

 6 Figure 6. (A): Training set and residuals (raw data and FFT). (B) Test set of 44 hours (raw data and forecast)

Table 1 :

 1 Forecast performance criteria.

			Forecasting Horizon = 44 h
		1st FFT-SETAR	2nd FFT-	FFT test set and	raw test set and FFT-
	Performance	(HCHO)	SETAR	FFT-SETAR	SETAR forecast
				forecast
	ME (ppb)	-0.96	0.53	0.43	0.48
	RMSE	1.15	1.65	1.58	2.57
	MAE (ppb)	1.06	1.29	1.18	1.80
	MPE (%)	-10.49	1.58	6.40	4.37
	MAPE (%)	11.21	176.73	12.43	16.75
	R 2	-	-	0.91	0.66


Rachid Ouaret
email: rachid.ouaret@cstb.fr

Anda Ionescu

Olivier Ramalho

Viorel Petrehus

Yves Candau

Forecasting indoor pollutant concentrations using Fast Fourier Transform and Regime Switching Models

Keywords: Forecasting, Threshold Autoregression (TAR), Fast Fourier Transform (FFT), high-frequency, formaldehyde

This paper considers the possibility that the high-frequency formaldehyde concentration follows a regime switching process. It explores the calibration of models built on Threshold Autoregression (TAR) combined with Fast Fourier Transform (FFT) for short-term forecasting of indoor formaldehyde fluctuation. The methodology uses signal decomposition with filtering of the original raw time series into different FFT components removing some high frequencies. Results from FFT-TAR model could predict HCHO concentration patterns with acceptable accuracy and suggest that there is benefit by taking choice of "cutoff" frequency. The mean absolute percentage error (MAPE) of the best scheme for 10-h (600); one day (1440) and 40h ahead out of sample forecasts, were obtained as follows: 3.54%, 14.36% and 15.64%.

Introduction

Due to high levels of indoor pollutants and usually 80% of our time spent indoors, significant health problems (eg. respiratory and cardiovascular) can be caused by exposure to indoor air. Therefore, estimating future values of indoor pollutant concentrations is of particular interest in the field of indoor air quality (IAQ). One of the major contaminant in indoor air is formaldehyde (HCHO). It originates from indoor materials such as particleboard, plywood, and combustion (e.g. tobacco smoke) and photochemistry. The dynamics of indoor formaldehyde concentration in real world environments and its relation to sources and climatic variables are complex and little studied. While the outdoor air research has much experience in long monitoring studies with advanced statistical modeling, the field of indoor air quality lacks monitoring studies of pollutant concentrations with a short time-step. The majority of papers published on air quality forecasting uses low-frequency data, i.e. weekly, monthly, quarterly, annual, etc. However, there are very few studies on shorter time interval environmental time series. Granger [START_REF] Granger | Extracting information from mega-panels and high-frequency data[END_REF] argued that high-frequency data pose a new set of forecasting problems, demanding new approaches.

In this paper, we propose a new approach in order to predict the indoor pollutant concentration of formaldehyde. The emission mechanisms of HCHO from materials are complex and depend on several parameters that vary over time like temperature and humidity. The temporal variability of these parameters is most often unknown and their physical models are often difficult to build for indoor air pollutants. In this work we develop a method to forecast shortterm formaldehyde indoor concentrations using nonlinear time series models and Fourier transform. The succession of phases with high and low variability and the presence of several mean shifts in HCHO concentration over time can be highlighted by regime switching models. In such a case, the dynamic behavior of HCHO depends on the regime that occurs at different time scales. It seems therefore natural to expect that Regime Switching Models may improve prediction of HCHO indoor concentrations. The forecasting model can also be used to provide valuable feedback for the physical processes. To our knowledge, this is the first study exploring HCHO concentration forecasting by nonlinear stochastic time series. The determinist part of HCHO concentration was modeled by multiple Fast Fourier Transform (FFT) which were incorporated within two types of Regimes Switching Models: Smooth Threshold Auto-Regressive (STAR) and Self Exciting Threshold Autoregressive (SETAR) models. In order to compare their performance, the root mean square error (RMSE) and the mean absolute percentage error (MAPE) were calculated for each forecast.

Different sub-samples of the given historical data were used to compare the different forecasting models. The predicted data were calculated up to 2 days ahead at most. The lack of references about forecasting in the field of indoor air quality is the main reason for us to use a multiple stage and hybrid approach. From that perspective, hybridization of many and varied models using complementary and common approaches and strategies is leading to a kind of robustness. Only two-regime switching models are considered in this study. There are mainly two reasons for this: on the one hand, it is difficult to check the stability of multiple regime parameters requiring several statistical hypothesis and, on the other hand, there is no fully elaborate forecasting methodology for these models.

All the computations involved in the present task were carried out on R (http://www.r-project.org/). The packages Nonlinear Time Series Analysis (RTISEAN) [START_REF] Hegger | Practical implementation of nonlinear time series methods: The TISEAN package[END_REF] and Nonlinear time series models with regime switching (ts-Dyn) -with personal very small changes - [START_REF] Di Narzo | tsDyn: Time series analysis based on dynamical systems theory[END_REF][START_REF] Matthieu | tsDyn: Threshold cointegration: overview and implementation in R[END_REF] were freely used; all the other complementary functions and simulations were coded by us.

Data, statistical proprieties and methodology

The indoor formaldehyde concentration was monitored every minute in an openplan office occupied by 6-8 persons over 21 days, i.e. 30,000 observations. The trajectory, the histogram, the periodogram and the autocorrelation function are provided in the following figures (Fig. 1 (A)). The observed mean is 7.85 ppb, the standard deviation equals to 3.64 and the kurtosis is not very high (-0.4). The fluctuation of the HCHO time series shows a relative stability during the first days and some small oscillations during the two days that followed. After that, the general behavior of the time series changed abruptly and became very variable toward the end of the period. The figure 1(C) presents the histogram and density estimates of the distribution of HCHO concentration which reveal three modes and two gaps around 7 ppb and 10 ppb. The data show an asymmetrical distribution with high density concentrations at the left side and a decreasing density towards the right side. The histogram shape reflects the presence of at least two different processes being "mixed" in the HCHO time dynamics. This indicates that the mean of the process shifted over the measuring period. We observe that the data are not second order stationary; there exists a trend with abrupt changes. This explains the form of the autocorrelation function (ACF) that does not quickly decrease towards zero and the spectral density explodes in many frequencies close to zero (fig 1(B and D)). Thus, according to [START_REF] Beran | Statistics for Long-Memory Processes[END_REF] and [START_REF] Guégan | Les chaos en finance: approche statistique[END_REF] definitions, the HCHO concentration seems to exhibit a long memory behavior in the covariance sense. Fig. 1 (B) shows the periodogram close to zero of HCHO time series; the frequencies corresponding to dominant periods were chosen to construct the first FFT component which represents the general trend, the following FFT components were identified by the same method. Once significant peaks have been isolated from the raw data and the successive residuals; the next step is to rebuild the HCHO in the time domain. The HCHO time series is a mix of several unknown source emissions that makes the extraction of dominant frequency quite difficult because of the random noise. As shown in the Fig. 1(B) (more clearly in Fig ( 3a)), the noise manifests itself across the entire spectrum, including within high frequencies where the time series should have little energy.

Fig. 2 illustrates the steps describing the research methodology used in this study. The process begins first with the identification of the non linear class model which is a TAR model. The next step involved identifying the adequate frequencies cutoff of FFTs components which are obtained successively by subtraction (5 time). The FFTs components were modeled using SETAR and LSTAR models by varying the parametersm and τ ; that are described later in this paper (see [START_REF] Beran | Statistics for Long-Memory Processes[END_REF]). In reality, using whether SETAR or LSTAR could change very slightly the accuracy of prediction, thus only the SETAR accuracy is presented. In model selection-prediction, both forecasting accuracy (in test set) and classical indicators -Akaike information and Schwarz information criterion -(in learning test) were used to evaluate the most appropriate models1 .

Models

Threshold models

The underlying idea behind the construction of a TAR model for IAQ is to represent the behavior of HCHO concentration by two separate regimes. Generally, the transitions between regimes in Threshold models (TAR) models are ruled by mechanisms performing with a transition variable, a threshold value and a function of transition. In this study, we restrict our attention to models for which the transition mechanism between regimes is assumed to be observed and, in each regime the dynamics behavior is described by an autoregressive model with an identical delay. Initially proposed extensively by Tong ([7,[START_REF] Tong | Threshold models in non-linear time series analysis[END_REF][START_REF] Tong | Non-linear Time Series: A Dynamical System Approach[END_REF]), the TAR model is defined as follows: for a k -dimensional time series {X t } t∈Z and for all t, J t is the transition variable relative to threshold value c; the pair (X t , J t ) is an TAR with two regimes T AR(2; p, p)

X t = Φ 1 0 + Φ 1 1 X t-1 + • • • + Φ 1 p X t-p (1 -I Jt>c ) + Φ 2 0 + Φ 2 1 X t-1 + • • • + Φ 2 p X t-p (I Jt>c ) + ε t (1)
where c ∈ R and ε t is an i.i.d white noise process conditional upon the

Ω t-1 = {X t-1 , X t-2 , . . . , , X t-p } with E (ε t | Ω t-1 ) = 0 and E ε 2 t | Ω t-1 = σ 2 .
The transition function is given by:

I Jt>c = 1 is J t > c 0 otherwise . (2)
A special case arises when the transition variable J t is taken as the lagged value of the time series itself, that is, it should be any of {X t-1 , X t-2 , . . . , , X t-p }. The resulting model for any lagged value of X t taken as J t in equation ( 1) is Self-Exciting TAR (SETAR). Because of the discrete nature of the threshold variable, the transition between regimes is carried out abruptly. In fact, there are some periods where the time series does not transit abruptly between regimes; therefore it is possible to replace the transition function in (2) by a continuous function G (J t ; γ, c) bounded between 0 and 1; the resulting model is called Smooth Transition AR (STAR) model and is given by:

X t = Φ 1 0 + Φ 1 1 X t-1 + • • • + Φ 1 p X t-p (1 -G (J t ; γ, c)) + Φ 2 0 + Φ 2 1 X t-1 + • • • + Φ 2 p X t-p (G (J t ; γ, c)) + ε t , (3)
where the parameter γ determines the smoothness of the change in the value of the transition function. In this paper, we choose as in many other works ( [START_REF] Teräsvirta | Specification, Estimation, and Evaluation of Smooth Transition Autoregressive Models[END_REF]), the logistic function:

G (J t ; γ, c) = 1 1 + exp {-γ (J t -c)} (4)
The function G (J t ; γ, c) in the STAR models switches instantaneously between the two regimes as soon as the quantity (J t -c) changes its sign and this, for γ sufficiently large. Consequently the function G (J t ; γ, c) approximates the indicator function I Jt in the model TAR (in equation 1) (or SETAR if J t = X t-d (in equation 3). When γ → 0, (G (J t ; γ, c) 1 /2) then the STAR model reduces to an AR linear model.

In our work several simulations were carried out to determine γ; the choice of γ is done according to the accuracy of prediction (we have avoided the two extremes mentioned above), and have set it at one (γ = 1).

In this work, we have chosen to reformulate the models (1) and ( 3) by introducing the delay time series of embedding dimension m and delay time τ , and the threshold delay δ in both equations system.

X t = Φ 1 0 + Φ 1 1 X t-τ + • • • + Φ 1 (m-1) X t-(m-1)τ (1 -G (J t ; γ, X t-δτ )) + Φ 2 0 + Φ 2 1 X t-τ + • • • + Φ 2 (m-1) X t-(m-1)τ (G (J t ; γ, X t-δτ )) + ε t (5)
All simulations were carried out by varying m, τ and δ for de-noising FFT components and the choice was done once the models gave us a good forecast. As we can see, doing all possible combinations is not an easy task due to the large data set, because it requires a lot of computational resources. As the frequency of the time series increases the required computational time for simulation become high, demanding a good combination in the models parameters that will make the most efficient best setup.

Removing noise and time series reconstruction using FFT

The first motivation for working in the frequency domain is that removing noise from a temporal signal can be accomplished, at least approximately, by removing the signal's "high frequency components". The second reason is that many operations are easier, both computationally and conceptually, in the frequency domain. In general way, the harmonic analysis decompose time series over oscillatory waveforms that reveal many properties and provide a path to sparse representations.

FFT is an efficient method for computing the discrete Fourier transform (DFT). The algorithm for computing the DFT was primarily attributed to Cooley and Tukey [START_REF] Cooley | An Algorithm for the Machine Calculation of Complex Fourier Series[END_REF]. Computationally, the DFT is of the order of O T 2 while the FFT is of the order of O (T log 2 T ). If {x t } t=0;...,(T -1) is the time series of HCHO, its DFT can be written as:

H k = T -1 t=0
x t e -2πj kt T ; with k = 0, . . . , T -1

As mentioned in the section (2) the random noise pollutes the signal. Thus a very simple approach to noise reduction is to determine the "cutoff" frequencies; this was performed by observing the different resolution after reconstruction i.e based in visual insight. This is done by the inverse FFT (IFFT)

x t = 1 T T -1 t=0 H k e 2πj kt T (7)
Once high frequencies removed from the spectrum, the associated signals can be reconstructed in the time domain using IFFT. Each reconstructed FFT component into the time domain corresponds to the real part of the HCHO Fourier term.

In Figure 3 (a-f) we show the reconstruction based on FFT of the HCHO time series; which represents different patterns and characteristics inside the raw time series. The first FFT component presents the general trend that captures the most important information in terms of representativeness. The rest of oscillatory components shows the different inter-day variability. The second FFT component exhibits roughly the same amplitude over 15 days, after that the shape of the oscillations changed near the end of the meausrement duration. The amplitude modulations of the 3rd, 4th and 5th FFT oscillations are similar, but with small differences in variance. Several simulations were carried out in order to determine the various plausible scenarios. In other word, we begin by whether fixing m and τ and varying the cutoff frequency, or vice versa. The cutoff varying frequencies of each FFT and the parameters of SETAR model for modeling the two different sets are shown in Tab 1. The combined models are compared among each other with corresponding data set for multiple step-ahead forecast.

Forecasting

The aim of the nonlinear time series considered in our study is to employ the models for forecasting future values of the HCHO variability. Furthermore, outof sample forecasting was considered as the way to evaluate the capability of regime switching model to capture some of the non-linearity within the data. The forecasting domain has been influenced, for a long time by linear models such as ARMA for which the forecasting function is not very difficult. For example, for AR(2) model, the one step-ahead forecast at time t is xt+1 = φ 1 x t + φ 2 x t-1 and it not difficult to show that in general it holds that: However, considering our case, the time series exhibits nonlinear patterns: asymmetric distribution, abrupt changes, breaks ...etc. Thus forecasting from nonlinear models is a more involved task than forecasting from linear models. That is, the optimal h-step-ahead forecast of x t+h at time t for a general nonlinear model

xt+h|t = E [x t+h | Ω t ] ≈ φ 1 xt+h-1|t + φ 2 xt+h|t ( 8
)
x t = f (x t-1 ; θ) + ε t , (9)
is given by:

xt+h|t = E [x t+h | Ω t ] = E f x t+(h-1) ; θ | Ω t . (10)
The main difference between the forecast in [START_REF] Teräsvirta | Specification, Estimation, and Evaluation of Smooth Transition Autoregressive Models[END_REF] and the general recursive forecast for an ARIMA (for example) is that the linear expectation operator E cannot be interchanged with f : the expected value of nonlinear function is not equal to the function evaluated at the value of its arguments [START_REF] Franses | Non-Linear Time Series Models in Empirical Finance[END_REF]. Overall, the relationship between forecasts at different horizons for the nonlinear model does not exist as shown for linear one in the [START_REF] Tong | Non-linear Time Series: A Dynamical System Approach[END_REF]. Forecasting 2-step ahead for an nonlinear model is given by

xt+2|t = E [f (f (x t ; θ) + ε t+1 ; θ) | Ω t ] = E f xt+1|t + ε t+1 ; θ | Ω t . (11)
Unlike the naive approach which considers that ε t+1 = 0 in [START_REF] Dijk | Smooth Transition Autoregressive Models-A Survey of Recent Developments[END_REF] or the possibility of interchanging between E and f , Monte Carlo simulation and bootstrapping methods are used to compute nonlinear forecasts and to approximate the conditional expectation [START_REF] Dijk | Smooth Transition Autoregressive Models-A Survey of Recent Developments[END_REF].

The 2-step-ahead Monte Carlo forecast is given by

x(mc) t+2|t = 1 k k i=1 f xt+1|t + ε i ; θ , ( 12
)
where k is some larger number and the ε t are simulated from prior information about the distribution of ε t+1 . While in the bootstrap forecast consider the residuals from the estimated model, {ε t } t=1,...,n in, the forecast function is given by:

x(boot) t+2|t = 1 k k i=1 f xt+1|t + εi ; θ . (13)

Results and discussion

The forecasting experiment is designed as follows. We use in-sample HCHO data to estimate the parameters of the models of interest for each simulation. Then we make multi out-of-sample forecasts up to 60 minutes, 10 hours, 1 day, 36 hours and 40 hours ahead and evaluate them. We carry out and evaluate forecasts in Fig. 6: The indoor HCHO forecating of 40h prediction time horizon from scenarion 5; 48 h forecast from scenario 7 and 12.5 h from scenario 8 (raw data and forecast) (here d is τ ).

terms of RMSE and MAPE, we repeat the procedures for each parameter m, τ in the equation [START_REF] Franses | Non-Linear Time Series Models in Empirical Finance[END_REF].

The forecast performance of applied FFT-LSTAR and FFT-SETAR models were evaluated against 60 minutes, 10 hours, 1 day, 36 hours and 40 hours out of sample values; the figures 4 and 5 shows the forecasting performances of the different scenarios models in terms of RMSE and MAPE. For a short-term forecasting 60-600 minutes, the scenarios 4, 5 and 6 outperform the other cases with RMSE around 0.3 ppb and 0.85 ppb respectively. The MAPE values for 60 and 600 out of sample forecasts were below 10%. This shows that the effect of cutoff frequencies in the FFTs have not a significant impact in short-term forecasts. However, if we add more high frequencies in the different FFTs; we reduce consistently the long-term accuracy of forecasting. Specifically, with the same m and τ in scenario 5, the periods included in each FFT are as follows: > 1 hour (FFT1); > 30 minutes (FFT2); > 15 minutes FFT3; > 7.5 minutes (FF4); > 5 minutes (FFT5). The RMSE values for one day out of sample (1440 min) in this configuration is 3,61 ppb, while it is 2.19 ppb for the scenario 5 (MAPE=28 %; 16% respectively). In general, for the model based on scenario 1, 2, 5, and 6; there is no significant difference with RMSE when the prediction time horizon change between 24 h, 36 h and 40 h.

The models studied in scenario 7 highlights a good performance until 1-day ahead (RMSE =1.1 ppb); and the model obtained from scenario 8 highlights very good accuracy until 12.5 h. The box plot of RMSE illustrates the difference between the prediction in the case of set 1 and set 2 as a function of the forecast horizon; as the prediction time interval increases, the actual HCHO fluctuations deviated from its assumed concentrations, which lead to poor forecasting accuracy. As shown in the (Fig [START_REF] Guégan | Les chaos en finance: approche statistique[END_REF]), the HCHO concentrations given by the prediction model of scenario 8 is closely aligned with the raw data. In fact this scheme is the best forecast we got for a 10 hours, with RMSE equal to 0.44 ppb and MAPE equal to 6.14 %. This performance could attributed to the basic underlying characteristic of the good compromise of the included period in FFTs components and the parameters of the SETAR model.

Conclusion

This study deals with forecasting indoor pollutant concentrations and in particular HCHO in an open place office, the focus goes beyond the mere estimation and take account of the forecast ability of Regime Switching Models coupled with the fast Fourier transform components. The results of forecast suggest that there is benefit by taking choice of "cutoff" frequency. We have shown that de-noise the raw data permits to identify some characteristics such as beaks. Furthermore, the non-linear models developed here capture the asymmetric behavior of the HCHO fluctuations.

The methodology developed here appears new in order to forecast indoor air quality, the practical implication of the kind of study is that to help decision-makers to carry out actions that attempt reduce indoor pollution. To do this, it depends on the expected values of indoor pollutants concentrations.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Descriptive statistics of HCHO: trajectory, raw periodogram, density and ACF.

  0.00006 ω included > 4.5h (f) Sum of the five FFTs.

Fig. 3 :

 3 Fig. 3: Multi-FFT reconstruction of original HCHO time series (scenario 6).

Fig. 4 :

 4 Fig. 4: RMSE of FFT-SETAR (2) for each scenario.

Fig. 5 :

 5 Fig. 5: MAPE of FFT-SETAR (2) for each scenario

Table 1 :

 1 Example of scenarios simulations.

	Defining sets scenarios Simulation conditions	FFT1 FFT2 FFT3 FFT4 FFT5 Periods (ωi) included in the FFTs
		scenario1	S1; m = 3, τ = 1	2.08 days 25 h 16.6 h 12.5 h 10 h
		scenario2	S2; m = 3, τ = 1	25 h 16.6 h 12.5 h 10 h 7.14 h
	Set1	scenario3 scenario4	S3; m = 3, τ = 1 S4; m = 3, τ = 1	16.6 h 12.5 h 10 h 8.33 h 6.25 h 12.5 h 10 h 8.33 h 7.14 h 5 h
		scenario5	S4s1*	12.5 h 10 h 8.33 h 7.14 h 5 h
		scenario6	S4s2*	12.5 h 10 h 7.14 h 5.5 h 4.5 h
	Set2	scenario7 scenario8	S6s3** S5; m = 3, τ = 1	12.5 h 10 h 7.14 h 5.5 h 4.5 h 10 h 8.33 h 7.14 h 6.25 h 4.5 h

* FFT1 was modeled with m = 2, τ = 1 and FFT2 to FFT5 with m = 3, τ = 1. ** FFT2 was modeled with m = 3, τ = 1 and others were modeled with m = 2, τ = 1. For Set 1: learning set=27600; test set=2399; and for Set 2: learning set=15000; test set=2880.

In diagnostic checking tests, we only consider (for now) the accuracy is important, even some models may exhibit the presence of unit root in both regimes.