Esperan Padonou

Olivier Roustant

Michel Lutz

Robust monitoring of an industrial IT system in presence of structural change

Keywords: Monitoring, Control charts, Holt-Winters, Robust statistics, Structural changes

This paper presents an original research initiated by the monitoring needs of a semiconductor production plant. The industrial operations rely on an Information Technology (IT) system, and several time series data are controlled statistically. Unfortunately, these variables often contain outliers, as well as structural changes due to external decisions in the IT activity. As a consequence, it has been observed that the monitoring results obtained with standard techniques could be severely biased. This paper presents some contributions to overcome such difficulties. A new monitoring method is proposed, based on robust Holt-Winters smoothing algorithm, and coupled with a relearning procedure for structural breaks detection. Such a method is flexible enough for a large-scale industrial application. We evaluate its performances through simulations studies, and show its usefulness in industrial real applications for univariate and multivariate time series. The scope of application deals with IT activity monitoring, but the introduced statistical methods are generic enough for being used in other industrial fields.

Introduction

Manufacturing activities efficiency relies increasingly on Information Technologies (IT). This is clearly true for wafer fabrication plants, where this research is ongoing through a partnership with the company STMicroelectronics. Several activities are useful for a proper IT management (see eg. Rudd et al. [START_REF] Rudd | Service Design, Itil, Version 3[END_REF] ). This paper is focused on the monitoring activity. A careful monitoring is based on a cautious observation of the IT system. The objective is to have a close look on all data recorded and stored to track the activity of a plant IT system.

As the size and complexity of IT systems strongly increase, the viability of monitoring can only be ensured by employing automated procedures as mentioned by Dugmore et al. [START_REF] Dugmore | Capacity management[END_REF] Usual monitoring is based on threshold exceeding detection (Rudd et al. [START_REF] Rudd | Service Design, Itil, Version 3[END_REF] ). However, facing the enormous diversity and expanse of the components of modern IT systems, IT experts do not have always enough knowledge or time to determine a priori the critical thresholds. Consequently, automated monitoring procedures that do not require any preliminary hypothesis about the monitored variables, are required. In this case, the monitoring activity requires the ability to distinguish normal system operations from exceptional events.

To quantitatively detect such abnormal behaviors in the fluctuation of a time series, a standard procedure is the following: 1) Model the time series to be monitored; 2) Monitor model residuals, through an appropriate control chart (see eg. Montgomery [START_REF] Montgomery | Statistical Quality Control: A Modern Introduction[END_REF] , Croux et al. [START_REF] Croux | Robust control chart for time series data[END_REF] ). This procedure is known in IT management studies. By way of example, Hellerstein [START_REF] Hellerstein | Quantitative Data Cleaning for Large Databases[END_REF] has developed an approach based on ARIMA modelling, whereas Brutlag [START_REF] Brutlag | Aberrant Behavior Detection in Time Series for Network Monitoring[END_REF] and Leikis [START_REF] Leikis | Consolidated Capacity and Performance Reporting[END_REF] have recourse to Holt-Winters smoothing. This paper focuses on the latter. Holt-Winters smoothing is well adapted in industrial contexts, where time series have to be modeled without any preliminary analysis. Moreover, it is a fairly good approximation for many kinds of time series, accurate enough for industrial applications (see eg. Makridakis et al. [START_REF] Makridakis | Forecasting, methods and applications[END_REF] , Croux et al. [START_REF] Croux | Robust control chart for time series data[END_REF] ). Furthermore, we observed that the Holt-Winters model is validated for most of the time series encountered in our research field.

However, a monitoring based on the usual Holt-Winters algorithm revealed several limitations at STMicroelectronics. That is why we developed some methods to overcome these difficulties. This paper provides further details about their implementation in such a real industrial background. First of all, we observed that the usual Holt-Winters monitoring procedure could be severely biased by the presence of outlier observations. Consequently, we implemented a robust approach, as proposed by Gelper et al. [START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] and Croux et al. [START_REF] Croux | Robust control chart for time series data[END_REF] Nonetheless, we noticed that the efficiency of this robust Holt-Winters control is seriously affected when structural changes occur. That is why we developed a dynamical procedure, which detects structural breaks in a time series and reinitializes the robust smoothing parameters accordingly. This dynamical procedure is an original contribution to References [START_REF] Croux | Robust control chart for time series data[END_REF][START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] , contributing to improve the flexibility for applications in industrial fields. Finally, it should be mentioned that an industrial monitoring implies a simultaneous control of several dozens of variables. There are numerous approaches in multivariate process control (see e.g. Bersimis et al. [START_REF] Bersimis | Multivariate statistical process control charts: An overview[END_REF] for a review), and a specific investigation of robust techniques in this context is beyond the scope of the paper. Alternatively, we show how the robust controlling procedure can be adapted to the framework of multivariate process control.

This paper is organized as follows. In Section 2, we introduce the main statistical methods that have been deployed to control the STMicroelectronics IT system: The robust Holt-Winters monitoring proposed by Gelper et al. [START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] and Croux et al.; [START_REF] Croux | Robust control chart for time series data[END_REF] The dynamical contribution developed through our researches, improving its flexibility in changing environments. In Section 3, some performance tests are introduced. They are grounded on quantitative simulations to: 1) Compare the usual and robust Holt-Winters monitoring procedures; 2) Evaluate the structural break detection capacity of the new dynamical procedure. Lastly, Section 4 introduces some example from real industrial cases: Univariate and multivariate applications are considered.

Investigation of a robust monitoring for trended time series with structural changes

We recall here the main principles of Holt-Winters based monitoring, and its robust version introduced by References 4, 6 . Then, we propose a new methodology in presence of structural changes.

Monitoring based on Holt-Winters smoothing

The Holt-Winters (HW) algorithm is a popular technique used to provide short-term forecasts of a given time-series (see e.g. Makridakis et al. [START_REF] Makridakis | Forecasting, methods and applications[END_REF] ). The predictions are built iteratively as a linear combination of the observed values and the prediction obtained at last step. For illustration, let us consider a time series y, observed at dates 1, 2, . . . , n -1. For the sake of simplicity, we assume that y is non-seasonal, though the methodology is similar. The HW algorithm is based on the assumption that y is a sum of two time-series α and β corresponding respectively to a local level (order of magnitude) and a trend. These auxiliary time series are estimated iteratively as averages of the last observation and the last predictions, weighted by two parameters λ 1 and λ 2 :

αt = λ 1 y t + (1 -λ 1 )ŷ t|t-1 , t = 1, . . . , n -1 (1) βt = λ 2 (α t -αt-1 ) + (1 -λ 2 ) βt-1 , t = 1, . . . , n -1 (2)
Logically, the one-step-ahead forecast done at date t -1 for date t is then given by:

ŷt|t-1 = αt-1 + βt-1 , t = 1, . . . , n (3)
which gives in particular the prediction at date n. In practice, λ 1 and λ 2 are estimated by minimizing a criterion (often the least-square criterion) based on the forecast errors:

E t = y t -ŷt|t-1 , t = 1, . . . , n -1 (4)
and the algorithm is initialized by a linear regression on the first m values.

As a second step, monitoring can be performed. While applying a control chart to y is not recommended, due to the violations of the usual assumption (identically and independently distributed data) especially if y is trended, the forecast errors E t may be close to satisfy it. Then, assuming furthermore that E t are normally distributed N (0, S 2 ) the upper and lower control limits for E t are then given by:

U CL = +q α/2 * Ŝ (5) LCL = -q α/2 * Ŝ (6)
where q α/2 is the quantile of a Student distribution at level α/2, and Ŝ2 is the usual variance estimator:

Ŝ2 = 1 n -m -1 n-1 t=m+1 E 2 t (7)
These limits intend to detect the dates that correspond to an anomaly: A value of E t outside the interval [UCL, LCL] should be a strong indication of an abnormal behavior (for a given confidence level α). However, the limits themselves are sensitive to outliers, since the variance estimator overestimates the true variance in presence of outliers. Furthermore, the predicted value ŷt|t-1 depends linearly on past values that may contain outliers. These problems are solved by the robust version of the Holt-Winters monitoring.

Robust monitoring based on Holt-Winters smoothing

The Robust Holt-Winters algorithm (RHW) introduced by Gelper et al. [START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] and Croux et al. [START_REF] Croux | Robust control chart for time series data[END_REF] considers two additional auxiliary time series: y * , representing a cleaned proxy of y after outliers treatment, and σ the expected prediction error, representing a robust estimate of the forecast error E t . To obtain a robust algorithm, large values are truncated when larger than a given threshold. More precisely, the expected errors σ t are computed recursively by:

σ2 t = λ σ ψ k E t σt-1 2 σ2 t-1 + (1 -λ σ )σ 2 t-1 (8)
where λ σ is a given weight, and ψ k is the Huber function with boundary value k:

ψ k (x) = x if |x| ≤ k, sign(x) × k if |x| > k (9)
The error E t is still given by E t = y t -ŷt|t-1 , with ŷt|t-1 = αt-1 + βt-1 , but the local level and trend are now estimated by using the cleaned time series y * :

αt = λ 1 y * t + (1 -λ 1 )ŷ t|t-1 (10) βt = λ 2 (α t -αt-1 ) + (1 -λ 2 ) βt-1 (11)
where y * t is given by:

y * t = ψ k E t σt × σt + ŷt|t-1 (12)
Notice that the role of the Huber function ψ is to truncate the forecast errors E t when larger than k times the expected prediction error σt :

ψ k E t σt × σt = E t if |E t | < kσ t , sign(E t ) × kσ t if |E t | > kσ t (13)
For instance, Equation ( 12) can be rewritten in a simpler way:

y * t = y t if |E t | < kσ t , sign(E t ) × kσ t + ŷt|t-1 if |E t | > kσ t (14)
Finally the parameters λ 1 and λ 2 , and the standard deviation S of the forecast errors are computed with robust procedures, while λ σ is fixed to 0.3 (see References [START_REF] Croux | Robust control chart for time series data[END_REF][START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] for more details). The algorithm initialization is also done in a robust way by repeated median regression. For that, we recommend using a period of length 7: This is short enough to assume a local linear trend and long enough to be resistant to 2 outliers.

Since the RHW smoothing is a full robust procedure, the control charts based on the errors E t (see Section 2.1) are now resistant to outliers, which is a clear improvement to the (non-robust) HW-based monitoring. However, an adaptation is necessary in presence of a structural change.

Robust monitoring with structural changes

