Zhiqiang Wang

Mathieu Ritou
email: mathieu.ritou@univ-nantes.fr

Catherine Da Cunha

Benoît Furet

Contextual classification for smart machining based on unsupervised machine learning by Gaussian mixture model

Keywords: Industry 4.0, machining, unsupervised machine learning, GMM, BIC

Intelligent machine-tools generate a large amount of digital data. Data mining can support decision making for operational management. The first step in a data mining approach is the selection of relevant data. Raw data must, therefore, be classified into different groups of contexts. This paper proposes an original contextual classification of data for smart machining based on unsupervised machine learning by Gaussian mixture model. The optimal number of classes is determined by the silhouette method based on the Bayesian information criterion. This method is validated on real data from four different machine-tools in the aerospace industry. Manual data mining and k-fold cross validation confirm that the proposed method provides good contextual classification results. Then, several key performance indicators are calculated using this contextual classification. They show the relevancy of the approach.

Introduction

High-Speed Machining (HSM) has significantly increased spindle speed compared to conventional machining. Production systems are very flexible, especially in the aeronautics industry, where hundreds of tools machine several thousands of parts with different references on the same site. It is not easy for the Manufacturing Department to identify problematic programs or tools. The main problems are chatter (instability during machining), tool breakage and over-vibrations. Therefore, a data mining system is necessary to protect the machine-tool and work-pieces (especially for aeronautical parts which present high value-added). [START_REF] Teti | Advanced monitoring of machining operations[END_REF] make an overview of the existing techniques of machining monitoring and list all measurable and studied defects. Maleki et al. (2017) proposed a FMEA-based approach for the choice of the sensors. Quintana and Ciurana (2011) focus on the state of the art on chatter and existing methods to detect and to avoid it. [START_REF] Zhou | Review of tool condition monitoring methods in milling processes[END_REF] make a summary of the methods used to monitor tool wear in milling processes, including sensors, feature extraction and monitoring models. Godreau et al. (2019) propose a vibration criterion for the detection of chatter, a source of non-quality.

In the general context of Industry 4.0, large volumes of manufacturing data are available on instrumented machine-tool. They are interesting to exploit not only to improve machine-tool performance but also to support the decision making of operational management. Lenz, Wuest, and Westkämper (2018) propose a holistic approach to machining data analysis in order to combine tasks and to group analysis objectives among different departments within the company. Morgan and ODonnell (2018) present the design and development of a cyber-physical process monitoring system for machinetool. Liu et al. (2018) propose a systematic development method for machine-tool 4.0. In order to extract data from databases more efficiently and accurately, contextual classification is necessary to understand the status of the machine-tool and to realize a more relevant calculation of Key Performance Indicators (KPI). Many researches propose the use of machine learning algorithms directly on the raw data (machine-tool signals). However, few researches focus on data pre-processing to identify the context. In this article, machine learning algorithms for contextual classification are developed. They will permit a more relevant calculation of Key Performance Indicators.

Machine Learning is a field of Artificial Intelligence that relies on statistical approaches. It can be applied to machine-tool signals analysis. Kim et al. (2018) list and summarize machining contributions using machine learning algorithms. There are two categories: supervised and unsupervised machine learning.

-Supervised machine learning learns to classify from already labeled output samples. It aims at making correct predictions on the data which is not present in the learning set. For example, El-Mounayri and Deng (2010) designed and implemented a neural networks (ANN) for force prediction to optimize the cutting parameters in optimizing 2 1 / 2 -axis milling; C. Zhang and H. Zhang (2016) show that the least squares support vector machine (SVM)-based tool wear model is suitable to predict tool wear at certain range of cutting conditions in milling operation; [START_REF] Rodríguez | A decision-making tool based on decision trees for roughness prediction in face milling[END_REF] propose a decision-making tool based on decision trees for roughness prediction in face milling. -Unsupervised machine learning learns to classify unlabeled data. In industry, it is challenging to collect labeled output data. As a result, most industrial applications use unsupervised machine learning methods. In a previous study, [START_REF] Wang | Comparison of K-means and GMM methods for contextual clustering in HSM[END_REF] compared the different methods of unsupervised machine-learning, and Gaussian Mixture Model (GMM) was the most effective for contextual classification. Another challenge of smart-machining comes from the management of time series; the state of the machine can change whenever and has to be identified every tenth of a second.

In this article, unsupervised machine learning by Gaussian Mixture Model (GMM) is used. GMM is a statistical tool used to parametrically estimate the distribution of random variables by modeling them as a mixture of Gaussian distributions (Reynolds 2015). The distributions' parameters are estimated from the learning data using the Maximization Expectation algorithm. The main advantage of GMM is that it can correctly classify data with unbalanced classes, while its disadvantage is that the results of GMM directly depend on the choice of parameters. An additional algorithm can therefore be useful to determine the necessary number of classes automatically [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. The Bayesian Information Criterion (BIC) is particularly recommended.

This paper proposes an original contextual classification of data for smart machining based on unsupervised machine learning by Gaussian mixture model. The method is applied to four machining databases collected in the aeronautics industry. The number of classes is determined by the BIC method. The classifications are then evaluated using manual data mining and k-fold validation. This paper is organized as follows. Section 2 introduces the framework of the study as well as the contextual classification method. Section 3 proposes the method based on Gaussian Mixture Model (GMM). The next section depicts the silhouette method, which determines the optimal number of classes. Section 5 applies the proposed method to vibration signals and power signals to know whether the tool is cutting materials or not. Section 6 verifies the obtained classification using manual data mining and k-fold cross-validation. Section 7 presents KPIs based on the obtained contextual classification. These indicators will aid decision-making within the company.

2. Data analysis process

Framework of the study

The work presented in this article is a part of the ANR SmartEmma project, which aims at developing intelligent and connected machine-tools. The machine learning methods will be developed on the one hand to perform contextual classifications, and on the other, to aggregate the data. The latter concept can be found in [START_REF] Ritou | Knowledge-based multi-level aggregation for decision aid in the machining industry[END_REF], which proposes a method to address the issue of data volume using multi-level aggregation. 1 GB of data is collected per day and per machine-tool. Contextual classification permits to classify raw data into ordered and meaningful forms. Using this contextual information, noise is reduced and relevant KPIs can be calculated. Those KPIs support decisions on the HSM process and the global management of the company.

Contextual classification method

As previously mentioned, contextual classification enables a more relevant evaluation of KPIs. In the context of our research work, it is interesting to know:

-if the spindle is stopped (S) or rotating at constant speed (CS) or varying speed (VS); -if the machine-tool is stopped or moving at constant or varying speed; -if the tool is cutting materials or not (N), with a constant (CC) or varying tool engagement (CV).

From these three contextual information, 17 potential states for the machine-tool can be defined (cf. Table 1). The objective is to determine the state of the machine-tool at each instant. Indeed, it will enable accurate data selection and analyses dedicated to the considered state. stant speed (N > 0 and ∆N ≤ T ∆N ); or varying speed (N > 0 and ∆N > T ∆N ).

In order to know if the tool is cutting materials or not, a classification according to the criteria Arms (root mean square-RMS of vibration acceleration value, in m/s 2 ) and the spindle power P (kW) is proposed, when the machine is moving (feedrate V f > 0). In machining, it is known that if the tool is cutting, its vibration (Arms) is greater than if the tool is not cutting. It is also known that when the tool is not cutting, the power (P ) will exceed the spindle power corresponding to idle rotation (P idle ). However, there is a subclass called spindle startup where the Arms is also greater due to the lubrication spindle startup. This subclass should be identified and assigned to the cluster tool not cutting. The proposed business rules are:

-if (P > P idle or Arms > T Arms ) and V f > 0, except the subclass spindle startup, the tool is cutting material; -if P ≤ P idle and Arms ≤ T Arms , including the subclass spindle startup, the tool is not cutting material.

The following procedure aims to check if the GMM classifications are correct or not:

-Modeling : The contextual classification thresholds are determined, using GMM on different descriptors (∆N , ∆V f , Arms, P and ∆P ). Each machine-tool has different characteristics, so it is necessary to perform an initial training for each machine-tool. One-day data of a machining database is used for this purpose. -Classifying: The classifications are carried out using GMM threshold and the business rules. -Validation: The quality of classifications is assessed by confusion matrix, using a vector labeled with actual classes of machine-tool states.

Gaussian Mixture Model

The Gaussian Mixture Model is a statistical model used to parametrically estimate the distribution of random variables by modeling them as a mixture of Gaussian distributions. It has been used in several areas such as audio recognition, measurement error modeling (Huang et al. 2005). The principle of GMM is to find an approximation of the probability density f from a sample of n observations of the random vector X = {x 1 , x 2 , x 3 ...x n } (where each observation x i is assumed independent) using a linear combination of several Gaussian components.

To evaluate the quality of the model parameter θ, the probability density associated with the observed data is calculated (Maugis, Celeux, and Martin-Magniette 2009):

f (X, θ) = k j=1 a j N j (x i |µ j , σ j ) (2)
where, k is the number of Gaussian; a j is the weight of the j th Gaussian. It verifies the probability conditions such as k j=1 a j = 1, and 0 ≤ a j ≤ 1; N j (x i |µ j , σ j ) expresses the probability that a measure x i is associated with the j th Gaussian. In this context, the parameters of the model are collectively represented as θ = {a j , µ j , σ j }.

In the literature, there are several methods for estimating these GMM parameters. In the present paper, the Maximum Likelihood method is used to determine the parameters θ.

The main idea of Maximum Likelihood Estimation (MLE) is to find a set of parameter estimates, so that the likelihood of the sample used is maximized. Since the n observations are independent, the likelihood function L(X|θ) can be written as (Bishop 2006):

L(X|θ) = f (x 1 , θ) × f (x 2 , θ) × • • • × f (x n , θ) = n i=1 f (x i , θ) (3)
Then substitute Eq. ( 2) into Eq. ( 3) leads to:

L(X|θ) = n i=1 k j=1 a j N j (x i |µ j , σ j ) (4)
It is easier to maximize log-likelihood than the likelihood function itself. The loglikelihood function is written as:

ln L(X|θ) = n i=1 ln ( k j=1 a j N j (x i |µ j , σ j )) (5)
The problem of estimation by the maximum likelihood method is therefore to find the roots of the following equation:

∂ ln L(X|θ) ∂θ = 0 (6)
In general, maximizing the likelihood function does not have an analytical solution. It is then necessary to use iterative methods. The most common is to use the EM algorithm (Expectation-Maximization) (Feng, Lai, and Kar 2016). However, the EM algorithm needs initial parameters; the most important of them is the number of classes k.

In general, the larger the number of parameters to be estimated into the model, the better the model approaches to the data. However, the complexity of the computation will also increase when there are too many parameters to be estimated. Therefore, a compromise must be found; this is the goal of the the Bayesian Information Criterion (BIC). The BIC is defined as follow [START_REF] Schwarz | Estimating the dimension of a model[END_REF]:

BIC(k) = 2 × M ax(ln(L(X|θ))) -N m ln(n) (7)
where, M ax(ln(L(X|θ)))is the maximized log-likelihood function, N m is the number of parameters to be estimated in this model (3 per Gaussian: weight, mean and standard deviation), n is the number of observations, and k is the number of classes.In the BIC, the term (N m ) penalizes the complexity of the model. Therefore, BIC is maximized for most parsimonious setting, i.e., with small N m , and larger likelihood. Consequently, BIC is a good indicator to determine the number of classes required.

Validation of the Gaussian model by silhouette method

In this section, the previously proposed machine learning procedure by GMM is applied to database from four different machine-tools in two factories producing aeronautical structural parts:

-The machine-tool #1 is a Forest Aeromill for aluminum machining, in company #A; -The machine-tool #2 is a Forest Aerostar for aluminum machining, in company #A; -The machine-tool #3 is a Mazak Mazatrol for aluminum machining, in company #B; -The machine-tool #4 is a Mazak i800T for titanium machining, in company #B.

An Emmatools device was installed on each machine-tool, and the machining data were collected for months during industrial productions. The classifications should determine:

-If the spindle is rotating (at constant or varying speed) or not, according to the classification on the descriptor ∆N ; -If the machine-tool is moving (at constant or varying speed) or not, according to the classification on the descriptor ∆V f ; -If the tool is removing materials (with constant or varying tool engagement) or not, according to the classification descriptor of Arms, P and ∆P .

For each machine-tool, the determination of the number of context classes is crucial and required for the Gaussian model that follows (in order to identify the classification thresholds for the business rules). The proposed method is:

-Apply GMM for each descriptor (∆N , ∆V f , Arms, P and ∆P ) for different numbers of classes (1, 2,..., k). For each number of classes, a BIC value can be determined according to Eq. ( 7), in order to evaluate the model. -Compute the BIC in relation to the number of classes for each descriptor. In the present case, each BIC value converges to a maximum. -Chose the minimal number of context classes required In this paper, a class number based on 90% of the BIC maximum is chosen. In this way, the class number is reduced and enables better interpretability of the results.

Silhouette method on the variation of spindle speed

The Gaussian model analyzed by BIC on the variation of spindle speed ∆N is done on the four machine-tools. Fig. 2 shows the results: the BIC value reaches 90% of the BIC maximum with 2 classes. Thus, ∆N was classified into 2 classes, 2 clusters that correspond to the constant and varying spindle rotation. The experimental data confirm what was assumed in Table 1.

Silhouette method on the variation of axis speed

The Gaussian model analyzed by BIC on the variation on ∆V f are performed on four machine-tools as shown in Fig. 3. The figure shows that the BIC value reaches 90% of the BIC maximum with 2 classes. Thus, ∆V f was classified into 2 classes, 2 clusters

Silhouette method on vibrations and power

The descriptors Arms (vibration) and P (spindle power) are used to determine if the tool is removing materials or not (cf. 1). To determine the number of classes required, Gaussian model analyses by BIC on Arms as shown in Fig. 4 and P Fig. 5. 4 different machine-tools are used in order to validate the results.

The figure shows that the BIC value reaches 90% of the BIC maximum with 3 classes. Thus, k=3 and the vibrations Arms and the power P were classified into 3 classes. The 3 classes are interpreted as: not cut, low material removal and high material removal.

Contextual classification of the material removal

Once the class number determined by BIC, the classification thresholds are identified from one-day data by GMM. Then, the business rules are applied, leading to the classification. This section presents the classification. GMM is applied to the collected data in the context of spindle rotation at constant speed and axis motion (V f > 0) see Table 1. Two datasets illustrate the approach: machine-tools #2 and #3 for aluminum machining.

Modeling and classification for machine-tool #3

The distribution of Arms in this class (N constant and V f > 0) is modeled by 3 Gaussians (see section 4).

-a Gaussian represents the class tool not cutting (Y1 green), -another the class tool cutting (Y3 blue), -the 3rd Gaussian represents 'severe machining vibrations' (Y2 red). The distribution of Y4 (in cyan) is the sum of these 3 Gaussians (Y4=Y1+Y2+Y3) as shown in Fig. 6(a). The threshold T Arms is defined as the intersection between the Gaussians Y1 and Y3, in order to minimize classification errors (false positives and negatives). The threshold is identified: T Arms =5.93 m/s 2 is obtained.

Applying the same method, P is modeled by 3 Gaussians:

-one Gaussian represents the class tool not cutting (Y1 green), -another represents the class tool cutting with low material removal (Y3 blue), -the 3rd Gaussian (Y2 red) models the class tool cutting with high material removal.

The distribution of Y4 (cyan) is the sum of these 3 Gaussians (Y4=Y1+Y2+Y3), as shown in Fig. 6(b). The threshold is chosen at the intersection between the Gaussian tool is not cutting and the other Gaussian tool cutting with high material removal, thus obtaining P idle =1.68 kW.

A combination of the two sensors measurements is performed. Thus, for data where spindle is rotating at constant speed and axes are moving (V f > 0) see Table 1, if Arms > T Arms or P > P idle , the data is assigned to the class tool cutting.

Fig. 7 presents data from machine-tool #3 as well as the result of the classification. The tool not cutting class is represented by black "+" signs on the P power curve, the tool cutting class is presented by red and blue dots. The tool cutting class can be classified into two sub-classes: tool cutting with constant engagement (in red points) and tool cutting with varying engagement (in blue points) by classification on ∆P (see Table 1). When the spindle has reached its target speed (N =11250 rpm in Fig. 7), the tool does not cut immediately. There is a subclass called spindle startup. This class is represented by green "+". It is defined by the following rule: from the reach of the target speed N until the power exceeds the no-load power P idle (which represents the beginning of machining). Identify this state permits to eliminate parasitic vibrations related to lubrication spindle startup, which would otherwise lead to classification errors. The subclass spindle startup is also assigned to the cluster tool not cutting.

Cyan line represents V f , blue line N , pink line Arms. The tool engagement is correctly detected.

Modeling and classification for machine-tool #2

The Arms distribution is modeled by 3 Gaussians (see section 4):

-a Gaussian represents the class tool not cutting (Y1 green), -another the class tool cutting (Y3 blue), -the 3rd Gaussian represents 'severe machining vibrations' (Y2 red).

The distribution of Y4 (in cyan) is the sum of these 3 Gaussians (Y4=Y1+Y2+Y3) as shown in Fig. 8(a). The threshold T Arms is defined as the intersection between the Gaussians Y1 and Y3, in order to minimize classification errors (false positives and negatives). This enables to identify T Arms =5.43 m/s 2 for the machine-tool #2.

Applying the same method, P is modeled by 3 Gaussians:

-one Gaussian represents the class tool not cutting (Y1 green), -another represents the class tool cutting with low material removal (Y3 blue), -the 3rd Gaussian (Y2 red) models the class tool cutting with high material removal.

The distribution of Y4 (cyan) is the sum of these 3 Gaussians (Y4=Y1+Y2+Y3), as shown in Fig. 8(b). The threshold is chosen at the intersection between the Gaussian tool is not cutting and the other Gaussian tool cutting with high material removal. This leads to P idle =2.57 kW. In summary, three classes of the tool engagement are defined.

-tool cutting when Arms > T Arms or P > P idle ; it consists of two subclasses : tool cutting with constant engagement and tool cutting with varying engagement; -the tool not cutting when Arms ≤ T Arms and P ≤ P idle by including the subclass spindle startup.

Fig. 9 plots the N (rpm), V f (m/min), Arms (m/s 2 ) and P (kW) as well as the classification results classes (dots and + sign) for machine-tool #2. The material removal is well classified.

Method validation

The estimation of the quality of the classification by machine learning isessential to choose an optimal descriptor or combining descriptors for a given data set [START_REF] Wolpert | Stacked generalization[END_REF]. Two criteria constitute the quality: accuracy and sensitivity.

Accuracy

To estimate the accuracy of the classification by GMM, manual data mining was used. In our approach, the classifications based on GMM model conduct to the predicted classes (tool cutting or tool not cutting). The actual classes correspond to 100% true classes that were constituted by manual data mining. The main steps are as follows: Table 2 reports the confusion matrix for 4 hours on the machine-tool #3 and Table 3 the confusion matrix for the machine-tool #2 . Accuracy of the Classification (ACC) is the ratio of correct predictions and is computed as follows: ACC = T P +T N T P +T N +F P +F N

The accuracy is obtained for the machine-tool #3: ACC #3 (4h)=99.94%, as well as for the machine-tool #2 is: ACC #2 (4h)=99.21%. The classifications obtained using GMM method are very accurate.

Sensitivity of the method

In order to test whether the method is sensitive to the training data set , the crossvalidation method can be used. Kohavi et al. (1995) present different cross-validation methods like 'Holdout', 'Leave-one-out', 'k-fold' and 'Stratification'. In this paper, the cross-validation method 'k-fold' is used. The original sample is divided into k samples, then one of the k samples is selected as the training set and the (k -1) other samples will constitute the validation set. The performance score is calculated. The operation is repeated k times so that each sub-sample is used once as a training set. The average of the performance scores is finally calculated to estimate the accuracy of the classifications. Here 10 days were chosen (k=10) as the original sample. The main steps are as follow. For each day:

-Apply GMM to the training set to its thresholds of Arms and P , -Apply these new thresholds to the validation set (nine other days) to obtain the predicted classes, -Compare these predicted classes with the actual classes and calculate the confusion matrix. The accuracy for each days is obtained.