In this section, we focus on time series that possibly contain outliers and structural changes. To detect structural changes in parametric models, three main classes exist (Zeileis 22 ): F statistics, fluctuation tests and maximum likelihood scores. We choose to use a common and simple F statistics, the Chow test, since it can be easily adapted to a robust framework. It consists in splitting the sample into 2 groups: The first one before the break date and the second one after. The model parameters are estimated for both of them so that an F test be performed to judge whether they are equal or not. The Chow test is easy to use but restricted by 2 limitations. The first one is mentioned by Hansen 7 : An important limitation of the Chow test is that the break date must be known a priori ". Moreover, the exact number of changes is unknown. The second problem is a question of robustness: A break may be missed or falsely detected because of outliers. In this section, we show how the robust Holt-Winters smoothing allows to deal with these problems and deduce a strategy for structural change monitoring.

An introducing example Consider the time series in Figure 1 with an outlier at date 10 and a break at date 26. As expected, the robust algorithm is not sensitive to the outlier at date 10 contrarily to the classical smoothing: This is its main advantage. But after that, the level changes suddenly at date 26. The robust algorithm does not admit this modification quickly and many false alarms follow. Nevertheless, this specificity can be used to detect break dates. Indeed, when a structural change happens, there is a quite long period (here 26-34) of successive false alarms corresponding to successive large errors in the RHW smoothing. During this period, the predicted values ŷt|t-1 seem to exhibit a deterministic pattern (Figure 1,d): Actually, we show below (Figure 2 and Proposition 1) that they match exactly with an exponential function. These two facts strongly suggest that the periods of successive large errors given by the robust Holt-Winters smoothing are useful for break dates detection. forewarning sign of structural change. More formally, we call t 1 a suspicious date for structural change if there exists an integer p ≥ 3 such that:

ψ k E t 1 σt 1 = ψ k E t 1 +1 σt 1 +1 = ... = ψ k E t 1 +p-1 σt 1 +p-1 = k (15)
or

ψ k E t 1 σt 1 = ψ k E t 1 +1 σt 1 +1 = ... = ψ k E t 1 +p-1 σt 1 +p-1 = -k (16)
The period [t 1 , . . . t 1 + p -1] is a suspicious period : a period when forecasting errors remain k times higher than their expected values.

Proposition 1 During a suspicious period [t 1 , . . . t 1 + p -1], the predictions ŷt|t-1 of the robust HW smoothing are given by a deterministic and monotonic function f which does not depend on any observation posterior to date t 1 . Its form is:

f(t) = r (t+1) (r -1) 2 + a.t + b ( 17
) with r = λ σ (k 2 -1) + 1 = 1.378405 (18)
a and b are given by ŷt 1 and ŷt 1 -1

Proof See Appendix 1.

A new methodology for structural change monitoring

The previous observations and result about the RHW smoothing in case of a structural change suggest the following methodology that we call RHW-SC:

1. Find the suspicious dates by looking at consecutive relative errors given by the robust HW smoothing (Equations 15 and 16)

2. Apply a robust version of the Chow test to the suspicious dates detected in 1. One robust version of the Chow test consists in replacing the usual linear regression by the repeated median regression (Siegel [START_REF] Siegel | Robust regression using repeated medians[END_REF] , Rousseeuw et al. 17

) and using a robust estimator of the residuals' variance (Croux et al. [START_REF] Croux | Robust control chart for time series data[END_REF] ) to compute the F statistic.

This strategy tackles the two main issues mentioned at the beginning of this section: All the possible break dates (number and locations) are automatically detected by the algorithm itself, and for a given date the statistical test for structural change is done in a robust way. In practice the methodology is applied dynamically, and when the robust Chow test is positive, the robust HW smoothing is reinitialized.

Performance of the RHW-SC methodology on the introducing example

To solve the problem raised by the introducing example, we perform dynamically the RHW-SC methodology with p=3. The methodology has detected only one suspicious date, namely day 26, which indeed corresponds to the structural change. Notice that only one Chow test has been used contrarily to exhaustive methods such as Quandt-Anderson that systematically test all dates. We reinitialized the robust HW smoothing at date 26, by using the repeated median regression coefficients. As a result, RHW-SC methodology remains resistant to the outlier at date 10 without generating a long sequence of false alarms after the break.

Global performance tests

Some ARIMA models correspond to exponential smoothing methods (see e.g. Hyndman et al. [START_REF] Hyndman | Local Linear Forecasts Using Cubic Smoothing Splines[END_REF] . In particular, a Holt-Winters forecasting with smoothing parameters (λ 1 , λ 2 ) is equivalent to an ARIMA (0, 2, 2) model of parameters (θ 1 , θ 2 ) if: -1 ≤ θ 2 ≤ 1 , θ 2 -θ 1 ≤ 1 and θ 2 + θ 1 ≤ 1, with :

λ 1 = 1 + θ 1 ( 19
)
λ 2 = 1 -θ 1 -θ 2 1 + θ 2 (20)
Thus, we first use ARIMA (0, 1, 1) time series with parameter θ = 0.5. Indeed, they correspond to the special case where θ 2 = 0. Next, we perform our simulations on ARIMA (0, 2, 2) models with (θ 1 , θ 2 ) = (1 , -0.25). These parameters have been estimated from an industrial time series by assuming that it comes from an ARIMA(0, 2, 2) model. In this section, we compare first the performances of the RHW and the HW smoothing statically and dynamically. Afterwards, we evaluate the ability of the RHW-SC method to detect a structural change.

Comparison tests in a static setting

We are interested in comparing 2 characteristics of the HW and the RHW methods:

1. The power : The probability that an outlying observation is detected 2. The false detection rate: The probability that a normal observation is detected.

This risk is called size of the control chart when there is no outlier. Below, we use 0.95 as a confidence level, so the false detection rate is expected to be 0.05. The RHW and HW monitoring are equivalent when there is no outlier (Reference 4 ). Their differences become significant when the data is contaminated.

On the one hand, Table 1 and Table 2 below show that the HW control chart is subject to two effects. The first one is the widening of its control limits by outliers, which tends to reduce abnormaly its false alarms rate. The second one concerns the dates following these outliers: the corresponding predictions are biased and this fact tends to raise the false alarms rate. Here, the first effect is predominant. This explains the too low false alarms risk and the poor detection rate of the HW method.

On the other hand, Table 3 and Table 4 show that the RHW methodology outperforms the HW smoothing especially when the number of outliers raises and when the magnitude of these outliers is not fixed: which is a realistic case. In fact,the RHW control chart doesn't lose neither its power, nor the stability of its false alarms rate that remains stable around the theoretical value 5%. The same strategy as in the precedent section is used. Only the first step is modified to generate ARIMA (0, 2, 2) models with (θ 1 , θ 2 ) = (1 , -0.25). The results are summarized below (Tables 5,6, 7 and 8). They show that the RHW control chart remains better in term of detection. Nevertheless, its false risk becomes higher than expected even if it remains stable enough for a fixed contamination rate contrarily to the HW smoothing. In practice, industrial variables are tracked daily; the smoothing parameters and control limits are re-estimated every day. So, we perform dynamical simulations that consist in using the same procedure as for the static setting but with updating the smoothing parameters and controls limits at each iteration. 11 and12 show how suitable the HW monitoring is to update its smoothing parameters and control limits contrarily to the RHW method. Indeed, the HW detection rate is drastically improved even if the power of the RHW smoothing remains the better. Let's remark that the HW false alarm risk has raised and approaches better the theoretical value 5%. As for the RHW false alarms risk, it has decreased to approach this same theoretical value; but this later improvement is less obvious. However, the RHW false alarms risk has become very stable and non dependant on the magnitude of outliers, which is not the case for the HW control chart. Anyway, the RHW control chart outperforms the HW one when the scale of the outliers is more realistically (randomly) chosen.

Structural change detection performance tests

To evaluate the RHW-SC methodology, we use the following strategy: Table 13: Detection rate

B = 0 B = ±50 A = 0 - 95 A = ±50 86 95
Table 14: False detection rate

B = 0 B = ±50 A = 0 - 6 A = ±50 4 5
The RHW-SC methodology is efficient enough for structural change detection especially when a shift happens. It looses in efficiency when there is no shift. But in these latest cases, the RHW period of successive false alarms is short. The RHW-SC false detection rate corresponds to the theoretical value expected since the confidence level used for the Chow test is 95%.

Applications

The RHW and the RHW-SC methods are performed daily on many indicators of the information system of a company. In general, the results are satisfactory. Here, we present some of these examples and a multivariate case.

Examples of univariate time series

This first example is illustrative of structural changes detection by the RHW-SC methodology. 5 points out the sensitivity of the HW method to outliers. The abnormal observation of the third week has deformed the seasonal component. Because of the overestimation of this seasonality, the 26th observation that is really aberrant is not detected. Moreover, it leads to four false alarms the two following weeks. These problems no longer exist with the seasonal RHW smoothing.

The multivariate case

In our company, several indicators are tracked daily. So, monitoring them separately leads to false alarms every day. This is foreseeable given that the theoretical false alarms risk tends to 100% when the number of independent variables approaches 100: Whence the necessity to perform a multivariate monitoring. Notice that there are a lot of possible approaches in the multivariate framework (see e.g. Bersimis et al. [START_REF] Bersimis | Multivariate statistical process control charts: An overview[END_REF] ), and our aim is neither to do a comparative study nor to propose a best one. Rather, we show how the new robust methodology can be adapted to the multivariate case. Thus, as an example, we constitute groups of variables. For each group of p variables we use the following strategy :

1. Perform the RHW smoothing for each variable of the group.

2. Use a robust Hotelling T2 control chart to analyze simultaneously the p vectors of residuals given by step 1.

Among existing robust Hotelling T2 control charts (see e.g. Rousseeuw, [START_REF] Rousseeuw | Least Median of Squares Regression[END_REF] Alfaro et al. [START_REF] Alfaro | A robust alternative to Hotelling's T2 control chart using trimmed estimators[END_REF] ) we consider here the computation of the confidence ellipsoid with the Minimum Covariance Determinant criterion introduced by Rousseuw et al. [START_REF] Rousseeuw | Least Median of Squares Regression[END_REF][START_REF] Rousseeuw | A Fast Algorithm for the Minimum Covariance Determinant Estimator[END_REF] . There are two identified difficulties: A poor orientation of the confidence ellipsoid and an underestimation of its size. To face these problems, we use recent solutions found in the literature. Firstly, the orientation is improved by choosing a subset of the 75 % best points for the MCD criterion instead of 50% (Huber et al. [START_REF] Hubert | Fast and robust discriminant analysis[END_REF] ). Secondly, the size estimation is improved by using two correcting factors: One asymptotic factor of consistency to the chi square distribution (Rousseeuw et al. [START_REF] Rousseeuw | Robust Regression and Outlier Detection[END_REF] ) and one empirical result for small samples (Pison et al. [START_REF] Pison | Small Sample Corrections for LTS and MCD[END_REF] ). Now, we present an example with a group of p = 4 variables. Each day, the chi square statistic is computed for this group. The result is plotted in Figure 6. This control chart has successfully detected the 6 most important outliers but does not say where the problems come from. This is the well known problem of multidimensional out of controls interpretation. As an example, let us focus on the alarm on date 56. Its cause can be known by looking at each variable separately. Then, Figure 7 shows that the variable BIN2 is mainly responsible for it. Nevertheless, this solution is not realistic for an industrial use because it produces too many graphs that need to be examined. Among the numerous existing criteria to interpret multivariate signals, the partial relative contributions mentioned by Montgomery [START_REF] Montgomery | Statistical Quality Control: A Modern Introduction[END_REF] are very popular for their efficiency. In Figure 7, these relative contributions at date 56 confirm the influence of BIN2.

Conclusion

In this paper, we have studied a monitoring based on robust Holt-Winters smoothing, as proposed by Gelper et al. [START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] and Croux et al. [START_REF] Croux | Robust control chart for time series data[END_REF] Based on an industrial application of this method at STMicroelectronics and simulation studies, we can confirm its high robustness. However, as poor results were noticed in case of structural changes, we have also proposed and evaluated an improved dynamical approach, for a better integration of changing environments. The efficiency of this robust and dynamical method has been demonstrated on real univariate and multivariate STMicroelectronics case studies. This contribution has been developed as an improvement for the STMicroelectronics IT system monitoring, but it could be applied to any other industrial applications, where time dependent variables has to be statistically controlled.

For further research, there are several interesting outstanding questions. Firstly, some tuning choices could be investigated. Hence, the truncation parameter was set to k = 2 and the smoothing parameter λ σ to 0.3 (still proposed by References [START_REF] Croux | Robust control chart for time series data[END_REF][START_REF] Gelper | Robust forecasting with exponential and Holt-Winters smoothing[END_REF] ). These values could be discussed: Do they always provide optimal results or should they be contextually adapted? Secondly, when structural changes occur, we decided to reinitialize the Holt-Winters parameters. An alternative could be to use dynamical smoothing parameters that may change over time, as proposed by Williams 21 or Taylor [START_REF] Taylor | Smooth transition exponential smoothing[END_REF] . Finally, in our methodology, we chose a Chow test to detect structural breakpoints, since it is easily adapted to the robust framework. Nonetheless, as mentioned earlier, other approaches exist (other F statistics, fluctuation tests, maximum likelihood scores). Their potential of application in our monitoring procedure and a comparative study or their performance may bring to interesting further insights.

  (a) Alarms given by the classical HW forecasting (b) Predicted values by the classical HW smoothing (c) Alarms given by the robust HW forecasting (d) Predicted values by the robust HW smoothing

Figure 1 :Figure 2 :

 12 Figure 1: A real time series with outlier and structural change

  (a) Alarms given by the RHW-SC methodology (b) Predicted values with the RHW-SC methodology

Figure 3 :

 3 Figure 3: Alarms and predicted values by the robust version with structural change detection

3. 1 . 1 1 .

 111 Simulation results for ARIMA (0, 1, 1) time series We use the following strategy: Generation of a time series: simulate an ARIMA (0, 1, 1) time series of length 160 with parameters θ = 0.5. The first 60 values serve as a training sample (including m = 7 values for HW and RHW initializing). 2. Generation of outliers: For a fixed contamination rate R, choose randomly 160R dates among the 160 dates. Contaminate these observations by adding or subtracting to them a value e. Whether it is an addition or a subtraction is chosen at random. Two cases are considered for e: (a) e is a fixed value among 10, 20 or 30 (b) "Mix": e is chosen at random uniformly among 10, 20, 30 3. Perform the HW and the RHW methods and estimate: (a) Their power : as the percentage of the contaminated observations that are really detected (b) Their false detection rate: as the percentage of non contaminated observations that are detected. 4. Repeat 100 times the steps 1 to 3.

1 . 3 . 4 .

 134 Generation of a time series: simulate an ARIMA (0, 2, 2) time series of length 160 with parameters (θ 1 , θ 2 ) = (1,-0.25). The first 60 values serve as a training sample (including m = 7 values for RHW initializing). 2. Generation of outliers: For a fixed contamination rate R = 0.05, choose randomly 160R dates among the 160 dates. Contaminate these observations by adding or subtracting to them a fixed value e = 5. Whether it is an addition or a subtraction is chosen at random. Generation of structural changes: choose randomly 1 date among the 100 last values. From this date to the last one, add a linear function: A * DAT E + B to the time series. Perform the RHW-SC methods with 0.95 as confidence level for the robust Chow test and compute: (a) The detection rate for structural changes (b) The false detection rate for structural changes. 5. Repeat 1000 times the steps 1 to 4.

Figure 4 :

 4 Figure 4: An industrial time series with a structural change

Figure 5 :

 5 Figure 5: An industrial and seasonal time series

Figure

  Figure5points out the sensitivity of the HW method to outliers. The abnormal observation of the third week has deformed the seasonal component. Because of the overestimation of this seasonality, the 26th observation that is really aberrant is not detected. Moreover, it leads to four false alarms the two following weeks. These problems no longer exist with the seasonal RHW smoothing.

Figure 6 :

 6 Figure 6: Chi squared daily statistics for a group of 4 variables

Figure 7 :

 7 Figure 7: RHW smoothing for the four variables of the group

Figure 8 :

 8 Figure 8: Partial contributions to variability at date 56

Table 1 :

 1

		HW detection rate	Table 2: HW false alarms rate
		R = 2% R = 5% R = 10%		R = 2% R = 5% R = 10%
	e = 10	93.5	70.0	31.0	e = 10	3.3	3.6	2.1
	e = 20	94	74	35	e = 20	1.7	3.0	7.5
	e = 30	91	72	36	e = 30	1.9	6.0	10
	Mix	72	55	33	Mix	3	5.3	4

Table 3

 3

		: RHW detection rate	Table 4: RHW false alarms rate
		R = 2% R = 5% R = 10%		R = 2% R = 5% R = 10%
	e = 10	100	100	99.2	e = 10	5.2	4.8	6.1
	e = 20	100	100	99.9	e = 20	5.2	5.2	6.4
	e = 30	100	100	100	e = 30	5.2	4.8	7.1
	Mix	100	100	99.8	Mix	5.1	4.9	6.1
	3.1.2 Simulation results for ARIMA (0, 2, 2) time series

Table 5

 5

		: HW detection rate	Table 6: HW false alarms rate
		R = 2% R = 5% R = 10%		R = 2% R = 5% R = 10%
	e = 10	97.6	91.3	67.0	e = 10	3.3	2.7	3.7
	e = 20	100	92.2	74.4	e = 20	2.1	4.1	8.0
	e = 30	97.1	91.0	73.3	e = 30	5.2	1.5	4.0
	Mix	92.7	77.3	57.5	Mix	2.2	3.9	6.4

Table 9

 9

		: HW detection rate	Table 10: HW false alarms rate
		R = 2% R = 5% R = 10%		R = 2% R = 5% R = 10%
	e = 10	99.6	95.4	80	e = 10	4.6	4.3	4.2
	e = 20	100	93.6	80	e = 20	3.3	2.2	2
	e = 30	100	97	78.6	e = 30	2.5	2.2	2.4
	Mix	95	70	52	Mix	3.0	3.2	2.0
	Table 11: RHW detection rate	Table 12: RHW false alarms rate
		R = 2% R = 5% R = 10%		R = 2% R = 5% R = 10%
	e = 10	100	100	99.7	e = 10	6.3	6.8	8.6
	e = 20	100	100	100	e = 20	6.0	7.1	9.1
	e = 30	100	100	100	e = 30	6.2	6.9	9.3
	Mix	100	100	100	Mix	6.1	7.0	7.2
	Tables

Acknowledgements

The authors would like to thank the company STMicroelectronics for giving the opportunity to work on this research project, and especially X. Ambrosioni and P. Palazon, for their support. We also acknowledge the participants of the ENBIS-12 conference (Ljubljana, 9-13 September 2012), for their relevant remarks on our first oral presentation of this work.

Appendix

We prove the Proposition presented in Section 2.3.1: If t 1 , t 1 + 2, . . . , t 1 + p -1 is a suspicious period, then the sequence ŷt 1 , ŷt 1 +1 , . . . ŷt 1 +p-1 is given by a deterministic and monotonic function f which does not depend on any observation posterior to date t 1 .

Proof

The period t 1 , t 1 + 1, . . . t 1 + p -1 is defined by either Equation 15 or Equation 16. Without lost of generality, consider Equation 15. By combining this definition 15 and Equation 8, we obtain the equation:

which shows that the predicted errors σt follow a geometric progression with common ratio,

Consequently, by denoting t 0 := t 1 -1, the date before t 1 , the predicted errors are given by: σt = r t-t 0 σt 0

With the common choice k = 2, λ σ = 0.3, we have r ≈ 1.378 > 1. Thus, the expected error goes increasing exponentially. Furthermore, the cleaned time series y * t becomes : y * t = k.r t-t 0 σt 0 + ŷt (24)

Now, relying on the equivalence with ARIMA(0,2,2) model (see e.g. Hyndman 10 ), the forecast values of the Holt-Winters smoothing with parameters (λ 1 , λ 2 ) follow the recursive scheme:

where θ 1 and θ 2 are the parameters of the corresponding ARIMA(0, 2, 2) model (See Equations 19 and 20). Given Equation 24, this scheme becomes:

with:

This is a linear equation, whose solutions are given by the sum of the solutions of the homogeneous linear equation and a particular solution:

The values of the constants a and b are imposed by ŷt 1 and ŷt 1 +1 . This shows indeed that the predicted values are purely deterministic and increasing exponentially during the suspicious period.


Dominique Guégan
email: dguegan@univ-paris1.fr

Bertrand K Hassani
email: bertrand.hassani@malix.univ-paris1.fr

Cédric Naud
email: cedric.naud@aon.fr

An efficient threshold choice for operational risk capital computation

Keywords:

Operational risk quantification requires dealing with data sets which often present extreme values which have a tremendous impact on capital computations (VaR). In order to take into account these effects we use extreme value distributions, and propose a two pattern model to characterize loss distribution functions associated to operational risks: a lognormal on the corpus of the severity distribution and a Generalized Pareto Distribution on the right tail. The threshold from which the model switches form a scheme to the other one is obtained using a bootstrap method. We use an extension of the Peak-over-threshold method to fit the GPD and the EM algorithm to estimate the lognormal distribution parameters. Through the VaR, we show the impact of the GPD estimation procedure on the capital requirements.

Besides, our work points out the importance of the building's choice of the information set by practitioners to compute capital requirements and we exhibit some incoherences with the actual rules. Particularly, we highlight a problem arising from the granularity which has recently been mentioned by the Basel Committee for Banking Supervision.

Introduction

The purpose of this paper is operational risks quantification following the Basel Advanced Measurement Approach (AMA) [START_REF] Bis | International convergence of capital measurement and capital standards[END_REF]). This proposal of the Basel committee binds banks to carry out their own models on internal data sets to evaluate amounts of capital necessary to face these risks.

In this paper, we focus on the computation of the loss distribution function (LDF) [START_REF] Frachot | Loss distribution approach for operational risk[END_REF], [START_REF] Guégan | A modified panjer algorithm for operational risk capital computation[END_REF]) which is currently used to evaluate this kind of risks [START_REF] Cruz | Operational Risk Modelling and Analysis[END_REF], [START_REF] Chernobai | Operational Risk: A Guide to Basel II Capital Requirements, Models, and Analysis[END_REF] and [START_REF] Shevchenko | Modelling Operational Risk Using Bayesian Inference[END_REF]). We compute it as a convolution of two distribution functions, one modeling the frequencies of the losses and the other one their severity. Following the literature, we restrict to a Poisson distribution for the frequencies, and in this paper we focus on the estimation of the severity distribution.