The thresholds of Arms and P for each day in the training sets as well as corresponding accuracy when applying these thresholds to the validation set for the machine-tool # 3 and the machine-tool # 2 are respectively reported in Fig. 10( a andb). The horizontal dotted lines are the thresholds of Arms and P as well as the corresponding accuracy obtained in the validation set in these two machine-tools. With 10 sets, the mean accuracy of the classification for the machine-tool #3: ACC(k-fold) #3 (4h)=97.60%, and for the machine-tool #2: ACC(k-fold) #2 (4h)=96.28%. This shows the robustness of the proposed method. The mean accuracy is lower than the one using the same set for training and validation. It is because the thresholds for materials removal or not (T Arms and P idle ) may slightly change as time goes by, due to the variation of the machine condition (e.g. wear of the spindle).

Nevertheless, the mean accuracy obtained is high enough that the thresholds learned by GMM during one-day can be used for any other day of production on the same machine-tool.

Application of the contextual classifications for KPIs evaluation on production data

According to the previous validation steps, the contextual classification made by the proposed method is accurate and robust. Thus, Key Performance Indicators can be calculated based on this contextual classification. For example, it is interesting for the industry to know:

-Smart Overall Equipment Effectiveness (OEE) of the machine-tool. (1) Complete business rules for contextual classification, (2) A methodology, based on GMM and BIC, to determine the different states of the machine-tool, Those proposition were validated on the raw data used in this paper was collected and stored by the device Emmatools which is installed on four machine-tools in two factories producing aeronautical structural parts.

-A classification of the tool cutting or not, by combining GMM thresholds of vibration measurements Arms and power measurements P , was validated on 2 different machine-tools. -k-fold validation enables to conclude that the thresholds learned by GMM during one-day can be applied to any other days for the same machine-tool. It was also found that the thresholds learned by GMM during one-day data may slightly change due to the variation of the machine condition. -this contextual classification permit to obatined precisel KPIs that facilitate decision making.

The learned values of thresholds given in this paper are specific to a given machinetool, but the proposed machine learning methodology is replicable to any machine-tool.

  Fig. 1(a) shows an overview of the road-map of this project. A device, called Emmatools, collects the measured data during machining and stores it in a database (De Castelbajac et al. 2014). It is installed on machine-tools in factories of aircraft manufacturers. Fig. 1(b) presents the schematic diagram of the Emmatools. It collects data every tenth of a second. The data comes from two information sources: the Computer Numerical Control (CNC) and the added sensors. By field bus, the CNC provides machining context data such as the tool name, the current program, machine, and spindle motions. The sensors integrated into the machine-tool collect the instantaneous power of the spindle, temperatures, as well as the vibrations from 4 accelerometers integrated into the spindle (radially at each bearing).

Figure 1 .

 1 Figure 1. Data mining process for machining (a) and schematic diagram of the data collect device (b)

Figure 2 .

 2 Figure 2. BIC on ∆N issued from the 4 machine-tools

Figure 3 .

 3 Figure 3. BIC on ∆V f issued from the 4 machine-tools

Figure 4 .

 4 Figure 4. BIC on Arms (when N constant et Vf>0) issued from the 4 machine-tools

Figure 6 .

 6 Figure 6. Distribution of Arms and P (when N constant and Vf>0) modeled by GMM for machine-tool #3

Figure 7 .

 7 Figure 7. Classifications of material removal for machine-tool #3

Figure 8 .

 8 Figure 8. Distribution of Arms and P (when N constant and Vf>0) modeled by GMM for machine-tool #2

Figure 9 .

 9 Figure 9. Classifications of material removal for machine-tool #2

Figure 10 .

 10 Figure 10. The new thresholds of Arms and P obtained in the training sets as well as these thresholds obtained in validation sets for (a) machine-tool #3 and (b) machine-tool #2

Fig. 11

 11 Fig.11stresses the real machine-tool use rate during the complete spindle lifetime (215 days) for machine-tool # 2. The contextual classification enables to determine at each instant if the machine-tool is stopped (in red), removing material (cyan and yellow) or not (in pink). Then, during machining, operations dedicated to high material removal (in cyan) or low material removal (i.e. finishing operations, in yellow) can be aggregated separately (contrary to classical OEE computations). Note that the factory where data were collected works in a three-shift system 5 days a week, and that stopped times of Fig.11include week-end.

Figure 11 .

 11 Figure 11. Machine-tool OEE during the complete spindle lifetime for the machine-tool # 2

Figure 12 .

 12 Figure 12. Pie charts representing the chatter occurrences per workpiece program (a) and per tool (b) during the complete spindle lifetime for machine-tool # 2

Table 1 .

 1 Classification of the different states of the machine-tool

				Feed rate of the machine-tool Vf
					Vf>0
		Vf=0	Arms≤T Arms	∆Vf >T ∆Vf Arms>T Arms or P>P idle	Arms≤T Arms	∆Vf ≤ T ∆Vf Arms>T Arms or P>P idle
			& P≤P idle	∆P≤T ∆P	∆P>T ∆P	& P≤P idle	∆P≤T ∆P	∆P>T ∆P
	N=0	Spindle: S
	Spindle speed N

Table 2 .

 2 Confusion matrix for 4 hours data from the machine-tool #3

	Confusion Matrix	Actual Classes Tool Cutting Tool Not Cutting
	Predicted Classes	Tool Cutting Tool Not Cutting	TP=51 476 FN=67	FP=18 TN=92 470

Table 3 .

 3 Confusion matrix for 4 hours data from the machine-tool #2

	Confusion Matrix	Actual Classes Tool Cutting Tool Not Cutting
	Predicted Classes	Tool Cutting Tool Not Cutting	TP=60 626 FN=938	FP=292 TN=93 707

Acknowledgment

The authors acknowledge the financial support of the French National Research Agency on the ANR SmartEmma project (ANR-16-CE10-0005). They would also like to thank the contributions made by the industrial partners.

Bhinge, Raunak, Nishant Biswas, David Dornfeld, Jinkyoo Park, Kincho H Law, Moneer Helu, and Sudarsan Rachuri. 2014. "An intelligent machine monitoring system for energy prediction using a Gaussian Process regression." In 2014 IEEE International Conference on Big Data (Big Data), 978-986. IEEE. Bishop, Christopher M. 2006. Pattern recognition and machine learning. springer. De Castelbajac, Côme, Mathieu Ritou, S Laporte, and Benoît Furet. 2014. "Monitoring of distributed defects on HSM spindle bearings." Applied Acoustics 77: 159-168. El-Mounayri, Hazim, and Haiyan Deng. 2010. "A generic and innovative approach for integrated simulation and optimisation of end milling using solid modelling and neural network." International Journal of Computer Integrated Manufacturing 23 (1): 40-60. Feng, Guodong, Chunyan Lai, and Narayan C Kar. 2016. "Expectation-maximization particlefilter-and Kalman-filter-based permanent magnet temperature estimation for PMSM condition monitoring using high-frequency signal injection." IEEE Transactions on Industrial Informatics 13 (3): 1261-1270. Godreau, Victor, Mathieu Ritou, Etienne Chové, Benoit Furet, and Didier Dumur. 2019. "Continuous improvement of HSM process by data mining." Journal of Intelligent Manufacturing 30 (7): 2781-2788. Huang, Yonghong, Kevin B Englehart, Bernard Hudgins, and Adrian DC Chan. 2005. "A Gaussian mixture model based classification scheme for myoelectric control of powered upper limb prostheses." IEEE Transactions on Biomedical Engineering 52 (11): 1801-1811. Kim, Dong-Hyeon, Thomas JY Kim, Xinlin Wang, Mincheol Kim, Ying-Jun Quan, Jin Woo Oh, Soo-Hong Min, et al. 2018. "Smart machining process using machine learning: A review and perspective on machining industry." International Journal of Precision Engineering and Manufacturing-Green Technology 5 (4): 555-568. Kohavi, Ron, et al. 1995. "A study of cross-validation and bootstrap for accuracy estimation and model selection." In Ijcai, Vol. 2, 1137-1145. Montreal, Canada. Lenz, Juergen, Thorsten Wuest, and Engelbert Westkämper. 2018. "Holistic approach to machine tool data analytics." Journal of manufacturing systems 48: 180-191. Liu, Chao, Hrishikesh Vengayil, Ray Y Zhong, and Xun Xu. 2018. "A systematic development method for cyber-physical machine tools." Journal of manufacturing systems 48: 13-24. Maleki, Elaheh, Farouk Belkadi, Mathieu Ritou, and Alain Bernard. 2017. "A tailored ontology supporting sensor implementation for the maintenance of industrial machines." Sensors 17 (9): 2063. Maugis, Cathy, Gilles Celeux, and Marie-Laure Martin-Magniette. 2009. "Variable selection for clustering with Gaussian mixture models." Biometrics 65 (3): 701-709. Morgan, Jeff, and Garret E ODonnell. 2018. "Cyber physical process monitoring systems." Journal of Intelligent Manufacturing 29 (6): 1317-1328. Quintana, Guillem, and Joaquim Ciurana. 2011. "Chatter in machining processes: A review." International Journal of Machine Tools and Manufacture 51 (5): 363-376. Reynolds, Douglas. 2015. "Gaussian mixture models." Encyclopedia of biometrics 827-832.

These classifications will be performed focusing on a variable X (e.g., feed rate V f , spindle speed N or spindle power P ) and also on the variation of X, noted ∆X (e.g., ∆V f or ∆N ), which is defined by the simple derivative:

where, X n is the n th measure; X n-1 is the measure at the previous instant and ∆t is the data recording period (0.1s). The objective is to calculate a classification threshold T ∆X by unsupervised machine learning methods. A previous study [START_REF] Wang | Comparison of K-means and GMM methods for contextual clustering in HSM[END_REF] compared the performance of GMM with the k-means method; GMM performs better in this context of research work. The classifications are then carried out based on the following rules:

-There are three clusters for machine-tool movements: stopped (V f = 0), moving at constant speed (V f > 0 and ∆V f ≤ T ∆V f ), or varying speed (V f > 0 and ∆V f > T ∆V f ); -There are three clusters for spindle rotation: stopped (N = 0); rotating at con-


Maximilian Aigner

Valérie Chavez-Demoulin

Armelle Guillou

Measuring and comparing risks of different types

Keywords: Asymptotically unbiased estimator, β-mixing, Convergence in distribution, Expected Proportional Shortfall. 1