The objective of the paper is to model the severities using a class of distributions which take into account the large losses: the Generalized Pareto Distribution (GPD). In order to capture extreme values on the first hand, and on the other hand the multimodality observed in operational risk data sets1 , we propose to model the severities as a mix of two distributions: one characterizing the most important losses denoted cauda using the GPD [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF], [START_REF] Davis | Tail estimates motivated by extreme value theory[END_REF] and [START_REF] Embrechts | Modelling Extremal Events: for Insurance and Finance[END_REF]), and another one on the remaining data corresponding to the central part of the distribution or corpus. A characteristic of the GPD stands in the determination of a threshold denoted u permitting to define the data set on which the distribution will be adjusted. Therefore, the estimation of this threshold is the core of this paper: we highlight the influence of the choice of the threshold on the computation of the capital requirement. Then, once the threshold has been identified we need to estimate the parameters of the two distributions, for instance the lognormal and the GPD. Using several estimation methods, we also exhibit the influence of these ones on the computation of univariate capital requirement.

Finally, considering the obtained results, we discuss the problem of granularity which is the main problem in the construction of operational risk data sets (or Basel matrix).

The paper is organized as follows: Section two introduces the data set, Section three presents the theoretical framework, Section four describes the estimation of the threshold and the construction of the LDF and Section five is devoted to capital charge estimations. A sixth Section concludes.

The Data Set

The data set we use is organized into the Basel Matrix (BIS ( 2001)). In its first level of granularity, this matrix is made up of 56 cases -8 business lines ("b") × 7 event types ("e")2 . Nevertheless, each event type might be decomposed in several elements. For example, the "external fraud" event may be shared in two items -"Theft and Fraud" and "Systems Security" (second level of granularity). In a third level, the element "Theft and Fraud" may be split in several components: "Theft/Robbery", "Forgery" and "Check kiting". After a deep analysis, we observe that the kind of losses expected from a fraud with a credit card does not correspond to losses caused by someone hacking the system for instance; nevertheless they are in the same cell. Therefore, considering the largest level of granularity, we could face multimodal empirical distributions.

Consequently, the methods used to model the losses depend on the granularity level choice. This choice might have a tremendous impact on capital requirement computations. Besides, we face a trade-off between quantity of data and robustness of the estimations: indeed, if the quantity of data is not sufficient, we cannot go lower in the granularity; on the another hand the ensuing empirical distribution is therefore an aggregate of various natures of data and the estimation of this last empirical distribution can be source of unusable or unrealistic results.

The previous paragraph dealt with the multimodality and the granularity problems. These points have a corollary: an extreme value is a relative concept, which require a referent. For example, a card fraud of 50000e is an extreme value compared to an average loss of 300e while a rogue trading of 50000e is not an extreme value when we know that the average loss represents several millions. Furthermore, we do not define an extreme value as an "extremely" rare event.

Therefore, one can use a GPD approach as the one presented below to model the severities.

The theoretical framework

The Loss Distribution Function (LDF) G (whose density is denoted g) is a mixture of two distributions, the frequency distribution p and the severity distribution F (whose density is f ):

G = ∞ γ=0 p(γ)F * γ (x; θ), x > 0, G = 0, x = 0, (3.1)
where * γ is the γ-order operator of convolution, θ a vector of parameters. The density is given by,

g = ∞ γ=0 p(γ)f ⊗γ (x; θ), x > 0. (3.2)
The severity distribution, on which we focus in the current paper, is defined as a mixture of a lognormal distribution on the corpus, and a GPD on the right tail (cauda). The GPD density is:

f GP D (x; u, β, ξ)) =      1 β 1 + ξ x-u β -1-( 1 ξ ) , if x ≥ u, 1 + ξ (x-u) β > 0, β > 0 1 β 1 -x-u β , if x ≥ u, ξ = 0 , ( 3.3)
where u ∈ R * + is the threshold, β ∈ R * is the scale parameter and ξ ∈ R the shape parameter.

Therefore, the severity density in our model is as follows:

f (x; u, β, ξ, µ, σ) =      f (corpus) (x; µ, σ), if x < u f (cauda) = 1 1- u 0 f (corpus) (x;µ,σ)dx × f GP D (x; u, β, ξ)), if x ≥ u , (3.4)
where, µ and σ are the lognormal distribution parameters, and f GP D is the density introduced in (3.3).

Remark 3.1. Modeling both the losses below and above the threshold has two major advantages:

• It enables specifying the threshold regarding a model, i.e. the quantile obtained empirically is different from the quantile we have theoretically.

• It permits creating a whole distribution we will be able to use in a multivariate approach [START_REF] Guégan | n-dimensional copula contributions to multivariate operational risk capital computations[END_REF]).

In order to estimate the parameters of the distribution (3.4), we need first to estimate the threshold u. To do so, we carry out a bootstrap method [START_REF] Hall | Using the bootstrap to estimate mean squared error and select smoothing parameter in nonparametric problems[END_REF], [START_REF] Danielsson | Using a bootstrap method to choose the sample fraction in tail index estimation[END_REF] for which we give practical solutions in the next Section. Once u has been found, ξ and β have to be estimated, therefore we implement a method based on the Anderson-Darling statistics [START_REF] Luceno | Fitting the generalized pareto distribution to data using maximum goodnessof-fit estimators[END_REF]). This method maximizes the following relationship to estimate the parameters:

A 2 n = n ∞ -∞ (F (GP D) (x; u, β, ξ) -S n (x)) 2 F (x; u, β, ξ)(1 -F (GP D) (x; u, β, ξ))) dF (GP D) (x; u, β, ξ), (3.5)
where, F (GP D) (x; u, β, ξ) is the cumulative distribution function of a GPD whose density is

given by (3.3), and n ∈ R + the number of observations. In the fifth Section, we compare this method with other traditional estimators and we show that this approach provides interesting results. Then, as we assumed a lognormal distribution to model the central part of the severity distribution (3.4), we implement an Expectation-Maximization algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF], [START_REF] Mclachlan | The EM algorithm and extensions[END_REF]) to estimate the parameters of the truncated distribution. Finally, the λ parameter of the Poisson distribution p(γ; λ) is estimated using maximum likelihood method. As soon as we have estimated the whole set of parameters, we can build the loss distribution function (3.1) using a convolution method based on the modified Monte Carlo Algorithm given in Appendix B [START_REF] Fishman | Monte Carlo[END_REF], [START_REF] Cruz | Modeling, measuring and hedging operational risk[END_REF], [START_REF] Peters | Simulation of the annual loss distribution in operational risk via panjer recursions and volterra integral equations for value at risk and expected shortfall estimation[END_REF] and Shevchenko ( 2010)), and compute a VaR [START_REF] Riskmetrics | Var[END_REF], [START_REF] Artzner | Coherent measures of risk[END_REF]),

V aR (1-α)% = inf(x ∈ R : P (X > x) ≤ (1 -α)), (3.6)
where, X is a random variable characterized by the previous LDF, and α = 0.1.

The Figure 1 illustrates the method introduced in this section through a simulation experiment.

We Remark 3.2. Drawing a parallel between this two-behavior approach with other methods such as using g-and-h distributions [START_REF] Peters | Bayesian inference, monte carlo sampling and operational risk[END_REF], [START_REF] Dutta | A tale of tails: an empirical analysis of loss distribution models for estimating operational risk capital[END_REF] and [START_REF] Degan | The quantitative modeling of operational risk: between g-and-h and evt[END_REF], we conclude that it does not exist a method which supplant the others. Trying the gand-h distribution to model some of our data sets, we had most of the time bad fitting results or unrealistic capital charges, nevertheless in some particular cases the flexibility offered by this distribution could be a reliable solution to model multimodal severity distributions. The method presented in this paper, especially when we have enough data above the threshold and as soon as the goodness-of-fit test validates the distribution, may be a powerful alternative.

4 Discussion around the threshold "u" of the GPD

In this section, we focus on the way to determine an efficient threshold u consistent with the specific data set we work with, that is to say operational risks losses. In order to solve this problem, we assume that we observe a data set X = (X 1 , ..., X n ), and we denote X = (X (1) , ..., X (n) )

its order statistics, We first estimate the parameter ξ of the GPD (3.3) with the whole sample using the traditional Hill estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]):

ξn (k) = 1 k k i=1 log X (n-i+1) -log X (n-k) , (4.1)
k being the number of losses above the threshold. Thus, by the way we introduce a direct relationship between k and u. Indeed, once we have k we have the threshold and vice versa.

We first determine k using the Hill estimate (4.1) and then implement a bootstrap method to obtain u. The different steps are the following:

1. We plot ξn (k) with respect to k, as in Figure 2.

2. As soon as we detect a change in the plot pattern ("break"3 ) as in Figure 2 we build an information subset X 1 (size n 1 ) containing this break and the values above it. For example, in Figure 2 we heuristically determine a break at the value 100 e which corresponds to the initial value of the subsample. Using the values in X 1 we determine the parameter k 1 minimizing the following criterion: 3. In the relation (4.2), the parameter ξ is unknown and we estimate it bootstrapping ξn 1 (k 1 ) as follows. From the set previous X 1 , we draw J (J ∈ N) subsamples with replacement of size n 2 (n 2 < n 1 ), (X

AM SE(k 1 ) := E( ξn 1 (k 1 ) -ξ)
(j) 1 , ..., X (j)
n 2 ), j = 1, ..., J. Their corresponding order statistics are (X

(1) , ..., X

(n 2 ) ) and the corresponding bootstrapped Hill estimate of ξn 1 (k 1 ) from X 1 is equal to,

ξn 2 (k 2 ) = 1 k 2 k 2 i=1 log X (n 2 -i+1) -log X (n 2 -k 2 ) , k 2 = 1, ..., n 2 , (4.3)
where

X l = 1 J J j=1 X (j) (l) , l = 1, ..., n 2 . (4.4)
4. Now, replacing ξ by ξn 2 (k 2 ) in (4.2), we minimize the following criterion with respect to

k 2 . AM SE(k 2 ) := E( ξn 1 (k 1 ) -ξn 2 (k 2 )) 2 .
(4.5) 5. Then, given (k 1 , k 2 ), and following [START_REF] Hall | On some simple estimates of an exponent of regular variation[END_REF][START_REF] Hall | Using the bootstrap to estimate mean squared error and select smoothing parameter in nonparametric problems[END_REF] we estimate4 the threshold u by,

u = k 2 * n 1 n 2 2 3 . (4.6)
In order to implement the previous methodology, we need to pay attention to the following points:

1. If an obvious pattern change (from erratic to linear for example) is observed in the representation of ξn (k) with respect to k (Step 1), we can directly read u on the graph.

2. In the step 4, k 1 is fixed. In practice, we can hesitate between several k 1 -values, denoted

k 1 1 , ..., k Λ 1 , Λ ∈ N, which means that for each k (j) 1 , j = 1, ..., Λ, we obtain a k 2 -value de- noted k (j) 2 . For each couple (k (j) 1 , k (j)
2 ), we choose the k (j) 2 -value, for which the relationship (4.5) is the "'most stable"', and therefore, we compute the criteria (4.5) for d successive values of k 2 taken around k (j) 2 . In that case the parameter k 2 is obtained minimizing:

φ(k (j) 2 , d) = d i=1 M (j) 1 -M * (j) i 2 , j = 1, ..., Λ, (4.7)
where, M

(j) 1 ≡ AM SE(k (j)
2 ) is obtained from (4.5) for each k

(j)
2 , and M * (j) i is computed using the i value around k 3. If we cannot interpret the Hill plot (Step 1), either because the data set is too unstable or too stable, we cannot use the previous algorithm, therefore, we use the whole information set X, and follow the method developed in [START_REF] Danielsson | Using a bootstrap method to choose the sample fraction in tail index estimation[END_REF] to estimate the ξ parameter (this method does not require an a priori value of the threshold). We estimate ξ by:

ξn (k) = 1 k k i=1 (log X (n-i+1) -log X (n-k) ) 2 2 ξn (k) . (4.8)
For a certain 5 n 2 , we bootstrap ξn (k) which appears in relationship (4.8) following the step 3 of the previous algorithm to create ξn 2 (k 2 ), where the parameter k 2 is obtained minimizing,

AM SE(k 2 ) := E( ξn 2 (k 2 ) -ξn 2 (k 2 )) 2 .
(4.9)

Then for a given k 2 , u is equal to:

u = (k 2 ) 2 k 3 (log k 2 ) 2 (2 log n 2 -log k 2 ) 2 log n 2 -log k 2 log n 2 , (4.10)
where k 3 is obtained exactly as k 2 with n 3 = n 2 2 n .

4. [START_REF] Shevchenko | Modeling operational risk data reported above a timevarying threshold[END_REF] suggest the use of a time varying threshold. In the current paper, we use a static methodology, nevertheless it would be interesting to integrate dynamics in the process of estimation. To do so, a possibility would be to use the dynamics of the severity data set given for instance by the incidents' amounts and their occurrence dates. In that case, the results would take into account the behavior changes in the risk faced through time and adjust with the state of nature they are currently in (crises, boom etc.).6

Capital Requirement

We applied the previous methodology to real operational risks data. We use two data sets denoted X (1) and X (2) .

• The set X (1) represents the severity of the business line "Retail Banking" and the event type "Business Disruption & System Failure and Execution" for the year 2007.

• The set X (2) represents the severity of the business line "Payment & Settlement" and the event type "Delivery, Execution and Process Management" for the year 2006.

Their statistics are given in Table 1. We deduce that their corresponding empirical severity distributions are asymmetric (right skewed) and leptokurtic. This statement strengthens our decision to implement a method which cares about extreme values in the distributions. 1: First four moments of the two data sets X (1) and X (2) . The third order empirical moment is non null and the fourth order empirical moment is far from 3.

Moments

In order to calibrate the appropriate loss distribution function on each data set, we use a Poisson distribution to model the frequency data. In the current paper, we did not challenge the choice of a Poisson distribution to model the frequency as the gain obtained from a continuous approximation or the use of another theoretical distribution has almost negligible impact on capital computation. Therefore, estimating the λ parameters of the Poisson distribution by maximum likelihood method, we get λ (1) = 2653 and λ (2) = 1292 respectively on the sets X (1) and X (2) .

To fit the severities, we tried several usual distributions suggested in the literature, but the Kolmogorov-Smirnov test rejected all (Table 6). Therefore, we continue our exercise and focus on the estimation of the threshold. To estimate the threshold u of the GPDs in the relation (3.4), we define two new sets X

(1) 1

an X

(2)

1 , following the notations introduced in the previous section. Their sizes are respectively n (2) 2 = 97. The resulting thresholds are u (1) = 1645.07 and u (2) = 179. They imply respectively 191 and 315 data above these thresholds.

We denote X (GP D 1 ) and X (GP D 2 ) the two corresponding data sets on which we fit the GPDs.

We estimate the parameters ξ and β of the GPDs defined in (3.4) using the method introduced by (Luceno ( 2006)) given in (3.5), denoting this method M1. We also consider three other alternative estimation methods to estimate these parameters in order to check their impact on VaR computations. These ones are respectively the Pickands method (M2) [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF]), the Hill method (M3) [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]), and the Maximum Likelihood method (M4).

We provide in Table 2 the estimations of GPDs parameters (with their standard deviation in brackets) obtained for both sets, using these four methods. Then we estimate the parameters of the lognormal distributions for the sets X (1) -X (GP D 1 ) and X (2) -X (GP D 2 ) and the results are provided in Table 3.

Method

X (1) X (2) β ξ β ξ 2: Estimations of the GPD's parameters ξ and β given in (3.3), using four methods respectively for u (1) = 1645.07 and u (2) = 179. We provide in brackets the standard deviations computed by bootstrapping.

M1
X (1) -X (GP D 1 ) X (2) -X (GP D 2 ) µ = 5.702144 σ = 1.103373 µ = 3.593098 σ = 1.510882 (0.07134) (0.03016) (0.04112) (0.02576)
Table 3: Estimations of the parameters of the lognormal distributions fitted on the data sets X (1) -X (GP D 1 ) and X (2) -X (GP D 2 ) . Standard deviation are provided in brackets.

As soon as all the parameters have been estimated, we implement a Monte Carlo algorithm with a million of iterations to build G, and we use it to compute a 99.9% VaR which provides the required amount of capital. Table 4 provides the VaRs of the set X (2) using the four estimation methods to estimate the threshold. We observe miscellanous amounts, varying from 5 725 341 euros to 538 480 990 euros. This last result is obtained from a distribution whose mean is infinite (the parameter ξ is greater than 1). This result points that the choice of the distribution and the estimation procedure are crucial for the VaR computations. The M1 method has the advantage to maximize the goodness-of-fit of the GPD and insures that the distribution cannot be statistically rejected. Therefore, one can argue that the VaR at 15 700 112 euros could be a reasonnable value considering the information set we have, and the assumptions done.

We do not provide the VaRs for the set X (1) as the results we obtained imply a simulated VaR multiplied by more than one hundred compared to the historical one. Here again, we can observe that the underlying model has an infinite mean.

Analyzing the results and controlling every step of the method, we concluded that the difficulty to obtain a "correct" and a "robust" model comes from the structure of the data set. Considering that different kinds of data are mixed in a single set it appears difficult to fit an accurate distribution on the severities. In order to analyze the impact of the granularity strategy, it would be interesting to split the data in multiple subsets whether the information sets sizes are large enough. In that latter case the computed VaRs would reflect risks created by specific kinds of risk, and finally aggregating them we could compare the results to the largest level of granularity. 4: VaR estimations for u = 179, µ = 5.681191 and σ = 1.081609 for the data set X (2) .

Estimations have a tremendous impact on the VaRs. In this table, these ones may differ by 9405%.

Finally we want to notice the influence of the threshold on the estimation of the GPD's parameters β and ξ. This effect is due to the fact that the information set changed. We illustrate this point providing estimates for five threshold values, for the set X (1) ( thresholds u for the data set X (2) . Parameters are very sensitive to the threshold value.

Conclusion

In this paper, we focus on the estimation of the LDF used to estimate operational risks. We highlight that the VaR computation depends on the method implemented. Therefore, a sharp analysis of the LDF is essential, and even if the LDF is a convolution of two distinct distributions, the capital requirements seem particularly sensitive to the choice and the estimation of the severity distribution. We present a new method to compute the severity distribution, splitting it in two parts in order to better take into account the large losses. We use a GPD on the right tail for which we provide innovative theoretical and practical solutions, and fit a lognormal distribution on the remaining data. Then, to build the final LDF, we apply a new adapted Monte Carlo algorithm. This method enables capturing at least two distinct behaviors in a severity distribution, and sometimes a bit more if the multimodality is not too strongly marked.

We point out the importance of the Basel matrix construction and the so-called granularity effect.