Being able to compare risk measures in practice is crucial in many applications such as in finance, insurance or environmental science. The difficulty is that the variables of interest are not always of the same nature, nor of the same type or scale. Thus the usual risk measures are often misleading and to solve this issue we propose to use the Expected Proportional Shortfall (EPS) which is scale invariant and thus, which does not depend on the unit of measurement.

To estimate the EPS, an estimator of the tail index γ is required. The main asymptotic properties of our EPS estimator are provided under very general assumptions in case of d-variate β-mixing processes with Pareto-type marginals. Then, we propose a test statistic based on the EPS estimators to compare different risks, whatever their nature/type/scale are. Since the performances of the test statistic are poor when a biased estimate of γ is used, we propose to perform our EPS estimation with an asymptotically unbiased estimator for γ.

The efficiency of our test statistic is illustrated in a simulation experiment and validated on an environmental dataset.

Introduction and motivation

In many fields where assessing risk matters, it is of main importance to compare several risk measures but often these measures are calculated from variables of different types or nature or have different unit scales. An insurance company, for instance, may be interested in comparing its claims linked to fire and its claims related to car accidents. In finance, one may want to compare market risks of assets of different currencies. In this paper, we consider an environmental application, where the question of interest is: what is the most risky climate extremes, floods or heatwaves, in a certain area? Answering this question is crucial for policy-makers and insurers to prepare for a future of climate change. Although humans can prevent the effects of floods by building dams, extreme rainfalls still strike Europe, destroying constructions, lands and, sometimes, killing people. Likewise, severe heatwaves can lead to fatalities. The European heatwave of 2003, for instance, resulted in at least 30 000 deaths in total, with more than 14 000 recorded in France. These two risks linked to heavy rainfall and heatwaves originate in phenomena of different natures and measured on different scales. To help in addressing the question above, namely, which is the higher risk, we use daily data of precipitation (hourly sum over day, in mm), of temperatures (daily minimum, in degree Celsius) and of humidity at 2m above ground (daily maximum, in %) observed in Lausanne, Switzerland, during three years (2016)(2017)(2018)(2019).

In this context of comparing risks linked to different types of data, the use of traditional risk measures such as the Value-at-Risk (VaR), commonly used in finance, or its equivalent used in environment (return level), can be misleading. To solve this issue, the Expected Proportional Shortfall (EPS) has been introduced by [START_REF] Belzunce | Comparison of risks based on the expected proportional shortfall[END_REF], and is defined for a non-negative random variable X with finite mean and continuous distribution function F X , as

EPS p = E X -VaR p VaR p + , (1)
for all p ∈ {p ∈ (0, 1) : VaR p > 0}, where VaR p = inf{x ∈ R : F X (x) ≥ 1 -p} and a + = max(a, 0). The EPS measures the thickness of the upper tail to a fixed VaR, proportional to the VaR, and thus it can be useful to measure risks of different nature. In particular, it has the nice property of being unaltered when all the values are multiplied by a positive scalar value.

In this paper, we consider d-variate (X (1) , . . . , X (d) ) T time series where each marginal distribution function is of Pareto-type with tail index γ j < 1, j = 1, . . . , d, respectively. We further assume that the vector is coming from a β-mixing process as we typically encounter in environmental applications, often dealing with auto-correlated and cross-correlated components.

In a recent paper, [START_REF] Hoga | Detecting tail risk differences in multivariate time series[END_REF] constructs tests for equality of "tail risk" in the multivariate context of β-mixing processes. The test is applied to financial data with an estimate of the test statistic based on the tail index estimator proposed by [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]. As is well-known in the literature, this estimator may suffer from bias in some circumstances, which affects the performance of the test. To illustrate this point, we simulate 1 000 independent bivariate samples of size 1 000 (corresponding to the size of observations in our application) where each marginal distribution function F j , j = 1, 2, follows a Burr(λ j , τ j ) distribution, namely:

F j (x) = 1 -(1 + x τ j ) -λ j , (2)
where γ j = 1 τ j λ j , j = 1, 2. In our simulations below, we set γ 1 = γ 2 = 0.2, mimicking the estimated values obtained in our application (Section 5). Our objectif is to test the null hypothesis: H 0 : γ 1 = γ 2 . To this aim, we use the test statistic proposed by [START_REF] Hoga | Detecting tail risk differences in multivariate time series[END_REF] with γ j , j = 1, 2, estimated by the Hill estimator (as done in the latter paper), but also with the bias-corrected estimator proposed by [START_REF] Chavez-Demoulin | Extreme quantile estimation for β-mixing time series and applications[END_REF]. The choice of the tuning parameters involved in these estimators is the one suggested in each of the two abovementioned papers.

Since we want to point out that the bias in the tail index estimation has an impact on the performance of the test, an important parameter in this illustration is the second order parameter, denoted subsequently by ρ j , and which is equal to ρ j = -1 λ j , j = 1, 2, for a Burr-type distribution as defined in (2). More precisely, the closer |ρ j | is to 0, the more biased is the Hill estimator.

Table 1 shows the rejection frequency of the null hypothesis H 0 : γ 1 = γ 2 . We choose, ρ 1 = ρ 2 = -1, corresponding to a situation where the Hill estimator does not suffer from bias. The columns of the table correspond to the test statistic proposed by [START_REF] Hoga | Detecting tail risk differences in multivariate time series[END_REF] in case where the tail index is estimated by the Hill estimator (left, Hill-estimate) and the Chavez-Demoulin and This was expected since the second order parameters ρ j , j = 1, 2, are far from 0, and thus the test statistic is not too much sensitive on the tail index estimation method used.

γ 1 = γ 2 in case γ 1 = γ 2 = 0.
2 and ρ 1 = ρ 2 = -1, based on 1 000 samples of size 1 000.

Similarly, Table 2 shows also the rejection frequency of the null hypothesis, but this time, when ρ 1 = -0.5 and ρ 2 ∈ {-0.7, -0.6, -0.5, -0.4, -0.3}.

0 : γ 1 = γ 2 in case γ 1 = γ 2 = 0.
2 and ρ 1 = -0.5, based on 1 000 samples of size 1 000.

Since the bias in estimating the tail index parameter clearly affects the performance of the test statistic, we only use throughout the paper the bias-corrected estimator proposed by [START_REF] Chavez-Demoulin | Extreme quantile estimation for β-mixing time series and applications[END_REF]. First, we extend the theoretical results of the latter paper to the multivariate context of d-variate β-mixing processes. Then, since a test based on equality of tail indices is difficult to interpret for practitioners, except if they are familiar with extreme value theory, we focus in this paper on a test statistic based on the EPS which, in addition to be easily interpretable, allows us to compare risks of different nature, and as such it provides an answer to our initial environmental question.

The remainder of the paper is organized as follows. In Section 2, an estimator of the Expected Proportional Shortfall is proposed and its asymptotic properties are established. Then, in Section 3, a test statistic based on this risk measure has been constructed in order to compare risks of different nature or different unit scales. The performance of our test statistic in terms of size and power is illustrated in Section 4 through a small simulation study. Finally, Section 5 is devoted to the analysis of the real data presented in the introduction in order to answer our environmental question. All the proofs are postponed to the appendix.

2 Estimation of the Expected Proportional Shortfall

Let {ℵ i = (X (1) i , • • • , X (d) i ) T } i=1,•••
,n be a d-variate time series, where the continuous distribution function of X (j) i , j ∈ {1, • • • , d}, are denoted by F j and is assumed to be of Pareto-type with index γ j > 0, that is 1 -F j (x) = x -1/γ j j (x), where j (.) is a slowly varying function at infinity. Define also the associated tail quantile function by U j (x) = inf{y : F j (y) ≥ 1 -1/x}. Our estimator of the Expected Proportional Shortfall is based on an intermediate sequence k, i.e., an integer sequence such that k = k n → ∞ with k/n → 0 as n → ∞. To establish its convergence towards some Gaussian process, we need to assume, as usual in extreme value theory, the following second order condition (SOC):

Assumption (SOC). For j ∈ {1, • • • , d}, there exist ρ j < 0 and a positive or negative function A j (•) with lim t→∞ A j (t) = 0, such that for all y > 0, lim t→∞ U j (ty)

U j (t) -y γ j A j (t) = y γ j y ρ j -1 ρ j , for y > 0.
The function A j (•) in Assumption (SOC) is regularly varying with index ρ j < 0, and the bias of the tail estimators, e.g., extreme value index or EPS estimators, can be explained by the rate of convergence of this function A j (•) towards 0: the closer ρ j is to 0, the larger the bias.

More precisely, our time series {ℵ i } i=1,••• ,n are assumed to be strictly stationary sequences of R d -valued random vectors, which are β-mixing, i.e., such that

β( ) := sup m∈N E sup A∈F ∞ m+ +1 |P(A|F m 1 ) -P(A)| -→ 0, as → ∞,
where F ∞ m+ +1 and F m 1 are the σ-algebra generated by (

ℵ m+ +1 , ℵ m+ +2 , • • • ) and (ℵ 1 , • • • , ℵ m ), respectively.
All our results require some additional assumptions which can be viewed as multivariate analogue of conditions required in extreme value theory for dependent data, as already used in, e.g., [START_REF] Drees | Weighted approximations of tail processes for β-mixing random variables[END_REF][START_REF] Drees | Extreme quantile estimation for dependent data, with applications to finance[END_REF], and defined as follows:

• Assumption (C 1 ). {ℵ i } i∈N is a strictly stationary β-mixing process with continuous marginals and mixing coefficients β(•) such that

lim n→∞ n r n β( n ) + r n √ k log 2 (k) = 0,
for sequences { n } n∈N ⊂ N, {r n } n∈N ⊂ N tending to infinity with n = o(r n ) and r n = o(n);

• Assumption (C 2 ). There exists a function r i,j (x, y), (i, j) ∈ {1, • • • , d} 2 , such that for some δ > 0, we have for all (x, y) ∈ (0, 1 + δ] 2 :

lim n→∞ n r n k Cov rn =1 1l {X (i) >U i ( n kx )} , rn m=1 1l {X (j) m >U j ( n ky )} = r i,j (x, y); • Assumption (C 3 ). For some constants C > 0, δ > 0, we have for i ∈ {1, • • • , d}: n r n k E rn =1 1l {U i ( n ky )<X (i) ≤U i ( n kx )} 4 ≤ C(y -x), ∀ 0 < x < y ≤ 1 + δ, n ∈ N.
Examples of time series satisfying Assumptions (C 1 ) -(C 3 ) can be found, for instance, in [START_REF] Drees | Weighted approximations of tail processes for β-mixing random variables[END_REF][START_REF] Drees | Extreme quantile estimation for dependent data, with applications to finance[END_REF], and include, among others, linear and ARCH models.