Second, we underline the fact that once the threshold of the GPD has been found, the method chosen to estimate GPD's parameters tremendously impacts the VaR: it seems that this fact has never been discussed before.

In our applications, we have been confronted to an important problem arising from the data bases construction. Following the regulator, who demands to take into account the right tail of the severity distribution to estimate the LDF, we sometimes obtained unrealistic VaRs (i.e.

equal to half the bank capitalization) due to infinite mean models. On the other hand, if we do not care about the severity tail, we have negligible VaRs. Then, even if it seems reasonable to follow regulator requirements, we point out two important facts: either the Bank are overexposed to these risks or the data sets are badly built. We privilege the last point and ask operational risks managers to pay attention to the large diversity of data origins. The correct question is to propose another way for building the data sets used to compute operational risks taking into account the difference which exists inside the different risk categories. It seems that it has no sense to mix some risks together. In our knowledge, this problem has not yet been discussed and the regulators need to reexamine risk typology and data sets granularity on which are computed operational risk capital requirements.

The method we provide enables to compute capital requirements with respect to Basel accords, so, we strictly stick to them. In this paper we focused on the effectiveness of the proposed method on an operational standpoint, however, other research tracks can extend the results developed in this paper. Even if in this paper we focus on the VaR computation to obtain the corresponding capital requirement, it will be interesting to use other risks measures, for instance the Expected Shortfall (ES) measure, and analyze their impact on the computation of the associated capital requirement. The use of the ES is not required by the regulator but it would be interesting because it is a coherent risk measure [START_REF] Artzner | Coherent measures of risk[END_REF], Extensions are proposed in [START_REF] Guégan | A mathematical resurgence of risk management: an extreme modeling of expert opinions[END_REF]). On the other hand there are alternatives to the Monte Carlo method to calculate the LDF [START_REF] Luo | Computing tails of compound distributions using direct numerical integration[END_REF], [START_REF] Guégan | A modified panjer algorithm for operational risk capital computation[END_REF]) permitting to improve the computation speed that could be investigated. Finally, to take into account the effect of the tail behavior, we focus in this paper on the Generalized Pareto Distribution, nevertheless other alternatives could be considered, e.g. the extreme value distributions defined in the Fisher-Tippett theorem [START_REF] Fisher | Limiting forms of frequency distributions of the largest or smallest member of a sample[END_REF]), the g-and-h distributions, or multivariate methods through the extreme value copulas [START_REF] Gudendorf | Extreme-value copulas[END_REF]). All these subjects will be discussed in companion papers.

B An adapted Monte Carlo algorithm

1. We simulate n realizations of the frequency distribution, q 1 , ..., q n , n ≥ 1000000.

2. For each q i , i = 1, • • • , n, we simulate q i values of the central distribution, (s 1 , ..., s q i ), and compare them to the threshold:

(a) If s j < u, ∀j ∈ [1, q i ] we keep them.

(b) If s j ≥ u, we count the number of exceedences, (z 1 , ..., z n ) for each frequency q i and for each i, we draw z i realizations of the GPD, (s

(GP D) 1
, ..., s

(GP D) z i
).

3. Finally, we have G i = q i j=1 (s j + s

(GP D) j
), a realization of the LDF, i = 1, ..., n.

4. The set (G 1 , ..., G n ) provides an empirical approximation of the LDF.

  exhibit the histogram of the Historical LDF. The black line corresponds to a LDF mixing a Poisson distribution and a lognormal severity. The associated 99.9% VaR is pointed out by ∆. The dashed line represents a LDF mixing a Poisson distribution and a multiple pattern severity distribution (lognormal-GPD given by (3.4)) using the Monte Carlo algorithm described in Appendix B. The corresponding VaR is represented by a +. The right graph focuses on the right tail of the LDFs described in this paragraph. This figure highlights the fact that this method enables thickening up the right tail of the LDF.

Figure 1 :

 1 Figure 1: Method Illustration: This figure presents the histogram (in grey) of the Historical LDF. The black line corresponds to a LDF mixing a Poisson distribution and a lognormal severity. The associated 99.9% VaR is pointed out by ∆. The dashed line represents a LDF mixing a Poisson distribution and a multiple pattern severity (lognormal-GPD given by (3.4)) using the algorithm discribed in Appendix B. The corresponding VaR is represented by a +. The right graph focuses on the right tail of the LDFs described above. This figure highlights the fact that this method enables thickening up the right tail of the LDF.

2 ,

 2 i = 1, ..., d.

  . When we bootstrap the data, we used respectively two sets of size n

Table 5 :

 5 Table5). All the estimated values are relevant despite a great instability. Thus the corresponding VaR and capital amount would be different. All steps of the previous methodology: estimation of u, estimation of the parameters ξ and β, methods of estimation play an important role in the computation of the This table presents GPD's parameters estimates and their standard errors for different

	capital requirement.
	Threshold u	β	Standard error β	ξ	Standard error ξ
	1162.12	1110.8538	87.77	0.6026	0.06877
	1608.27	946.6787	94.68	0.7924	0.09329
	1645.07	904.087	92.31	0.827	0.09713
	2021.48	1188.2956	155.48	0.9149	0.1397
	2177.58	1385.5295	199.56	0.8879	0.14

We used data from BPCE on the French Caisse d'Epargne perimeter.

The business lines are corporate finance, trading & sales, retail banking, commercial banking, payment and settlement, agency services, asset management and retail brokerage. The event types are internal fraud, external fraud, employment practices & workplace safety, clients, products & business practices, damage to physical assets, business disruption & system failures and execution, delivery & process management.

The break is determined empirically by finding a change in the slop.

We are only working on a part of the initial sample, it is necessary to adjust k2 to the whole data set.

A short algorithm providing the optimal n2 is given in[START_REF] Danielsson | Using a bootstrap method to choose the sample fraction in tail index estimation[END_REF].

Nevertheless in the current paper, this methodology cannot be considered because we do not have the corresponding information set. Besides, in[START_REF] Guégan | A mathematical resurgence of risk management: an extreme modeling of expert opinions[END_REF] we dealt with a dynamic threshold using a extreme distribution[START_REF] Fisher | Limiting forms of frequency distributions of the largest or smallest member of a sample[END_REF]) with another data set.

We do not provide those results as we did not have the corresponding data sets.

A Traditional distributions

X (2)

Table 6: Goodness-of-fit of the presented theoretical distributions to the sets X (1) and X (2) : the parameters have been estimated by maximum likelihood and we have chosen the Kolmogorov-Smirnov statistic. The p-values are all lower than 1%. It's not satisfactory. Considering goodness-of-fit test results, it is not possible to state that a distribution is better than another to model the losses.


Michèle Sebag

Olivier Teytaud

Combining Myopic Optimization and Tree Search: Application to MineSweeper

Many reactive planning tasks are tackled by optimization combined with shrinking horizon at each time step: the problem is simplified to a non-reactive (myopic) optimization problem, based on the available information at the current time step and an estimate of future behavior, then it is solved; and the simplified problem is updated at each time step thanks to new information. This is in particular suitable when fast off-the-shelf components are available for the simplified problemoptimality stricto sensu is not possible, but good results are obtained at a reasonnable computational cost for highly untractable problems. As machines get more powerful, it makes sense however to go beyond the inherent limitations of this approach. Yet, a brute-force solving of the complete problem is often impossible; we here propose a methodology for embedding a solver inside a consistent reactive planning solver. Our methodology consists in embedding the solver in an Upper-Confidence-Tree algorithm, both in the nodes and as a Monte-Carlo simulator. We show the mathematical consistency of the approach, and then we apply it to a classical success of the myopic approach: the MineSweeper game.

Introduction

Solving a reactive planning problem is a very hard task; e.g. playing optimally at the MineSweeper game is very hard, and planning power plants maintenances [START_REF] Roadef-Challenge | A large-scale energy management problem with varied constraints[END_REF] is usually done in a similar approximate manner: a move is decided based on current knowledge, with no tree search. This is quite reasonnable as too many uncertainties make tree search unreliable; however, it makes sense to look for the best of both worlds. As computers get stronger and stronger, it makes sense to work on the "real" problem using the approximate model as a preliminary guide. For this, and following [START_REF] Rolet | Optimal active learning through billiards and upper confidence trees in continous domains[END_REF][START_REF] Rolet | Optimal robust expensive optimization is tractable[END_REF] (combining heuristics with Upper Confidence Tree), we propose Upper Confidence Tree (UCT) combined with solvers, for progressively boosting a fast solver (constraint satisfaction problems) to an optimal behavior.

UCT is well known for its ability to work on problems with no or almost no expert information; in particular, it does not need any evaluation function. But is quite convenient for using prior knowledge encoded in heuristic strategies.

These heuristic strategies can be used both in the tree part of the algorithm or in the leaf evaluation; in particular, in the leaf evaluation, it takes into account long-term effects by simulating complete runs (as in T D(1) methods, see [START_REF] Bertsekas | Neuro-dynamic Programming[END_REF] for more on this).

MineSweeper is solved until 4x4 [START_REF] Nakov | Minesweeper, #minesweeper[END_REF], and is NP-complete in the general case [START_REF] Ben-Ari | Minesweeper as an np-complete problem[END_REF]. NP-completeness is not that impressive in games, as EXP or 2EXP complete games are not rare; but, in the partially observable case, we have not so many algorithms which give good results in practice. It is widely known that counting the number of solutions to a MineSweeper problem such that a mine appears in a given location is an efficient tool for playing the game; it is directly proportional to the probability of death if playing in this location. Heuristically, one can use such a method (playing a move with minimum probability of death) for approximately solving MineSweeper; however, it is known [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF] that such a strategy is suboptimal, i.e. this is not an optimal strategy, whatever may be the computation time.

We use the following classical terminology:

-A Markov Decision Process (MDP) is a set of states, with edges labelled with actions and probabilities; when a run is in state s and an action a is chosen, the next state is s ′ with probability p if there is an edge from s to s ′ labelled with a and p. The sum of the probabilities on out-edges of a state s labelled with a same action a must be equal to 1. A strategy for the MDP maps a finite sequence of states to an action (possibly in a stochastic manner). There are leaf-states (with no out-edge), and one state distinguished as the initial state; there are also rewards on edges, so that the overall Markov Decision Process, together with a strategy, leads to a probability distribution on rewards. -A Partial Observable Markov Decision Process (POMDP) is similar to a MDP, except that each state is equipped with an observation. A strategy is not a function from a finite sequence of states to an action, but a function from a finite sequence of observations to actions. -The state in a POMDP can refer to the "real" (partially hidden) state, or to the observed part only; we will refer to the "real" partially hidden state as the "complete state" in order to avoid confusion, in particular because in the MineSweeper literature the "state" usually refers to the observations. -The feasible set is, given a sequence of observations, the set of possible complete states. -A consistent belief state estimation is an algorithm which exactly samples the uniform probability distribution on the feasible set, given the observation. -An asymptotically consistent belief state estimation is an algorithm which samples the probability distribution of the hidden state with a distribution converging to the real one as the computational effort runs to infinity.

Section 2 presents the MineSweeper game. Section 3 will present the algorithms analyzed or proposed in this paper (Section 3.4 is a newly designed algorithm for combining belief state estimation and tree search). Section 4 will discuss mathematical aspects. Section 5 will explain why the MineSweeper game is interesting and present experimental results on it.

As UCT is now a well established part of the game literature, we refer to [START_REF] Coulom | Efficient Selectivity and Backup Operators in Monte-Carlo Tree Search[END_REF][START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF][START_REF] Lee | The Computational Intelligence of MoGo Revealed in Taiwan's Go Tournaments[END_REF] for more information on it; we only briefly summarize UCT in Fig. 1. Our -Initialize the memory to one single node, representing the current state; this is the root of the tree that will grow in memory. -While ( there is time left), simulate an episode as follows:

• First part of the episode: * Increase (by r) the counter of the total reward. -Play the action which is played most often from the root in episodes simulated above.

*
Fig. 1. Overview of the original UCT algorithm. There are other variants of the rule for choosing the actually made action (last line). There are also plenty of improvements over the initial UCB formula in the episodes; we also use progressive widening and an improved Monte-Carlo part.

UCT is as specified in Table 1. We use double progressive widening, as explained in [START_REF] Couetoux | Continuous Upper Confidence Trees[END_REF].

The MineSweeper game

The MineSweeper game exists on many platforms. At each turn, the player chooses a location in a grid. After the first move, M mines are randomly put on the h×w board (anywhere except on the chosen location); the player is informed of the number of mines in the 8-neighborhood of his move. From now on, when a player plays in location l, there are two cases:

-There is a mine: then, the game is over.

-There is no mine: then, the player is informed of the number of mines in the 8-neighborhood. Otherwise, with probability 0.7

Move with minimum probability of mine (computed by CSP)

Otherwise

Randomly draw a hidden state and play a move with no mine for this hidden state (computed by CSP) Bandit formula Upper Confidence Bound with variance estimates[1] Forced moves [START_REF] Teytaud | On the Huge Benefit of Decisive Moves in Monte-Carlo Tree Search Algorithms[END_REF] Points with null probability of mines (estimated by CSP) Table 1. Our UCT implementation in one table. Our work is on the belief estimation part and our UCT is just a classical variant. The single point strategy (from PGMS) is detailed later (Section 3.4). CSP refers to Constraint Satisfaction Problems, and SPS refers to the Single Point Strategy.

The game is a win if the player plays the h×w-M locations with no mine without ever playing a move on a mine. In usual implementations, when a move is played with no mine in the 8-neighborhood, then all the neighbours are automatically played as they are for sure secure moves; this just saves up time.

All our experiments are performed on the MineSweeper game; the MineSweeper game is more complicated than expected at first view [START_REF] Kaye | Minesweeper is np-complete[END_REF][START_REF] Ben-Ari | Minesweeper as an np-complete problem[END_REF][START_REF] Nakov | Minesweeper, #minesweeper[END_REF]. MineSweeper has motivated a lot of research, sometimes for pedagogical reasons [START_REF] Becker | Teaching with games: the minesweeper and asteroids experience[END_REF] (including a nice version aimed at teaching quantum mechanics [START_REF] Gordon | Quantum computer games: quantum minesweeper[END_REF]), or because it is a widely spread game (anyone here who has never played a MineSweeper game ?), or as a model for real problems [START_REF] Hein | Minesweeper for sensor networks-making event detection in sensor networks dependable[END_REF], and as a challenge for machine learning [START_REF] Castillo | Learning minesweeper with multirelational learning[END_REF] or genetic programming [START_REF] Koza | Genetic Programming II: Automatic Discovery of Reusable Programs[END_REF].

Minesweeper artificial intelligences fall in two categories:

-Algorithms based on the understanding of the underlying Markov Decision Process, which is quite hard to analyze and solve [START_REF] Nakov | Minesweeper, #minesweeper[END_REF]; -Algorithms based on heuristically choosing the move with lowest probability of immediate death by Constraint Satisfaction Problems; these algorithms are probably the most impressive for real board size and time settings.

The classical approach for playing MineSweeper consists in using Constraint Satisfaction Problems [START_REF] Studholme | Minesweeper as a constraint satisfaction problem[END_REF] (CSP). It provides a consistent estimate of the probability of immediate death in case of play in a given location; minimizing this is a quite good heuristic. However, this is not optimal [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF]; we here propose an approach which is asymptotically optimal by using UCT (a consistent MDP solver) combined with a consistent belief state estimation by a CSP algorithm. Importantly, this work in based on [START_REF] Nakov | Minesweeper, #minesweeper[END_REF], which shows the relevance of CSP as a consistent belief state estimator. Using CSP inside UCT will require more than just the computation of the probability of presence of a mine in each location -we will need the complete probability of a given complete state given observations.

The best published results are usually those of the CSP approach, with 80% (beginner: 9x9, 10 mines), 45% (intermediate: 16x16, 40 mines); and 34% in expert mode (16x30, 99 mines). [START_REF] Pedersen | The complexity of minesweeper and strategies for game playing[END_REF] proposed a so-called limited search method which gets a earthshaking 92.5% (beginner), 67.7% (intermediate).

Following most previous works, we work on a version of the game in which the first move can not be on a mine (the mines are randomly drawn until they are not at the first chosen location). We do not use the variant in which the initial move is on a "zero" location (as the Linux Gnomine version); this looks like a very good idea for reducing the probability of unlucky death, but it is not widely used and therefore we prefer to be consistent with the body of work around MineSweeper.

We point out the importance of the initial move; the CSP approach does not say anything on it, whereas it is critical; we believe that this explains the success of some implementations -even if the algorithm is not better than the CSP approach, and if the authors have implemented a good initial move or a good opening book, then they might get better results than the CSP approach. However, we point out that the initial move is not the only suboptimality of CSP; as shown in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF], many small patterns lead to suboptimal moves with non-zero probability.

Algorithms

Many algorithms for solving POMDP are based on two tools, more or less separated, namely:

the estimation of the belief state (i.e. the probability distribution on the complete state, including the hidden part); -the choice of a move, based on this estimated belief state.

These two components are more or less separated, depending on the approaches. In CSP for MineSweeper, the second part is trivial (play the move with lowest probability of being a mine) whereas the first part is exact. Our goal is to do a non-trivial second part. For the second component, we will use a Monte-Carlo Tree Search; following [START_REF] Rolet | Optimal active learning through billiards and upper confidence trees in continous domains[END_REF], we will generate transitions conditionally to past observations (we need a complete forward model for a UCT implementation).

We here focus on the first component, namely belief state estimation, and discuss several approaches:

-A simple rejection method (Section 3.1); -Constraint Satisfaction Problem (CSP, Section 3.2); this is an optimal belief state estimation algorithm, but to the best of our knowledge it has never been used inside a tree search; as a consequence, existing CSP algorithms are not optimal; we will here use it as a belief state estimation rather than as a myopic heuristic; -Naive Constraint Solving plus Markov-Chain Monte-Carlo (Section 3.3); this was tested in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF] but the authors do not report better results than the simple rejection method; -BSSUCT (Belief State Solver + UCT), a new method presented in Section 3.4.

Rejection method

The rejection method is the simplest tool for estimating the complete state, given the observable part. It is consistent, i.e. it proposes an exact belief state estimation; it can be slow, but it will use, by definition, a computation time sufficient for ensuring a correct estimate. The pseudo-code is given in Alg. 1.

In spite of its simplicity, it was reported to be the best method in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF]; the interpretation of this result in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF] is that rejection was the only method which was consistent (The Markov-Chain Monte-Carlo method is asymptotically consistent only). Following this comment, we propose a method which is consistent, as well as the rejection method, but faster, cf Section 3.4.

Algorithm 1

The rejection method for consistent belief state estimation. C is a parameter such that C × Likelihood observation (s) ≤ 1; if the likelihood of observations, given hidden state, is binary (0 or c), as in MineSweeper (c = 1 for MineSweeper), then C = 1/c is the optimal choice; it just means that we accept the state if and only if it is consistent with observations. Input: observations. Output: a (uniformly sampled among consistent states) complete state. s ←random complete state while random ≥ C × Likelihood observations (s) do s ←random complete state end while Return s

Constraint Satisfaction Problem (CSP)

We here refer to the many implementations based on counting all complete states which are consistent with observations, i.e. the entire feasible set. This provides the exact probability, for each location, to be a mine. The CSP algorithm can be used for exactly estimating the probability of transition to a given state; yet, to the best of our knowledge, it has only be used for estimating the probability for each location to be a mine and for playing this move. Such an algorithm is presented in Fig. 2.

Algorithm 2

The CSP algorithm for playing MineSweeper. The algorithm is intrinsically myopic: it only optimizes the 1-step result (minimum probability of immediate death).

CSP-function for choosing a move on a partially visible board with M mines for each location l of the board do N b(l) ← 0 end for for each positioning p of the M mines on the board consistent with observations do for each location l with a mine for p do N b(l) ← N b(l) + 1 end for end for Play move l uniformly chosen among the set of uncovered locations l such that N b(l) is minimum.

Naive Constraint Solving plus Markov-Chain Monte-Carlo

Markov-Chain Monte-Carlo (MCMC) is a tool for asymptotically consistently estimating the belief state. Metropolis-Hastings (MH) is a classical version of MCMC. In the case of MineSweeper, or more generally a case in which observations are deterministic and binary as a function of the belief state, and with a symmetric transition kernel (i.e. the probability of mutating from state x to state x ′ is equal to the probability of mutating from state x ′ to state x), MH boils down to an algorithm as shown in Alg. 3.

Algorithm 3 A version of Metropolis-Hastings with symmetric transition kernel in the case of deterministic binary observations (as in MineSweeper). Please note that removing the initial constraint solving does not change the asymptotical properties of MH; this improved initialization is here for saving up time.