Our aim is to estimate the Expected Proportional Shortfall defined in (1). To this goal, observe that the EPS based on X (j) , j ∈ {1, • • • , d}, can be rewritten, for γ j < 1, as follows

EPS p,j = ∞ VaR p,j x VaR p,j -1 dF j (x) = p ∞ 1 1 -F j (zVaR p,j ) 1 -F j (VaR p,j ) dz = p ∞ 1 z -1/γ j dz + A j 1 p ∞ 1 z -1/γ j z ρ j /γ j -1 γ j ρ j dz +A j 1 p ∞ 1   1-F j (zVaR p,j ) 1-F j (VaR p,j ) -z -1/γ j A j 1 p -z -1/γ j z ρ j /γ j -1 γ j ρ j   dz    = p γ j 1 -γ j + O A j 1 p , (3)
since, according to de Haan and Ferreira (2006, p. 161), under our Assumption (SOC), we have the following inequality: for any ε, δ > 0, there exists s 0 = s 0 (ε, δ) > 1 such that for all s, sz ≥ s 0 and j ∈ {1, • • • , d}:

1-F j (zs) 1-F j (s) -z -1/γ j A j (1/(1 -F j (s))) -z -1/γ j z ρ j /γ j -1 γ j ρ j ≤ εz -1 γ j + ρ j γ j max z δ , z -δ . (4)
Thus, to estimate the Expected Proportional Shortfall, we need to estimate the extreme value index γ j .

Chavez-Demoulin and Guillou (2018) introduced a class of estimators for the extreme value index in the case of a univariate β-mixing time series, say X. This class can be viewed as statistical tail functionals, T (Q n ), where Q n is the tail quantile function defined as Q n (u) :=

X n-ku ,n , 0 < u < n/k, and T (•) is a suitable functional, see also [START_REF] Goegebeur | Asymptotically unbiased estimation of the coefficient of tail dependence[END_REF].

We propose to extend this class of estimators in the case of d-variate β-mixing time series, as follows:

γ K j,k (t) := 1 0 log Q (j) n,t (z) Q (j) n,t (1) d(zK(z)), (5)
where

Q (j) n,t (z) := X (j) nt -kt z , nt , for 0 < z < n/k, t ∈ [0, 1], and j ∈ {1, • • • , d}, (6)
and where K satisfies the following assumption:

Assumption (K). Let K be a function such that 1 0 K(t)dt = 1. Suppose that K is continuously differentiable on (0, 1) and that there exist M > 0 and τ ∈ [0, 1/2) such that |K(t)| ≤ M t -τ .

These conditions on the function K(.) are not restrictive, in particular they are satisfied by the usual power kernel function K(u) = (1 + ν)u ν , ν ≥ 0, and the log-weight function K(u) = (-log u) ν /Γ(1 + ν), ν ≥ 0. Note also that the usual tail quantile function Q n has been slightly generalized in (6) since, later in Section 3, we use the self-normalization approach in order to construct the test statistic, and thus we need to evaluate this function for different values of t.

Clearly, the estimator proposed in Chavez-Demoulin and Guillou ( 2018) is a particular case of (5) corresponding to the specific value t = 1, whereas the well-known Hill estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] corresponds to the case K(z) = 1 and t = 1.

Before stating our first theorem on the convergence in distribution of our kernel-type estimator for the index, we first define, for any t 0 ∈ (0, 1), the space D[t 0 , 1] of càdlàg functions on [t 0 , 1] equipped with the Skorohod metric, and by

D d [t 0 , 1] = D[t 0 , 1] × • • • × D[t 0 , 1]
the product space endowed with the product metric.

Theorem 1. Let {ℵ i } i=1,••• ,n be a time series satisfying Assumptions (C 1 ) -(C 3 ) together with the second order condition (SOC) with √ kA j ( n k ) → λ j ∈ R for j ∈ {1, • • • , d}, as n → ∞.
Assume that the distributions F j are continuous, for all j ∈ {1, • • • , d}, and that the kernel K satisfies Assumption (K). Then we have

t √ k      γ K 1,k (t) -γ 1 -A 1 n k 1 0 s -ρ 1 -1 ρ 1 d(sK(s)) . . . γ K d,k (t) -γ d -A d n k 1 0 s -ρ d -1 ρ d d(sK(s))      d -→ Σ 1/2 d B d (t) in D d [t 0 , 1],
where B d (•) denotes a standard d-dimensional Brownian motion and

Σ d := γ i γ j 1 0 1 0 σ i,j (x, y) xy - σ i,j (x, 1) x - σ i,j (1, y) y + σ i,j (1, 1) d(xK(x))d(yK(y)) (i,j)∈{1,••• ,d} 2 , with σ i,j (x, y) := 1 2 [r i,j (x, y) + r i,j (y, x)] , (i, j) ∈ {1, • • • , d} 2 .
This theorem can be viewed as a generalization of Theorem 1 in Chavez-Demoulin and Guillou (2018) to the multivariate context. Recall that our aim is to construct asymptotically unbiased estimators. To reach this goal, we have to use a mixture between two kernels, of the form

∆K 1 (•) + (1 -∆)K 2 (•)
, where ∆ ∈ R is selected in such a way that the bias of our kernel-type estimator computed with this mixture is 0. This yields to a class of asymptotically unbiased estimators. Although the question of finding the asymptotically unbiased estimator with minimal variance is open in the case of a β-mixing time series, in the independent and identically distributed (iid) framework it has been solved by [START_REF] Goegebeur | Asymptotically unbiased estimation of the coefficient of tail dependence[END_REF], who found the following "optimal" kernel

K ρ (z) := 1 -ρ ρ 2 - (1 -ρ)(1 -2ρ) ρ 2 z -ρ , z ∈ (0, 1).
In the sequel, we propose to use this advocated kernel for our β-mixing time series, as already done in Chavez-Demoulin and Guillou (2018) in the case of a univariate β-mixing time series.

From (3), a natural estimator of the Expected Proportional Shortfall in case p → 0 is given by EPS

ρ j p,j (t) = p γ Kρ j j,k (t) 1 -γ Kρ j j,k (t) , for j ∈ {1, • • • , d}. (7)
We are now able to show the convergence in distribution of our EPS estimator defined in (7).

Corollary 1. Let {ℵ i } i=1,••• ,n be a time series satisfying Assumptions (C 1 ) -(C 3 ) together with the second order condition (SOC) with √ kA j ( n k ) → λ j ∈ R for j ∈ {1, • • • , d}, as n → ∞.
Assume that the distributions F j are continuous with γ j < 1, for all j ∈ {1, • • • , d}, and that

p = p n is such that k/(np) → ∞. Then we have t √ k      EPS ρ 1 p,1 (t) EPS p,1 -1 . . . EPS ρ d p,d (t) EPS p,d -1      d -→ Σ 1/2 d B d (t) in D d [t 0 , 1],
where

Σ d := 1 (1-γ i )(1-γ j ) 1 0 1 0 σ i,j (x,y) xy - σ i,j (x,1) x - σ i,j (1,y) y + σ i,j (1, 1) d(xK ρ i (x))d(yK ρ j (y)) (i,j)∈{1,••• ,d} 2 ,
and σ i,j (x, y), (i, j) ∈ {1, • • • , d} 2 , defined as in Theorem 1.

Our Expected Proportional Shortfall estimator is thus asymptotically unbiased in the sense that the expectation of the limit is zero whatever λ j , j ∈ {1, • • • , d}, is.

Unfortunately the second order parameters ρ j , j ∈ {1, ..., d}, are unknown in practice and thus we can either replace them in ( 7) by canonical values or estimators. First, we consider the case where ρ j is replaced by some canonical value ρ j , j ∈ {1, • • • , d}. In that case, the convergence in distribution of our EPS estimators given in Corollary 1 is modified as follows.

Corollary 2. Let {ℵ i } i=1,••• ,n be a time series satisfying Assumptions (C 1 ) -(C 3 ) together with the second order condition (SOC) with √ kA j ( n k ) → λ j ∈ R for j ∈ {1, • • • , d}, as n → ∞.
Assume that the distributions F j are continuous with γ j < 1, for all j ∈ {1, • • • , d}, and that

p = p n is such that k/(np) → ∞. Then we have t √ k       EPS ρ 1 p,1 (t) EPS p,1 -1 -A 1 n k (1-ρ 1 )( ρ 1 -ρ 1 ) ρ 1 (1-ρ 1 )(1-ρ 1 -ρ 1 )(1-γ 1 )γ 1 . . . EPS ρ d p,d (t) EPS p,d -1 -A d n k (1-ρ d )( ρ d -ρ d ) ρ d (1-ρ d )(1-ρ d -ρ d )(1-γ d )γ d       d -→ Σ 1/2 d B d (t) in D d [t 0 , 1],
where

Σ d := 1 (1-γ i )(1-γ j ) 1 0 1 0 σ i,j (x,y) xy - σ i,j (x,1) x - σ i,j (1,y) y + σ i,j (1, 1) d(xK ρ i (x))d(yK ρ j (y)) (i,j)∈{1,••• ,d} 2 ,
and σ i,j (x, y), (i, j) ∈ {1, • • • , d} 2 , defined as in Theorem 1.

One clearly loses the bias correction, except if ρ j corresponds to the true value of ρ j , j ∈ {1, • • • , d}. However, the Expected Proportional Shortfall estimators are not very sensitive to such a mis-specification and thus our estimator EPS ρ j p,j (t), j ∈ {1, ..., d}, can still outperform the estimators that are not corrected for bias. If we want nevertheless to keep the property of being asymptotically unbiased, we need to estimate ρ j in a consistent way. This leads to the following general result.

Corollary 3. Let {ℵ i } i=1,••• ,n be a time series satisfying Assumptions (C 1 ) -(C 3 ) together with the second order condition (SOC) with √ kA j ( n k ) → λ j ∈ R for j ∈ {1, • • • , d}, as n → ∞.
Assume that the distributions F j are continuous with γ j < 1, for all j ∈ {1, • • • , d}, and that p = p n is such that k/(np) → ∞. Let ρ j,kρ := ρ j be an external estimator for ρ j , j ∈ {1, • • • , d}, consistent in probability, which depends on an intermediate sequence k ρ . Then we have

t √ k       EPS ρ 1 p,1 (t) EPS p,1 -1 . . . EPS ρ d p,d (t) EPS p,d -1       d -→ Σ 1/2 d B d (t) in D d [t 0 , 1],
where Σ d is defined as in Corollary 1.

A possible choice for ρ j,kρ , j ∈ {1, • • • , d}, is that proposed by [START_REF] Gomes | Semi-parametric estimation of the second order parameter in statistics of extremes[END_REF] and defined as ρ j,kρ := -4 + 6 S

(2)

j,kρ + 3 S (2) j,kρ -2 4 S (2) j,kρ -3 provided S (2) j,kρ ∈ 2 3 , 3 4 ,
where

S (2) j,kρ := 3 4 M (4) j,kρ -24 M (1) j,kρ 4 M (2) j,kρ -2 M (1) j,kρ 2 M (3) j,kρ -6 M (1) j,kρ 3 2 with M (α) j,kρ := 1 k ρ kρ =1 log X (j) n-+1,n -log X (j) n-kρ,n α , α ∈ N.
According to Theorem 2.1 in [START_REF] Gomes | Semi-parametric estimation of the second order parameter in statistics of extremes[END_REF], this estimator ρ j,kρ is consistent in probability as soon as the intermediate sequence k ρ satisfies k ρ A j (n/k ρ ) → ∞ and our second order condition (SOC) hold.

Test statistic based on the Expected Proportional Shortfall

We have now all the ingredients to move to our final goal of the paper, i.e., the construction of our test statistic, based on the Expected Proportional Shortfall, in order to compare risks of different nature or of different unit scales. Our null hypothesis is formulated as follows

H EPS 0 : lim p→0 EPS p,j EPS p,d = 1, or equivalently lim p→0 log EPS p,j EPS p,d = 0, ∀j = 1, • • • , d -1,
and allows us to compare the EPS of β-mixing series as p goes to zero.

To construct our test statistic we propose to use the self-normalization approach, originally introduced by [START_REF] Kiefer | Simple robust testing of regression hypotheses[END_REF], [START_REF] Lobato | Testing that a dependent process is uncorrelated[END_REF] and [START_REF] Shao | A self-normalized approach to confidence interval construction in time series[END_REF]. See also [START_REF] Shao | Self-normalization for time series: a review of recent developments[END_REF] for a recent review. The idea behind this approach is to use a ratio-type statistic in order to cancel the unknown asymptotic variance Σ d appearing in both the numerator and denominator. Thus a pivotal limiting distribution results, whose critical values may be computed through simulations. The main advantages of this approach is the simplicity of implementation together with the fact that it does not depend on any tuning parameters.

To this aim, define the two (d -1) × 1 matrices:

LREPS p (t) :=         log EPS ρ 1 p,1 (t) EPS ρ d p,d (t) . . . log EPS ρ d-1 p,d-1 (t) EPS ρ d p,d (t)         and LREPS p (t) :=         1 d d j=1 log EPS ρ j p,j (t) EPS ρ d p,d (t) . . . 1 d d j=1 log EPS ρ j p,j (t) EPS ρ d p,d (t)         .
Then, our test statistic for H EPS 0 is given by

T EPS p := LREPS p (1) -LREPS p (1) T × 1 t 0 t 2 LREPS p (t) -LREPS p (t) -LREPS p (1) -LREP S p (1) × LREPS p (t) -LREPS p (t) -LREPS p (1) -LREPS p (1) T dt -1 × LREPS p (1) -LREPS p (1) , (8)
and its behaviour, under the null hypothesis and its alternative, is summarized in the following theorem.

Theorem 2. Under the same conditions as in Corollary 3, if Σ d is a positive-definite matrix, then we have:

• under H EPS 0 : T EPS p d -→ B T d-1 (1) 1 t 0 {B d-1 (t) -tB d-1 (1)} {B d-1 (t) -tB d-1 (1)} T dt -1 B d-1 (1), • under H EPS 1 : T EPS p k = O P (1).

Simulation experiment

In this section, we perform a simulation study in order to illustrate the performance in terms of size and power of our test statistic defined in (8). We limit ourselves to the case d = 2, although similar results can be obtained in higher dimensions. Note that our EPS estimator ( 7) is defined in case p tends to 0, but since our test statistic involves ratios between two EPS estimators, it no longer depends on p. In practice this means that the conclusion of the test is stable over levels p (small enough) or equivalently over time horizons if the interpretation of the EPS is in terms of the return period. This seems realistic in the environmental context: If, in a given area and over a certain (long-term) horizon, flood is more risky than heatwave, this is likely to be true for higher time horizon, under same underlying conditions of climate or human intervention. We choose throughout this section p = 1/n corresponding to the boundary of observations. Now, we need to select t 0 . This parameter cannot be too large or too small. Indeed, a small value of t 0 means that the EPS estimator ( 7) is based on samples that are too small for t close to t 0 , whereas a too large value, implies that the integral-term in our test statistic ( 8) is based on too few differences. An extensive simulation study (not reported in the paper) shows that the value t 0 = 0.3 provides a reasonable trade-off and performs well in all the models considered and whatever the sample size is. Then, we have to compute the critical values of the limiting distribution in Theorem 2 for t 0 = 0.3. This is done by a Monte Carlo method using 100 000 paths of standard Brownian motions, each path coming from 100 000 standard normal random variables. Table 3 gives these critical values.

Next, the intermediate sequence k used in our test statistic needs to be selected. To this aim,

Critical value 67 41

Table 3: Critical values of the limiting distribution in Theorem 2 for t 0 = 0.3 and d = 2, for two levels.

the following algorithm is used. For a given series j, j = 1, ..., d, since the Chavez-Demoulin and Guillou (2018) estimator for γ is asymptotically unbiased, we expect to see a plateau when we look at this estimate as a function of the intermediate sequence. Thus, the latter is selected as the smallest value from which the median of the estimates of γ over all the replications exhibit a stable part. We denote this value by k j , j = 1, ..., d. Finally, the retained value of k used in our test statistic corresponds to k = min 1≤j≤d k j . We use here the minimum since a too large value of this parameter might induce a bias in the estimation procedure.

Concerning the choice of the sequence k ρ required in the estimation of the second order parameter In this simulation study, we use two data-generating β-mixing processes. One corresponds to the iid case of two Burr-type random variables as defined in (2) and the second one, to a VMA process, which is a more realistic process encountered in real environmental context.

More precisely, it is a β-mixing process that allows both auto-correlation and cross-correlation, and thus it can naturally appear in environmental data, auto-correlation reflecting persistence of weather, whereas cross-correlation reflects the dynamics of weather. Our VMA(1) process considered is of the following form

   X (1) i = W (1) i + W (1) i-1 + νW (2) i-1 X (2) i = W (2) i + W (2) i-1 , (9)
i = 1, . . . , n, where (W (1) , W (2) ) T , is a bivariate random vector with margins W (j) ∼ Burr(λ j , τ j ), j = 1, 2 and with dependence structure coming from a Gumbel copula with parameter θ = 2.

The choice of a Gumbel copula is to exhibit upper tail dependence. We take ν = -0.45 to mimic the (negative) cross correlation we observe between temperature and humidity in our application. with γ 1 = 0.2, mimicking the estimated values obtained in our application and with ρ 1 = -1

and ρ 2 = -0.25. The index γ 2 varies from 0.2 to 0.9. We run 1 000 samples of size 500, 1 000 and 2 000.

From these two tables, one can draw the following conclusions:

• If γ 2 = γ 1 , the frequencies of rejection, which correspond to the type I error, are rather close to the level of the test, especially in Table 4 when the sample size is sufficiently large;

• If γ 2 is different from γ 1 , but without being too far away from it (e.g., γ 1 = 0.2 and γ 2 = 0.25), the frequencies of rejection are already high, especially in Table 5. It underlines the quality of the test which is able to detect deviations, even small, with the null hypothesis;

• When γ 2 deviates significantly from γ 1 , the frequencies of rejection drastically increase for both levels and all sample sizes, illustrating the high power of the test.

Application to an environmental dataset

The past few decades have seen extreme climate events with catastrophic impacts on human society and the economy. In this study, our aim is to compare the different risks of observing extreme precipitation, extreme humidity and extreme daily maximum temperature in Lausanne. The three being measured on different scales, it does not make sense to compare their VaR whereas the EPS is precisely an adequate risk measure in this context. Comparing the tail index γ of the three time series would be an alternative but unfortunately it is not easily interpretable for practitioners or people who are not familiar with extreme value theory. On the other hand, the EPS can claim to be a universal measure for risks of different natures.

Our dataset concern 1003 observations measured daily in Lausanne, from July 9, 2016 to April 12, 2019. The data consist of humidity at 2m above ground, in %, daily maximum (Figure 1); precipitation, in mm, hourly sum over day (Figure 2); and temperature, in degree Celsius, daily maximum (Figure 3). ACF plots (not shown here) of the series show some auto-correlation, so we apply ARMA models to the series individually and keep the residuals series to work with.

According to the CCF plots of Figure 4, the residuals look cross-correlated, especially when the residual series of precipitation is involved, which means that a VMA process can be adjusted to our data. Figure 5 represents the estimators for γ, with, on the left panel, the Chavez-Demoulin and Guillou (2018) estimators and, on the right panel, the Hill estimators, for temperature (solid line), humidity (dotted line) and precipitation (dashed line). In the three cases, the estimates of γ based on the Chavez-Demoulin and Guillou (2018) method look more stable as a function q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 40 50 60 70 get EPS H = 2.79 × 10 -4 , for precipitation we have EPS P = 5.38 × 10 -4 , and for temperature, EPS T = 7.78 × 10 -5 . As expected, and in accordance with climate science literature, the risk of extreme rainfall in this region is higher than the risk of high humidity or high daily temperature. Lausanne has been struck by huge spring precipitation in June 2018 and 2019.