According to [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF] and for the mutation used there, the initial state chosen by constraint solving unfortunately leads to a too big bias and makes the overall algorithm weaker than the simple rejection method.

Input: observations. Output: a (nearly uniformly sampled) complete state. Find an initial state s by constraint solving, consistent with observations. Select a number T of iterations by heuristic methods // in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF], T depends on the number of UCT-simulations.

for t ∈ [[1, T ]] do Let s ′ be a mutation of s if s ′ is consistent with observations then s ← s ′ end if end for
Variants of this version were tested in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF], with disappointing results; the authors conclude that the asymptotic consistency of this approach (instead of "real" consistency of the rejection method) was the reason for the moderately good results. Obviously, it is hard to conclude on this point without testing rigorously the infinitely many possible approaches (the many parameters) for MCMC; maybe, there is a version which would make Alg. 3 extremely efficient. However, we want to have a simple methodology, as far as possible independent of the problem; therefore, we will propose in Section 3.4 a rigorously nonasymptotically consistent belief state estimation.

BSSUCT: Belief State Solver + Upper Confidence Trees

Our variant of UCT is detailed in Table 1. Our UCT uses several modules:

the single point strategy as in PGMS (http://www.ccs.neu.edu/home/ ramsdell/pgms/). The basic idea is that if k of the neighbours of a location are mines, and there are n neighbours, m of which are non-mines and p of which are mines and nmp are uncovered, then:

• If k = n -m -p,
then all uncovered neighbours are mines;

• If k = p, then all uncovered neighbours are not mines. Propagating this is usually done in MineSweeper interfaces in the case k = 0 (all neighbours are then automatically played without human action for fastening the gameplay). Implementing the general case is easy and can't hurt (this plays riskless moves only).

-CSP. Constraint Satisfaction is applied for finding all possible hidden states.

We use least constraint variables first -whereas most constrained variables should be used first when looking for one single solution, it is better to use least constrained variables first when looking for all solutions. CSP, in reasonnable time, can provide many information:

• the probability, for each location, that there is a mine;

• as a consequence, it can point out sure moves, to be played with no risk;

• the exact probability distribution of the next state (i.e. how will be the board after next move). With this, the POMDP becomes a MDP.

The use of these modules is detailed in Table 1.

Theoretical analysis

We show the consistency of our approach (section 4.1) in the generic case and the non-consistency of the classical approach (constraint satisfaction problems) for the MineSweeper game (section 4.2).

Consistency of the approach

We first discuss the consistency of our UCT approach as a MDP solver. UCT (Upper Confidence Trees) is a consistent MDP solver [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF] in the finite case. Our approach uses heuristics in the evaluation of leafs by Monte-Carlo simulations, but it does not matter for the consistency proof in [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF]. We also use heuristics in the tree part of UCT; this is not detailed in [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF], but the upper-confidence-bound proof from [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF][START_REF] Auer | Using confidence bounds for exploitation-exploration trade-offs[END_REF], used in [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF], is also independent of heuristics. A last component which differs from [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF] is progressive widening [START_REF] Couetoux | Continuous Upper Confidence Trees[END_REF]; however, progressive widening here has only a non-asymptotic impact because the number of actions and the number of random outcomes is finite. So, if the problem is a MDP, our UCT-based approach is consistent. This shows that our approach is consistent for solving a MDP. We then discuss the consistency for the application to MineSweeper. The difficulty is that MineSweeper is not a MDP: there is partial information, and therefore it is a POMDP. However, [START_REF] Nakov | Minesweeper, #minesweeper[END_REF] has shown that the partial information can be exactly estimated by CSP. The classical approach by CSP is not optimal for playing when just used for the heuristic "play the move with lowest probability of immediate death" (as shown in [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF] there are simple situations in which this is not optimal), but it is an exact tool for rephrasing a POMDP as a MDP.

As a consequence, our algorithm is a consistent MDP solver, and combined with CSP it is an asymptotically optimal MineSweeper player.

We do validation experiments in the GnoMine Custom mode, as the GnoMine rule (the first move is a zero) makes the problem easier to analyze in a mathematical manner so that we can check the optimality of our algorithm. Then, we switch to the classical MineSweeper rule (the first move is not a mine). We compare our approach to the CSP-PGMS results.

A Gnomine Custom mode: 15 mines on a 5x5 board

The GnoMine version has this special rule that the location played first has no mine in the neighborhood. This implies that playing in the center leads to a sure win in the configuration 15 mines on a 5x5 board (see Fig. 3). Interestingly, our algorithm finds this optimal strategy (that humans easily find as well); but the important point here is that the CSP method does not say anything about the choice of the first move. If we play randomly and uniformly the first move, then we have a probability 4 25 of playing in (2,2) (or a symmetry of it), with then a probability of loosing the game at least 1 2 ; and we have a probability 4 25 of playing (2,3) (or a symmetry of it), with then a probability of loosing the game at least 1 4 . This leads to an overall probability of loosing the game at least

1 2 × 4 25 + 1 4 × 4 25 = 3
25 ; the real probability is indeed much bigger, as (1,1), (1,2), (1,3) are also very bad moves, but we do not compute the detailed probabilitythe lower bound 3 25 on the probability of loosing is enough for showing that the CSP can not reach the optimal play (that our BSSUCT provably reaches).

Given the rules only (no expert knowledge added), our algorithm asymptotically finds the best move, and in 5x5 with 15 mines and gnomine rules it finds it in 5 seconds (in all runs among 500 runs).

Windows version, various board sizes

We compare our algorithm (BSSUCT) with CSP on various board sizes; the CSP implementation we use is the CSP-PGMS implementation from http:// www.ccs.neu.edu/home/ramsdell/pgms/ which plays in corners for the initial moves. Results are presented in Table 2 The winning rate for CSP-PGMS is estimated on 100000 games; we provide the standard deviation for our program which is much slower and therefore has non-negligible confidence intervals. Importantly, CSP-PGMCS does not benefit from additional time. The games are played with 10s per move, except for 10 mines in 5x5 (300 seconds per move) and 10 mines in 3x7 (30s per move). All experiments on a 16-cores Xeon 2.93GHz Ubuntu 2.6-32-34 using 1 core only per experiment.

Conclusion

We have proposed the use of UCT for difficult decision tasks, with CSP for two simultaneous purposes:

estimating the belief state, so that a POMDP becomes a MDP; -suggesting good moves, in a heuristic manner. This provides a generic methodology for improving a myopic solver by including it into a UCT. We strongly believe that such approaches, using consistent asymptotic behaviors and myopic-guided heuristics, are important for future artificial intelligence, benefiting both from relevant heuristics (which provide tractability) and strong computational power and mathematical consistency analysis (which provide asymptotic optimality). Starting at a reasonnably good strategy and improving it towards optimality by UCT looks like a simple efficient idea.

We have shown the mathematical consistency of the approach, as well as experimental results. The experimental results are on MineSweeper, which is particularly challenging as the CSP approach, which completely neglects long term effects, is very effective as uncertainties are so big that long-term effects are very unreliable -we have chosen a very hard test case and we nonetheless get significant improvements. When the computation time becomes big, CSP does not improve anymore (it has inherent limitations), whereas our approach performs better and better.

For getting impressive results in the game community, we will work on a faster implementation; our CSP module is much slower than e.g. CSP-PGMS. It should be easy to work on expert modes, provided that we have an implementation as fast as CSP-PGMS for the CSP module. We have already clear improvements on small boards, compared to [START_REF] Couetoux | Consistent belief state estimation, with application to mines[END_REF]. As a future work, we will also compare with results of [START_REF] Pedersen | The complexity of minesweeper and strategies for game playing[END_REF].

A crucial good piece of news for our algorithm is that we do not enter any opening move manually (as necessary for good performance with CSP); and we can switch easily to variants of MineSweeper just by changing the implementation of the rules.

Fig. 3 .

 3 Fig.3. The Gnomine version of the 5x5 board with 15 mines is a sure win: on each of these figures, one can deduce the position of all mines (there is only one non-mine location). This covers all possibles cases by rotation.

Table 2 .

 2 . Winning rate of our implementation, versus the winning rate of CSP-PGMS.

	Format	CSP-PGMS	BSSUCT
	4 mines on 4x4	64.7 %	70.0% ± 0.6%
	1 mine on 1x3	100 %	100% (2000 games)
	3 mines on 2x5	22.6%	25.4 % ± 1.0%
	10 mines on 5x5	8.20%	9% (p-value: 0.14)
	5 mines on 1x10 12.93%	18.9% ± 0.2%
	10 mines on 3x7	4.50%	5.96% ± 0.16%
	15 mines on 5x5	0.63%	0.9% ± 0.1%

Acknowledgements

The authors are grateful to NSC for funding NSC100-2811-E-024-001, to ANR for funding COSINUS program (project EXPLO-RA ANR-08-COSI-004), to L. Simon for fruitful discussions. We are very grateful to L. Simon for very fruitful discussion, and to the authors of PGMS and PGMS-CSP for making their implementations freely available. We also acknowledge gratefully the support by European Program FP7, through the MASH project.

A simple case in which our algorithm is better than CSP-based approaches

Our algorithm is asymptotically consistent, i.e. finds the optimal strategy if given a sufficient thinking time (it is asymptotically consistent because UCT is asymptotically consistent for solving MDPs and the rejection method provides the right MDP as explained in Section 4.1). Therefore, it finds the optimal strategy for MineSweeper in 3x3 with 7 mines. We here show that CSP, playing uniformly the initial move, can not be optimal in such a case. Let's see the success rate, depending on the first move, in this small version of MineSweeper: -If we play in the center, we necessarily observe 7; the probability of winning is then 1 8 as there are 7 mines in the 8 neighbours (see Fig. 2). -If we play on a side, then we observe:

• a 5 with probability 3 8 ; then the probability of winning is 1 3 ; • a 4 with probability 5 8 ; then the probability of winning is 1 5 . This means an overall probability of winning at most

when starting on the side. -If we play in a corner, then we observe:

• a 2 with probability 3 8 ; then the probability of winning is 1 3 ; • a 3 with probability 5 8 ; then the probability of winning is 1 5 . This means an overall probability of winning at most

when starting in a corner.

An optimal strategy then consists in playing either on a side or in a corner, in order to get a probability 1 4 of winning. The CSP algorithm plays uniformly the first move, because all moves have the same probability of death; therefore, it has a probability of winning , which is less than the probability 1 4 of winning obtained by an optimal strategy (and asymptotically reached by BSSUCT).