On June 11th, 2018, 41.0 mm fell within 10 minutes, making that period of three days (three crossed points in Figure 2) an extreme climate event. Although the risk of humidity is lower than the risk of high precipitation, it is linked to precipitation by physical law. Table 6 provides the values of the test statistic based on the comparison of the three series together and then of all combinations of pairs. According to Table 3, the critical value corresponding to a level of 5% is 67 in dimension 2 and a similar computation yields to 168 in dimension 3. Since we reject the null hypothesis when using the three series, it is important to perform pairwise tests. We notice that we reject the null hypothesis based on the EPS of the pair (Precipitation, Temperature) whereas we do not reject the null hypothesis when we use the pairs (Humidity, Precipitation) and (Humidity, Temperature). As shown in the left panel of Figure 5, the estimator for γ based on the precipitation variable is larger than the ones for humidity and temperature. This means that we can expect more rainfall extremes than high temperature or humidity extremes in Lausanne. The fact that we get significantly much higher risk (EPS) of extreme precipitations than temperature is in accordance with the increase of extreme rainfalls recently observed in Switzerland and planned by the climate extremes projections. Although there is an increase of heatwaves observed in Switzerland, note that the daily maximum temperature does not really define an heatwave but just a punctual peak, the definition of an heatwave including also the heat duration. For the pair (Humidity, Precipitation), the non-rejection of the null hypothesis may be due to the physical link between these two variables, whereas for the pair (Humidity, 6 Appendix: Proofs of the results

Before proving our Theorem 1, we first remark that our kernel-type estimator of the index, γ K j,k (t), can be rewritten as where

γ K j,k (t) = 1 0 log Q (j) n,t (z) Q (j) n,t (1)
d(zK(z)) = kt i=1 i kt i-1 kt log Q (j) n,t (z) Q (j) n,t (1)
d(zK(z)) = kt i=1 log X (j) nt -i+1, nt X (j) nt -kt , nt i kt i-1 kt d(zK(z)) = kt i=1 log X (j) nt -i+1, nt X (j) nt -kt , nt i kt K i kt - i -1 kt K i -1 kt = kt i=1 log X (j) nt -i+1, nt i kt K i kt - i -1 kt K i -1 kt -K(1) log X (j) nt -kt , nt = kt -1 i=0 kt -i kt K kt -i kt log X (j) nt -kt +i+1, nt X (j) nt -kt +i, nt = kt -1 i=0 X (j) nt -kt +i+1, nt X (j) nt -kt , nt X (j) nt -kt +i, nt X (
F n,j (t, y) := 1 kt nt i=1 1l {X (j) i >yX (j) nt -kt , nt } =        kt -i kt if y ∈ X (j) nt -kt +i, nt X (j) nt -kt , nt , X (j) nt -kt +i+1, nt X (j) nt -kt , nt 0 if y ≥ X (j) nt , nt X (j) nt -kt , nt
.

Thus to study the asymptotic properties of γ K j,k (t), we need first to look at the behaviour of F n,j (t, •). This preliminary result is given and proved in Section 6.1 below.

Preliminary result

Theorem A. Let {ℵ i } i=1,••• ,n be a time series satisfying Assumptions (C 1 ) -(C 3 ) together with the second order condition (SOC) with √ kA j ( n k ) → λ j ∈ R for j ∈ {1, • • • , d}, as n → ∞. Assume that the distributions F j are continuous, for all j ∈ {1, • • • , d}. Then, for any t 0 ∈ (0, 1) and ν ∈ [0, 1/2), under a Skorohod construction, we have, as n → ∞ sup

t∈[t 0 ,1],z≥1        z ν γ 1 t √ k F n,1 (t, z) -z -1 γ 1 -A 1 n k z -1 γ 1 z ρ 1 γ 1 -1 γ 1 ρ 1 -W 1 t, z -1 γ 1 -z -1 γ 1 W 1 (t, 1) . . . z ν γ d t √ k F n,d (t, z) -z -1 γ d -A d n k z -1 γ d z ρ d γ d -1 γ d ρ d -W d t, z -1 γ d -z -1 γ d W d (t, 1)        a.s.
-→ 0, where W(., .)

= (W 1 (•, •), • • • , W d (•, •)) T is
(i,j)∈{1,••• ,d} 2 .
Proof of Theorem A. We follow the lines of the proof of Theorem 1 in [START_REF] Hoga | Detecting tail risk differences in multivariate time series[END_REF], in particular, we consider random variables U

(j) i , i ∈ {1, • • • , n}, j ∈ {1, • • • , d}, uniformly distributed on [0, 1], such that • Assumption (U 1 ). U i = (U (1) i , • • • , U (d) i ) T i∈N
is a strictly stationary β-mixing process with uniform[0, 1] marginals and mixing coefficients β(•) such that

lim n→∞ n r n β( n ) + r n √ k log 2 (k) = 0
for sequences { n } n∈N ⊂ N, {r n } n∈N ⊂ N tending to infinity with n = o(r n ) and r n = o(n);

• Assumption (U 2 ). There exists a function r U i,j (u, v), (i, j) ∈ {1, • • • , d} 2 , such that for some δ > 0, we have for all (u, v) ∈ (0, 1 + δ] 2 :

lim n→∞ n r n k Cov rn =1 1l {U (i) >1-k n u} , rn m=1 1l {U (j) m >1-k n v} = r U i,j (u, v);
• Assumption (U 3 ). For some constants C > 0, δ > 0, we have for

i ∈ {1, • • • , d}: n r n k E rn =1 1l {1-k n v<U (i) ≤1-k n u} 4 ≤ C(v -u), ∀ 0 < u < v ≤ 1 + δ, n ∈ N.
Then, as n → ∞, we have, for any ν ∈ [0, 1/2)

y -ν √ k       1 k nt =1 1l {U (1) >1-k n y} -k n y . . . 1 k nt =1 1l {U (d) >1-k n y} -k n y       d -→ y -ν      W 1 (t, y) . . . W d (t, y)      in D d ([t 0 , 1] × (0, 1 + δ]). ( 11
)
Note that Assumptions (C 1 ) -(C 3 ) for X (j) , j ∈ {1, • • • , d}, imply Assumptions (U 1 ) -(U 3 ) for j) ). Also

U (j) = F j (X (
X (j) > U j n ky ⇔ U (j) > 1 - k n y
which implies by ( 11) and since nt -1 < nt ≤ nt, that, as n → ∞,

y -ν √ k       1 k nt =1 1l {X (1) >U 1 ( n ky )} -ty . . . 1 k nt =1 1l {X (d) >U d ( n ky )} -ty       d -→ y -ν      W 1 (t, y) . . . W d (t, y)      in D d ([t 0 , 1] × (0, 1 + δ]).
By the Skorohod's representation theorem, this convergence holds a.s. on a suitable probability space, i.e., sup t∈[t 0 ,1],y∈(0,1+δ]

y -ν √ k       1 k nt =1 1l {X (1) >U 1 ( n ky )} -ty . . . 1 k nt =1 1l {X (d) >U d ( n ky )} -ty       -y -ν      W 1 (t, y) . . . W d (t, y)      a.s.
-→ 0,

the limit being continuous, the convergence is thus uniform. Now, we need to show that U j ( n ky ) can be replaced by U j ( n k )y -γ j , j ∈ {1, • • • , d}. To this aim, we use the inequality (4) by setting s = U j ( n k ), z = y -γ j , from which we deduce that, as n → ∞,

sup y∈(0,1+δ] n k [1 -F j (y -γ j U j ( n k ))] -y -A j ( n k ) y y -ρ j -1 γ j ρ j yA j ( n k ) = o(1). ( 13
) Set y n,j := n k [1 -F j (y -γ j U j ( n k ))]. In the j-th row of the matrix in (12), j = 1, • • • , d, inserting y n,j for y, implies that sup t∈[t 0 ,1],y∈(0,1+δ]        y -ν n,1 √ k 1 k nt =1 1l {X (1) >U 1 ( n ky n,1
)} -ty n,1 -W 1 (t, y n,1 ) . . .

y -ν n,d √ k 1 k nt =1 1l {X (d) >U d ( n ky n,d )} -ty n,d -W d (t, y n,d )        a.s. -→ 0, i.e., sup t∈[t 0 ,1],y∈(0,1+δ]       y -ν n,1 √ k 1 k nt =1 1l {X (1) >y -γ 1 U 1 ( n k )} -ty n,1 -W 1 (t, y n,1 ) . . . y -ν n,d √ k 1 k nt =1 1l {X (d) >y -γ d U d ( n k )} -ty n,d -W d (t, y n,d )       a.s.
-→ 0.

Using ( 13), the fact that

√ kA j ( n k ) → λ j ∈ R and a uniform continuity argument for W j (t, •), j ∈ {1, • • • , d}, we have sup t∈[t 0 ,1],y∈(0,1+δ]       y -ν √ k 1 k nt =1 1l {X (1) >y -γ 1 U 1 ( n k )} -ty -A 1 ( n k )t y y -ρ 1 -1 γ 1 ρ 1 -W 1 (t, y) . . . y -ν √ k 1 k nt =1 1l {X (d) >y -γ d U d ( n k )} -ty -A d n k t y y -ρ d -1 γ d ρ d -W d (t, y)       a.s.
-→ 0, or by changing y into z

-1 γ j (j-th component) sup t∈[t 0 ,1],z≥1-δ        z ν γ 1 √ k 1 k nt =1 1l {X (1) >zU 1 ( n k )} -tz -1 γ 1 -A 1 n k tz -1 γ 1 z ρ 1 γ 1 -1 γ 1 ρ 1 -W 1 t, z -1 γ 1 . . . z ν γ d √ k 1 k nt =1 1l {X (d) >zU d ( n k )} -tz -1 γ d -A d n k tz -1 γ d z ρ d γ d -1 γ d ρ d -W d t, z -1 γ d        a.s. -→ 0, (14)
for some δ > 0. The next step consists to replace U j ( n k ), j ∈ {1, • • • , d}, by suitable order statistics. This can be done since from ( 14) and the generalized Vervaat Lemma of Einmahl et al. (2010, Lemma 5), we have for all j ∈ {1,

• • • , d}, as n → ∞, sup t∈[t 0 ,1] √ k        X (j) nt -kt , nt U j ( n k )   -1 γ j -1      + W j (t, 1) t a.s.
-→ 0.

Finally, replacing in ( 14) z by z n = zX

(j) nt -kt , nt /U j ( n k ) (j-th component) and since kt k = t + O 1
k uniformly in t, Theorem A follows. Note that the covariance structure of the Gaussian process W(•, •) comes from Lemma 1 in [START_REF] Hoga | Detecting tail risk differences in multivariate time series[END_REF].

We have now all the ingredients in order to prove our main results.

Proofs of the main results

Proof of Theorem 1. To prove this theorem, we make use of the Cramér-Wold device (see, e.g., Severini, 2005, p. 337). More precisely, we show that

β T t √ k      γ K 1,k (t) -γ 1 -A 1 n k 1 0 s -ρ 1 -1 ρ 1 d(sK(s)) . . . γ K d,k (t) -γ d -A d n k 1 0 s -ρ d -1 ρ d d(sK(s))      d -→ β T      γ 1 1 0 s -1 W 1 (t, s) -W 1 (t, 1) d(sK(s)) . . . γ d 1 0 s -1 W d (t, s) -W d (t, 1) d(sK(s))      for all β ∈ R d .
To this aim, first remark that, using (10), and since, under Assumption (K), we have

γ j = ∞ 1 y -1/γ j 0 d(zK(z)) dy y , direct computations lead to β T t √ k      γ K 1,k (t) -γ 1 -A 1 n k 1 0 s -ρ 1 -1 ρ 1 d(sK(s)) . . . γ K d,k (t) -γ d -A d n k 1 0 s -ρ d -1 ρ d d(sK(s))      = d =1 β t √ k γ K ,k (t) -γ -A n k 1 0 s -ρ -1 ρ d(sK(s)) = d =1 β t √ k γ 1 0 F n, (t,z -γ ) 0 d(sK(s)) - z 0 d(sK(s)) dz z -A n k 1 0 s -ρ -1 ρ d(sK(s)) .
Then, applying the mean value theorem, we get for some ξ z ∈ (0, 1)

F n, (t,z -γ ) 0 d(sK(s)) = z 0 d(sK(s)) + F n, (t, z -γ ) -z d(sK(s)) ds | s=z+ξz[F n, (t,z -γ )-z] .
This implies that, uniformly in t ∈ [t 0 , 1]:

β T t √ k      γ K 1,k (t) -γ 1 -A 1 n k 1 0 s -ρ 1 -1 ρ 1 d(sK(s)) . . . γ K d,k (t) -γ d -A d n k 1 0 s -ρ d -1 ρ d d(sK(s))      = d =1 β t √ k γ 1 0 F n, (t, z -γ ) -z d(sK(s)) ds | s=z+ξz[F n, (t,z -γ )-z] dz z -A n k 1 0 s -ρ -1 ρ d(sK(s)) d -→ d =1 β γ 1 0 s -1 W (t, s) -W (t, 1) d(sK(s))
by our Theorem A. This achieves the proof of Theorem 1.

Proof of Corollary 1. Combining (3) with ( 7), together with the fact that |A j (.)| is regularly varying at infinity with index ρ j < 0, we deduce that (1 -ρ j )( ρ j -ρ j )

t √ k EPS ρ j p,j (t) EPS p,j -1 = t √ k   γ Kρ j j,k (t) 1 -γ Kρ j j,k (t) 1 -γ j γ j + O(A j ( 1 p )) -1   = 1 γ j 1 -γ Kρ j j,k (
ρ j (1 -ρ j )(1 -ρ j -ρ j )(1 -γ j )γ j   = 1 γ j (1 -γ j ) t √ k γ K ρ j j,k (t) -γ j -A j n k
(1 -ρ j )( ρ j -ρ j ) ρ j (1 -ρ j )(1 -ρ j -ρ j )

+ o P (1), uniformly in t ∈ [t 0 , 1]. Now, combining Theorem 1 with the fact that for all j ∈ {1, • • • , d} 1 0 s -ρ j -1 ρ j d(sK ρ j (s)) = (1 -ρ j )( ρ j -ρ j ) ρ j (1 -ρ j )(1 -ρ j -ρ j ) Corollary 2 follows.

Proof of Corollary 3. Again, following the lines of proof of Corollary 1, we have The convergence in distribution of T 1 follows from Theorem 1. Concerning T 2 , we have to show that each component is uniformly negligible as n → ∞. To this aim, remark that, using (10) combining with the mean value theorem as in the proof of Theorem 1, we have for some F n,j (t, z -γ j ) -z d(s[K ρ j (s) -K ρ j (s)]) ds | s=z+ξz[F n,j (t,z -γ j )-z] dz z since the term in ( 15) is equal to 0. Now, using again the mean value theorem, we have d(s[K ρ j (s) -K ρ j (s)]) ds = ( ρ j -ρ j ) -2(1 -ρ j ) ρ 3 j + 2(1 -ρ j )(1 -ρ j -ρ 2 j ) ρ 3 j s -ρ j +

t √ k   EPS ρ j p,j (t) EPS p,j -1   = t √ k   γ K ρ j j,k (t) 1 -γ K ρ j j,k (t) 1 -γ j γ j + O(A j ( 1 p )) -1   = 1 γ j 1 -γ K ρ j j,k (t) t √ k γ K ρ j j,k (t) -γ j + O P √ kA j 1 p = 1 γ j (1 -γ j ) t √ k γ K ρ j j,k (
ξ z ∈ (0, 1) t √ k γ K ρ j j,k (t) -γ Kρ j j,k (t) = t √ k ∞ 1 F n,j (t,
(1 -ρ j ) 2 (1 -2 ρ j ) ρ 2 j s -ρ j ln s with ρ j an intermediate value between ρ j and ρ j . Since ρ j is a consistent estimator for ρ j , Theorem A achieves the proof of Corollary 3.

  ρ, as in Chavez-Demoulin and Guillou (2018) we use k ρ := sup k : k ≤ min m -1, 2m log log m and ρ k exists , with m the number of positive observations in the sample.

Figure 2 :Figure 3 :

 23 Figure 2: Precipitation in Lausanne (in mm, hourly sum over day) from July 9, 2016 to April 12, 2019.

Figure 4 :

 4 Figure 4: Cross Correlation plot (CCF) between Precipitation-Humidity, Precipitation-Temperature and Humidity-Temperature.

  Chavez-Demoulin and Guillou estimator (left panel) and the Hill estimator (right panel) for temperature (solid line), humidity (dotted line) and precipitation (dashed line). Temperature), it shows an equivalent and limited risk for each variable, in Lausanne. Test statistic based on the three series (Humidity, Precipitation, Temperature) and the pairs (Humidity, Precipitation), (Humidity, Temperature), (Precipitation, Temperature).

  (t) -γ j + o P (1), uniformly in t ∈ [t 0 , 1].Now, using Theorem 1 and the fact that for all j ∈ {1, • • • , d},

  t) -γ j + o P (1), uniformly in t ∈ [t 0 ,

F

  sK ρ j (s)) -d(sK ρ j (s))] sK ρ j (s)) -d(sK ρ j (s))] sK ρ j (s)) -d(sK ρ j (s))] n,j (t, z -γ j ) -z d(s[K ρ j (s) -K ρ j (s)]) ds | s=z+ξz[F n,j (t,z -γ j )-z]

d

  the (d -1) × d matrix defined as can remark that under H EPS 0 , according to (3), we have all the indices equal, i.e., γ 1 = • • • = γ d . This implies that, under H since |A (.)| and |A d (.)| are regularly varying with negative parameters and k/(np) → ∞. B d (s) = min(s, t)AΣ d A T ,

Table 1 :

 1 Burr distribution. Rejection frequency of the null hypothesis H 0 :

		Hill-estimate CG-estimate
	5%	0.055	0.047
	10%	0.12	0.092

Table 2 :

 2 Burr distribution. Rejection frequency of the null hypothesis H

	ρ 2	-0.7	-0.6	-0.5	-0.4	-0.3
	5%
	Hill-estimate 0.036 0.060 0.091 0.146 0.274
	CG-estimate 0.046 0.047 0.052 0.052 0.049
	10%
	Hill-estimate 0.101 0.123 0.160 0.236 0.421
	CG-estimate 0.088 0.100 0.099 0.096 0.110

As |ρ 2 | becomes closer to 0, the rejection frequency based on the test statistic computed with the Hill-estimate increases drastically, leading to an important type I error. On the contrary, with the CG-estimate, the rejection frequency stays close to the level.

Table 4 ,

 4 respectively Table5, shows the rejection frequency of the null hypothesis H EPS

	0	in the

Table 4 :

 4 Rejection frequency of the null hypothesis H EPS

	In addition to global warming, many climate projections suggest an
	increase in climate extremes in Europe (IPCC, 2012). With the warming of the oceans, periods
	of excessive rainfall will appear considerably more often. The European Environment Agency
	(www.eea.europa.eu) supports this projection, endorsing, in November 2019, the statement that
	"Heavy precipitation events are likely to become more frequent in most parts of Europe" and

0 in the iid case of two Burr-type random variables based on 1 000 samples, for different values of γ 2 , levels and sample sizes. The parameters are γ 1 = 0.2, ρ 1 = -1 and ρ 2 = -0.25.

year, what is the highest risk between a heatwave or an heavy precipitation? These two events have alternatively stricken Lausanne several times in late spring season.

  (y 1 , y 2 ) + r i,j (y 2 , y 1 )]

		a d-variate continuous centered Gaussian process
	with covariance function given by
	Cov(W(t 1 , y 1 ), W(t 2 , y 2 )) =	min(t 1 , t 2 ) 2	[r i,j

  Proof of Theorem 2. First remark thatt √ k LREPS p (t) -LREPS p (t) =: A × t

	√	k	      log log	EPS ρ 1 p,1 (t) EPS p,1 . . . EPS ρ d p,d (t) EPS p,d	     	+ t √	k	     log log EPS p,d-1 EPS p,1 EPS p,d -1 d EPS p,d -1 . . . d	d =1 log d =1 log EPS p, EP S p,d EPS p, EPS p,d

Acknowledgement

The research of Armelle Guillou was supported by the French National Research Agency under the grant ANR-19-CE40-0013-01/ExtremReg project and an International Emerging Action (IEA-00179). The authors would like to thank Yannick Hoga for letting the R-code of its paper published in 2018 available.

where AΣ d A T is a (d -1) × (d -1) Hermitian positive-definite matrix, which thus can be decomposed into LL T , where L is a lower triangular matrix with positive diagonal entries. Thus

The first part of Theorem 2 then follows from an application of the continuous mapping theorem.

The second part of Theorem 2 is a direct consequence of our Corollary 3 combining with the fact that, under