Caroline Bernard-Michel

Laurent Gardes

Stéphane Girard
email: stephane.girard@inrialpes.fr

Gaussian Regularized Sliced Inverse Regression

Keywords: Inverse regression, regularization, sufficient dimension reduction. methods

published or not. The documents may come

Introduction

Many methods have been developed for inferring the conditional distribution of an univariate response Y given a predictor X in R p , ranging from linear regression [START_REF] Draper | Applied regression analysis (3rd edition)[END_REF] to support vector regression [START_REF] Cristianini | An introduction to support vector machines and other kernel-based learning methods[END_REF]. When p is large, sufficient dimension reduction aims at replacing the predictor X by its projection onto a subspace of smaller dimension without loss of information on the conditional distribution of Y given X. In this context, the central subspace, denoted by S Y |X plays an important role. It is defined as the smallest subspace such that, conditionally on the projection of X on S Y |X , Y and X are independent. In other words, the projection of X on S Y |X contains all the information on Y that is available in the predictor X.

The estimation of the central subspace has received considerable attention since the past twenty years. Without intending to be exhaustive, we refer to Sliced Inverse Regression (SIR) [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], sliced average variance estimation [START_REF] Cook | Discussion of "Sliced Inverse Regression for dimension reduction" by K.C. Li[END_REF], and graphical regression [START_REF] Cook | Regression graphics. Ideas for studying regressions through graphics[END_REF] 2 Inverse regression without regularization Consider X a R prandom vector, Y the real response variable and let us denote by S Y |X the central subspace. In the following, for the sake of simplicity, we assume that K = dim(S Y |X ) = 1, the case 1 < K < p being discussed in Section 3. We thus introduce β ∈ R p such that S Y |X = span(β). In Subsection 2.1, the considered inverse regression model is presented. The estimation of the unknown parameters is discussed in Subsection 2.2 and the links with the SIR method are established in Subsection 2.3.

Single-index inverse regression model

The following inverse single-index regression model is considered (see [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF], equation [START_REF] Antoniadis | Bayesian estimation in single-index models[END_REF] for the multi-index model):

X = µ + c(Y )V b + ε, (1)
where µ and b are nonrandom R pvectors, ε is a centered R p -Gaussian random vector, independent of Y , with covariance matrix Cov(ε) = V and c : R → R is a nonrandom function. Under this model, E(X|Y = y) = µ + c(y)V b and thus, after translation by µ, the conditional expectation of X given Y is a random vector located in the direction V b.

From [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF], Proposition 1, b corresponds to the direction β of the central subspace. In the sequel, it will appear that, under appropriate conditions, the maximum likelihood estimator of b is (up to a scale parameter) the SIR estimator of β. Moreover, note that, under [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF], one has

c(y) = E(b t (X -µ)|Y = y) b t V b . (2)
Now, restricting ourselves to the single-index case, the forward model of SIR asserts that there exists a univariate link function g such that E(Y |X) = g(b t X) or equivalently, b t X = g -1 (E(Y |X)). Thus, replacing in (2) yields c(y) = g -1 (y)b t µ b t V b , i.e. the coordinate function is, up to an affine transformation, the inverse of the link function in the single-index forward model of SIR.

Maximum likelihood estimation

We now address the estimation of the coordinate function c(.), the direction b, the covariance matrix V and the location parameter µ in model [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF]. To this end, we focus on a projection estimator of the unknown function c(.). More precisely, it is expanded as a linear combination of h basis functions s j (.), j = 1, . . . , h:

c(.) = h j=1 c j s j (.), (3)
where the coefficients c j , j = 1, . . . , h are unknown whereas h is supposed to be known. Introducing c = (c 1 , . . . , c h ) t and s(.) = (s 1 (.), . . . , s h (.)) t , model (1) can be rewritten as

X = µ + s t (Y )cV b + ε, ε ∼ N (0, V ), (4)
where N (0, V ) is the multivariate centered Gaussian distribution with covariance matrix V . In the sequel we denote by

ρ = var(b t (X -ε)) var(b t ε) = b t Σb -b t V b b t V b ,
the signal to noise ratio in the direction b. Let (X i , Y i ), i = 1, . . . , n be a sample of independent random variables distributed as (X, Y ). Clearly, estimating (µ, V, b, c) by maximization of the likelihood in the Gaussian model ( 4) consists in minimizing

G(µ, V, b, c) = log det V + 1 n n i=1 (µ + s t (Y i )cV b -X i ) t V -1 (µ + s t (Y i )cV b -X i ), (5)
with respect to (µ, V, b, c). Note that G(µ, V, b, c) can also be interpreted as a discrepancy functional, see equation ( 5) in [START_REF] Cook | Sufficient dimension reduction via inverse regression: A minimum discrepancy approach[END_REF]. Up to our knowledge, the introduction of such functional in the inverse regression framework is due to [START_REF] Cook | Testing predictor contributions in sufficient dimension reduction[END_REF]. Let us introduce the h × h empirical covariance matrix W of s(Y ) defined by

W = 1 n n i=1 (s(Y i ) -s)(s(Y i ) -s) t , the h × p matrix M defined by M = 1 n n i=1 (s(Y i ) -s)(X i -X) t ,
and Σ the empirical p × p covariance matrix of X

Σ = 1 n n i=1 (X i -X)(X i -X) t , where X = 1 n n i=1 X i and s = 1 n n i=1 s(Y i ).
Lemma 1 Using the above notations, G(µ, V, b, c) can be rewritten as

G(µ, V, b, c) = log det V + tr( ΣV -1 ) + (µ -X + st cV b) t V -1 (µ -X + st cV b) + (c t W c)(b t V b) -2c t M b.
The maximum likelihood estimators of µ, V , b and c are closed-form.

Proposition 1 Under (4), if W and Σ are regular, then the maximum likelihood estimator of (µ, V, b, c) is defined by:

• b is the eigenvector associated to the largest eigenvalue

λ of Σ-1 M t W -1 M , • ĉ = 1 bt V b W -1 M b, • μ = X -st ĉ V b, • V = Σ - λ bt Σb Σbb t Σ.
As a consequence of the above equation, one also has λ = 1 -bt V b/ bt Σb , which provides an estimation of the signal to noise ratio in the direction b through

ρ = λ/(1 -λ). (6)
Let us now show that the SIR method corresponds to the particular case of piecewise constant basis functions s j (.), j = 1, . . . , h.

Sliced Inverse Regression (SIR)

Suppose the range of Y is partitioned into h + 1 non-overlapping slices S j , j = 1, . . . , h + 1 and consider the h basis functions defined by

s j (.) = I{. ∈ S j }, j = 1, . . . , h (7)
where I is the indicator function. Let us denote by n j the number of Y i in slice j = 1, . . . , h + 1, define the corresponding proportion by f j = n j /n and the empirical mean of

X given Y ∈ S j by Xj = 1 n j Y i ∈S j X i .
The covariance matrix of the regression curve Γ = Cov(E(X|Y )) is then estimated by the p × p empirical "between slices" covariance matrix

Γ = h+1 j=1 f j ( Xj -X)( Xj -X) t .
In this context, the following consequence of Proposition 1 can be established.

Corollary 1 Under ( 4) and [START_REF] Cook | Discussion of "Sliced Inverse Regression for dimension reduction" by K.C. Li[END_REF], if Σ is regular, then the maximum likelihood estimator b of b is the eigenvector associated to the largest eigenvalue of Σ-1 Γ.

To summarize, the maximum likelihood estimator of b is (up to a scale factor) the SIR estimator of the direction β spanning the central subspace. The next section is dedicated to the introduction of a regularization in the inverse regression problem in order to avoid the inversion of Σ.

Regularized inverse regression

First, we present in Subsection 3.1 how the introduction of a prior information on the unknown direction b can overcome the SIR limitations due to the ill-conditioning or singularity of Σ. Second, some links with the existing SIR regularizations are highlighted in Subsection 3.2. Finally, three new regularizations of the SIR method based on our framework are introduced in Subsection 3.3.

Introducing a Gaussian prior

We propose to introduce a prior information on s t (Y )cb appearing in the projection of X on b under the inverse regression model. More precisely, we focus on

E(b t (X -µ)|Y )b b t V b = c(Y )b = s t (Y )cb, from (2)
and (3), assuming that

(1 + ρ) -1/2 (s(Y ) -s) t cb ∼ N (0, Ω). ( 8
)
The role of the matrix Ω is to describe which directions in R p are the most likely to contain b. Some examples are provided in the next two paragraphs. The scalar (1 + ρ) -1/2 is introduced for normalization purposes, permitting to preserve the interpretation (6) of the eigenvalue in terms of signal to noise ratio. As a comparison, in [START_REF] Antoniadis | Bayesian estimation in single-index models[END_REF], a Bayesian estimation method is proposed using a B-splines approximation of the link function g in the forward model and a Fisher-von Mises prior on the direction b. In our approach, working on the inverse regression model allows to obtain explicit solutions, see Lemma 2 and Proposition 2 below.

Lemma 2 Maximum likelihood estimators are obtained by minimizing

G Ω (µ, V, b, c) = G(µ, V, b, c) + (b t Ω -1 b)(b t V b)(c t W c) b t Σb (9)
with respect to (µ, V, b, c).

Comparing to [START_REF] Bernard-Michel | A note on Sliced Inverse Regression with regularizations[END_REF], the additional term due to the prior information can be read as a regularization term in Tikhonov theory, see for instance [START_REF] Vogel | Computational methods for inverse problems[END_REF], Chapter 1, penalizing large projections. The following result can be stated:

Proposition 2 Under ( 4) and [START_REF] Cook | Regression graphics. Ideas for studying regressions through graphics[END_REF], if W and Ω Σ + I p are regular, then the maximum likelihood estimator of (µ, V, b, c) is defined by:

• b is the eigenvector associated to the largest eigenvalue λ of (Ω Σ+I p ) -1 ΩM t W -1 M ,

• ĉ = 1 (1 + η( b)) bt V b W -1 M b, where η( b) = bt Ω -1 b bt Σb , • μ = X -st ĉ V b, • V = Σ - λ bt Σb Σbb t Σ.
The inversion of Σ in Proposition 1 is replaced by the inversion of Ω Σ+I p in Proposition 2. Thus, for a properly chosen prior matrix Ω, the numerical instabilities in the estimation of b disappear. Note that, since the estimation of V is formally the same as in Proposition 1, the interpretation (6) of λ still holds. As previously, this result can be applied to the particular case of the SIR method.

Corollary 2 Under ( 4), ( 7) and ( 8), if Ω Σ + I p is regular, then the maximum likelihood estimator of b is the eigenvector associated to the largest eigenvalue of (Ω Σ + I p ) -1 Ω Γ.

In the following, the above estimator of the direction b will be referred to as the Gaussian Regularized Sliced Inverse Regression (GRSIR) estimator. Let us emphasize that the GRSIR estimator can be extended to the multi-index situation by considering the eigenvectors b1 , . . . , bK associated to the

K largest eigenvalues λ1 > • • • > λK of (Ω Σ + I p ) -1 Ω Γ.
For instance, one can show that b2 maximizes G Ω under the orthogonality constraint bt 1 ( Σ + Ω -1 ) b2 = 0.

Links with existing methods

In all the next examples, a non-negative parameter τ is introduced in the prior covariance matrix in order to tune the importance of the penalty term in [START_REF] Cook | Testing predictor contributions in sufficient dimension reduction[END_REF]. Consequently, in the sequel, τ is called a regularization parameter.

Classical SIR approach. It is easily seen from Corollary 2 that choosing the prior covariance matrix Ω 0 = (τ Σ) -1 in GRSIR gives back SIR, and this for all τ > 0. This prior matrix indicates that directions corresponding to small variances are most likely, i.e the SIR method favors directions in which Σ is close to singularity. In practice, this choice yields instabilities in the estimation.

Ridge approach. The simplest choice for the prior covariance matrix is Ω 1 = τ -1 I p . In this case, the identity matrix indicates that no privileged direction for b is available. Following Corollary 2, the corresponding GRSIR estimator of b is the eigenvector of ( Σ + τ I p ) -1 Γ associated to its largest eigenvalue, which is the ridge estimator introduced independently in [START_REF] Zhong | RSIR: Regularized Sliced Inverse Regression for motif discovery[END_REF] and [START_REF] Scrucca | Regularized Sliced Inverse Regression with applications in classification[END_REF][START_REF] Scrucca | Class prediction and gene selection for DNA microarrays using regularized Sliced Inverse Regression[END_REF].

PCA+SIR approach. As already seen, a popular technique to overcome the singularity problems of Σ is to use PCA as a preprocessing step [START_REF] Chiaromonte | Dimension reduction strategies for analyzing global gene expression data with a response[END_REF][START_REF] Li | Dimension reduction methods for micro-arrays with application to censored survival data[END_REF]. The principle is the following. Let d ∈ {1, . . . , p} be fixed and denote by δ1 ≥ • • • ≥ δd the d largest eigenvalues of Σ (supposed to be positive), q1 , . . . , qd the associated eigenvectors and S d = span(q 1 , . . . , qd )

the linear subspace spanned by q1 , . . . , qd . The first step consists in projecting the predictors on S d . The second step is to perform SIR in this subspace. The next result shows that this method corresponds to a particular prior covariance matrix.

Proposition 3 PCA+SIR corresponds to GRSIR with prior covariance matrix

Ω 2 = 1 τ d j=1 1 δj qj qt j ,
where τ > 0 is arbitrary.

Let us note that although Ω 2 depends on τ , the GRSIR estimator does not, since there is no regularization parameter in the PCA+SIR methodology.

L 2regularization. The regularization proposed in [START_REF] Li | Sliced Inverse Regression with regularizations[END_REF] consists in estimating (b, c) by minimization of

H τ (b, c) = h j=1 f j (µ + Σc j b -Xj ) t N (µ + Σc j b -Xj ) + τ b t b,
the matrix N being either N = I p or N = Σ-1 . In our opinion, this approach suffers from a lack of invariance since the functional H τ does not penalize the same way two different axes (b and 2b for instance) defining the same direction. As a consequence, one can show ( [START_REF] Bernard-Michel | A note on Sliced Inverse Regression with regularizations[END_REF], Proposition 1) that the only possible solution b of the minimization problem is b = 0. In view of this result, the proposed alternating least squares algorithm ( [START_REF] Li | Sliced Inverse Regression with regularizations[END_REF], Section 2) cannot be justified theoretically. As a comparison, our method does not yield this kind of problem thanks to the invariance property: G Ω (µ, V, tb, c/t) = G Ω (µ, V, b, c) for all real numbers t = 0. We now propose some alternative choices of the covariance matrix Ω yielding new regularizations of the SIR method.

Three new SIR regularizations

Tikhonov regularization. An alternative choice of the prior covariance matrix is Ω 3 = τ -1 Σ. Comparing Ω 3 to the matrix Ω 0 associated to the SIR method, it appears that the underlying ideas are opposite. Here, directions corresponding to large variances are most likely. The associated GRSIR estimator of the direction b is the eigenvector of ( Σ2 + τ I p ) -1 ΣΓ associated to its largest eigenvalue. In the following, this estimator will be referred to as the Tikhonov estimator. Indeed, let us recall that SIR estimator is obtained by a spectral decomposition of Σ-1 Γ. For all k = 1, . . . , p, denote by x k the k-th column of this matrix. Computing x k is equivalent to solving with respect to x the linear system Σx = Γk where Γk is the k-th column of Γ. The associated Tikhonov minimization problem (see (1.34) in [START_REF] Vogel | Computational methods for inverse problems[END_REF]) can be written as

x k = arg min x Σx -Γk 2 + τ x 2 = ( Σ2 + τ I p ) -1 ΣΓ k .
Thus, in this framework, ( Σ-1 Γ) k is estimated by ( Σ2 +τ I p ) -1 ΣΓ k and consequently Σ-1 Γ is estimated by ( Σ2 + τ I p ) -1 ΣΓ .

Dimension reduction preprocessing. It has been seen in Proposition 3 that the PCA+SIR approach is equivalent to using the prior covariance matrix Ω 2 in GRSIR. The following result is an extension to more general covariance matrices.

Proposition 4 For all real function ϕ let

Ω(ϕ) = d j=1 ϕ( δj )q j qt j .
Then, the associated GRSIR estimator can be obtained by first projecting the predictors on S d = span(q 1 , . . . , qd ) and second performing GRSIR on the projected predictors with prior covariance matrix

Ω(ϕ) = p j=1 ϕ( δj )q j qt j .
The dimension d plays the role of a "cut-off" parameter, since when computing b, all directions qd+1 , . . . , qp are discarded. Three illustrations of this result can be given:

• Choosing ϕ(t) = 1/(τ t), we obtain Ω(1/(τ Id)) = Ω 2 and Ω(1/(τ Id)) = (τ Σ) -1 = Ω 0 ,
where Id is the identity function. It appears that Proposition 3 is a particular case of Proposition 4. As already discussed, since the choice of Ω 0 as a prior covariance matrix seems not very natural, we thus propose two new choices.

• First, ϕ(t) = 1/τ yields

Ω 4 def = Ω(1/τ ) = 1 τ d j=1 qj qt j ,
and Ω(1/τ ) = I p /τ = Ω 1 . Consequently, this new method consists in applying the ridge approach [START_REF] Zhong | RSIR: Regularized Sliced Inverse Regression for motif discovery[END_REF] on the projected predictors, the interpretation being that, in the subspace S d , all directions share the same prior probability. This method will be referred to as PCA+ridge.

• Second, ϕ(t) = t/τ yields

Ω 5 def = Ω(Id/τ ) = 1 τ d j=1 δj qj qt j ,
and Ω(Id/τ ) = Σ/τ = Ω 3 . This new method consists in applying Tikhonov approach on the projected predictors. In this context, directions of S d carrying a large fraction of the total variance of X are more likely. This method will be referred to as PCA+Tikhonov.

4 Simulation study

Experimental setup

GRSIR methods associated to the prior covariance matrices Ω 0 (SIR), Ω 1 (ridge), Ω 2 (PCA+SIR), Ω 3 (Tikhonov), Ω 4 (PCA+ridge) and Ω 5 (PCA+Tikhonov) are compared on data simulated from a random pair (X, Y ) where X ∈ R p with p = 50 and Y ∈ R. The random vector X is Gaussian, centered, with covariance matrix Σ = Q∆Q t where ∆ is the diagonal matrix containing the eigenvalues of Σ defined by δ j = (p + 1j) θ , j = 1, . . . , p and Q is a matrix drawn from the uniform distribution on the set of orthogonal matrices. Following [20], the orthogonal matrix Q is obtained by performing a QR-decomposition on a p × p matrix whose coefficients are independent standard Gaussian random variables. Note that the condition number of Σ is given by p θ and is thus an increasing function of θ.

The random variable Y is defined by the forward single-index model Y = g(β t X/σ X ) + ξ where σ X is the standard deviation of the projection of X on β i.e. σ X = (β t Σβ) 1/2 and ξ is a centered Gaussian random variable independent of X with standard deviation σ ξ = 0.03. Two link functions are considered: g 1 (t) = sin(πt/2) and g 2 (t) = |t -1/2|. The true index β = 5 -1/2 Q(1, 1, 1, 1, 1, 0, . . . , 0) t is located on the subspace spanned by {q 1 , . . . , q 5 } and representing about 27% of the total variance when θ = 2. In our experiments, N = 100 samples of size n = 100 are considered. More precisely, for all i = 1, . . . , n and r = 1, . . . , N , the pairs (X i , Y

i ) are simulated following the model Both criteria take their values in [0, 1] and rely on the computation of the squared cosine between two axes. We refer to [START_REF] Ferré | Determining the dimension in Sliced Inverse Regression and related methods[END_REF] for an adaptation to multi-index models. PC is a proximity measure between the estimator b and the true direction β: A value close to 0 implies a low proximity ( b is nearly orthogonal to β) whereas a value close to 1 implies a high proximity ( b is approximatively collinear with β). SC is a measure of the estimator stability: A value close to 0 implies a low stability (two realizations of b are approximatively orthogonal) whereas a value close to 1 implies a high stability (two realizations of b are approximatively collinear).

Y (r) i = g(β t X i /σ X ) + ξ r , with g ∈ {g 1 , g 2 }.

Results

In the sequel, the number of slices is fixed to h + 1 = 10.

Influence of the regularization parameter. The two previously defined criteria are computed as functions of the regularization parameter τ . A logarithmic scale is adopted, 150 values of log(τ ) regularly distributed in [-5, 25] were considered. Here, we limit ourselves to θ = 2. Moreover, we choose d = 20 in the PCA+SIR, PCA+ridge and PCA+Tikhonov methods. Results obtained with link functions g 1 and g 2 are respectively displayed on Figures 1,2 and Figures 3,4. It appears that, for all the considered methods, the shape of the link function does not influence much the results. For both link functions, SIR gives very poor results, especially in terms of proximity. Ridge and Tikhonov regularizations can bring a significant improvement provided τ is large enough. PCA+SIR obtains reasonable results compared to SIR, with the advantage of not requiring the selection of τ . The selection of d is addressed in our third experiment. Note that PCA+ridge and PCA+Tikhonov methods are less sensitive to the choice of τ than ridge and Tikhonov methods. Besides, they both outperform PCA+SIR for sufficiently large values of the regularization parameter. Unsurprisingly, the stability of GRSIR methods (ridge, Tikhonov, PCA+ridge and PCA+Tikhonov) increases with the regularization parameter.

For one of the N = 100 simulated datasets, the pairs (β t X i , bt X i ), i = 1, . . . , n are plotted on Figure 5 for the link function g 1 and on Figure 6 for the link function g 2 . Two estimators are considered for b: SIR and GRSIR with PCA+Tikhonov prior. GRSIR estimator is computed with the optimal regularization parameter, i.e. the regularization parameter maximizing PC. The horizontal shape of the SIR cloud of points indicates that the estimated projection bt X is almost independent of the true one β t X. Conversely, the GRSIR points are concentrated along the first diagonal, indicating a high correlation between bt X and β t X.

Influence of the condition number. The robustness with respect to the condition number is investigated by varying θ in {0, 0.1, 0.2, . . . , 3}. For each value within this set, the optimal regularization parameter is selected for each method and the corresponding PC is displayed on Figures 7,8 in case of link function g 1 . Similar results have been observed with g 2 . Clearly, SIR is very sensitive to the ill-conditioning of the covariance matrix. For all the other considered methods, results are getting better while the condition number increases. Note that ridge and Tikhonov methods as well as PCA+ridge and PCA+Tikhonov yield very similar results.

Influence of the "cut-off " dimension. This experiment is dedicated to the illustration of the role of d in PCA+SIR, PCA+ridge and PCA+Tikhonov methods. Here, the condition number is fixed by choosing θ = 2. For each value of d in {0, 1, . . . , p}, the optimal regularization parameter is selected for each method and the corresponding PC is displayed on Figure 9 in case of link function g 1 . Similar results have been observed with g 2 . One can see that the PCA+SIR method is very sensitive to d. Indeed, if d is large, then this approach reduces to SIR, whose accuracy is low for large dimensions. At the opposite, PCA+ridge and PCA+Tikhonov results remain stable as d increases, since these methods get close to ridge and Tikhonov methods respectively.

Retrieval of Mars surface physical properties from hyperspectral images

We propose to compare SIR and GRSIR on a nonlinear inverse problem in remote sensing. Hyperspectral remote sensing is a promising space technology regularly selected by agencies with regard to the observation of planets. It allows the collection for each pixel of a scene, the intensity of light energy reflected from materials as it varies across different wavelengths. Hundreds of spectels in the visible and near infra-red are recorded for each image cell. The analysis of these spectral signatures allows the identification of the physical, chemical or mineralogical properties of the surface that may help to understand the geological history of planets. Our goal is to evaluate the physical properties of surface materials on the planet Mars from hyperspectral images collected by the OMEGA instrument onboard the Mars express spacecraft. Our approach is based on the estimation of the functional relationship G between some physical parameters Y and observed spectra X. For this purpose, a database of synthetic spectra is generated by a physical radiative transfer model and used to estimate G. The high dimension of spectra (p = 184 wavelengths) is reduced with SIR and GRSIR.

Data

We focus on an observation of the south pole of Mars at the end of summer, collected during orbit 61 by the French imaging spectrometer OMEGA onboard the Mars Express Mission. A detailed analysis of this image [START_REF] Douté | South pole of Mars: Nature and composition of the icy terrains from Mars Express OMEGA observations[END_REF] revealed that this portion of Mars mainly contains water ice, carbon dioxide and dust. This has led to the physical modeling of individual spectra with a surface reflectance model. This model allows the generation of n = 12, 000 synthetic spectra with the corresponding parameters that constitute a learning database. Here, we study the terrain unit of strong CO 2 concentration determined by a classification method based on wavelets [START_REF] Schmidt | Wavanglet: An efficient supervised classifier for hyperspectral images[END_REF]. It contains approximately 9,000 spectra to analyze. The most important parameters characterizing the morphology of these spectra are the proportions of CO 2 ice, water ice and dust as well as the grain sizes of water ice and CO 2 ice. In the sequel only one parameter, the grain size of CO 2 ice is presented. A detailed analysis for the four other parameters can be found in [START_REF] Bernard-Michel | Retrieval of Mars surface physical properties from OMEGA hyperspectral images using Regularized Sliced Inverse Regression[END_REF].

Methodology

Two methods are compared in order to reverse the hyperspectral image: SIR and GRSIR with the PCA + ridge prior. Both methods aim at finding a lower dimensional predictor space, sufficient to describe the relationship of interest. In practice, it appears that a unidimensional predictor space gives satisfactory results [START_REF] Bernard-Michel | Retrieval of Mars surface physical properties from OMEGA hyperspectral images using Regularized Sliced Inverse Regression[END_REF] leading to the single-index model G(X) = g(β t X). As a consequence, the estimation of the relationship G reduces to computing one direction β, the univariate function g being estimated by piecewise linear regression, see [START_REF] Bernard-Michel | Retrieval of Mars surface physical properties from OMEGA hyperspectral images using Regularized Sliced Inverse Regression[END_REF] for further details. The "cut-off" dimension is fixed to d = 20 in GRSIR. The regularization parameter, fixed to τ = 10 -4.2 , is chosen to minimize the mean squared error when estimating the grain size of CO 2 ice on the learning database itself. The estimated relationship between reduced spectra β t X and Y , the grain size of CO 2 , is presented on Figure 10. It appears that the grain size of CO 2 can be accurately modeled by a nonlinear function applied to a linear combination of the 184 wavelengths.

Results

The inversion of the image from orbit 61 using GRSIR shows a smooth mapping of the grain size of CO 2 (see Figure 12) making it possible to distinguish some areas with large grain sizes of CO 2 ice on the boundaries and some areas with small values inside the cap. On the opposite, SIR estimates assume a small number of values that seem to be randomly distributed (see Figure 11) and correspond to the extremal values of the parameter in the learning database. These poor results can be explained by the high condition number, about 10 14 , of the empirical covariance matrix. Since no ground data is available, it is quite difficult to quantify the accuracy of GRSIR estimates. However, comparisons with other approaches or with estimates from other hyperspectral images of the same portion of Mars [START_REF] Bernard-Michel | Retrieval of Mars surface physical properties from OMEGA hyperspectral images using Regularized Sliced Inverse Regression[END_REF] give consistent results. GRSIR then appears promising for model inversion in remote sensing.

Concluding remarks

A new framework has been presented to regularize the SIR method. It provides new interpretations of the existing methods as well as the construction of new regularization techniques. Among them, it appears that the PCA+ridge and PCA+Tikhonov methods are interesting alternatives to the existing PCA+SIR and ridge methods. The use of a "cut-off" dimension d serves to limit the sensitivity to the choice of the regularization parameter τ . The choice of this dimension itself seems not crucial since for large value the above methods are close to the ridge and Tikhonov techniques. In our experiments, the choice d ≃ p/2 appears as a good heuristic in most situations. Of course, an automatic selection of (τ, d) would be of interest. To this end, the construction of a generalized cross-validation criterion is under investigation. We also plan to study the introduction on non-Gaussian priors to obtain L 1penalizations and thus sparse estimates of β.

Proof of Proposition 1. Annulling the gradients of G(µ, V, b, c) yields the system of equations

1 2 ∇ µ G = V -1 (μ -X + st ĉ V b) = 0, (11) 1 2
∇ b G = (ĉ t s) 2 + ĉt W ĉ V b -M t ĉ + (ĉ t s)(μ -X) = 0, (12) 1 2
∇ c G = ( bt V b)(ss t + W )ĉ -M b + (μ -X) t bs = 0, (13)
∇ V G = V -1 + bb t (ĉ t s) 2 + ĉt W ĉ -V -1 (μ -X)(μ -X) t + Σ V -1 = 0. ( 14
)
From ( 11), we have μ = Xst ĉ V b.

Replacing in ( 12) and ( 13) yields the simplified system of equations

(ĉ t W ĉ) V b = M t ĉ, (16)
( bt V b)W ĉ = M b. ( 17
)
Assuming W is regular, ( 17) entails ĉ = W -1 M b/ bt V b and replacing in ( 16) yields

(ĉ t W ĉ)( bt V b) V b = M t W -1 M b. ( 18
)
Now, multiplying ( 14) on the left and on the right by V and taking account of ( 15) entail

Σ = V + (ĉ t W ĉ) V bb t V . (19)
As a consequence of ( 19), we have

Σb = 1 + (ĉ t W ĉ)( bt V b) V b, (20)
which means that Σb is proportional to V b. Consequently, one also has

V b = θ( b) Σb , (21)
where we have defined θ( b) = bt V b/ bt Σb . Substituting ( 21) in ( 18) yields

(ĉ t W ĉ)( bt V b)θ( b) Σb = M t W -1 M b
and thus b is an eigenvector of Σ-1 M t W -1 M . Let us denote by λ the associated eigenvalue

λ = (ĉ t W ĉ)( bt V b)θ( b). (22)
Collecting ( 20) and ( 21) yields 1 + (ĉ t W ĉ)( bt V b) = 1/θ( b) and thus the eigenvalue can be rewritten as λ = 1θ( b). Moreover, we have, from [START_REF] Ferré | Smoothed functional inverse regression[END_REF],

ĉt M b = λ/θ( b), (23) tr
( Σ V -1 ) = p + λ/θ( b), (24)
log det V = log det Σlog det

I p + (ĉ t W ĉ) V bb t = log det Σ -log 1 + λ/θ( b) , (25)
entailing

G(μ, V , b, ĉ) = p + log det Σ -log 1 + λ/θ( b) = p + log det Σ + log(1 -λ).
As a consequence, to minimize G, λ should be the largest eigenvalue. Finally, let us consider

V 0 = Σ - λ bt Σb Σbb t Σ, leading to V 0 + (ĉ t W ĉ)V 0 bb t V 0 = Σ + (ĉ t W ĉ)θ 2 ( b) - λ bt Σb Σbb t Σ = Σ + λθ( b) bt V 0 b (ĉ t W ĉ)( bt V 0 b)θ( b) -λ Σbb t Σ = Σ,
in view of ( 22) and thus V 0 verifies equation [START_REF] Gannoun | An asymptotic theory for SIR α method[END_REF].

Proof of Corollary 1. Let us remark that, under [START_REF] Cook | Discussion of "Sliced Inverse Regression for dimension reduction" by K.C. Li[END_REF], the coefficients W ij of W are given by W ij = f i I{i = j}f i f j for all (i, j) ∈ {1, . . . , h} 2 . The inverse matrix of W is

W -1 = diag 1 f 1 , . . . , 1 f h + 1 f h+1 U,
where U is the h × h matrix defined by U ij = 1 for all (i, j) ∈ {1, . . . , h} × {1, . . . , h}. Since the jth line of M is given by f j ( Xj -X) t for all j = 1, . . . , h and taking account of U 2 = f h+1 U , we have

M t W -1 M = h j=1 f j ( Xj -X)( Xj -X) t + 1 f h+1 M t U M = h j=1 f j ( Xj -X)( Xj -X) t + 1 hf h+1 (M t U )(M t U ) t . (26)
Now, remarking that all the columns of M t U are equal to

h j=1 f j ( Xj -X) = h j=1 f j ( Xj -X) -f h+1 ( Xh+1 -X) = -f h+1 ( Xh+1 -X), it follows that (M t U )(M t U ) t = hf 2 h+1 ( Xh+1 -X)( Xh+1 -X) t
and thus replacing in [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF] yields

M t W -1 M = h+1 j=1 f j ( Xj -X)( Xj -X) t = Γ.
The result is then a consequence of Proposition 1.

and thus b is an eigenvector of (Ω Σ + I p ) -1 ΩM t W -1 M . Let us denote by λ the associated eigenvalue

λ = (ĉ t W ĉ)( bt V b)(1 + η( b))θ( b). (33)
Collecting ( 31), ( 32) and ( 33) yields λ = 1θ( b). Now, let us remark that ( 23), ( 24) and ( 25) still hold in this context entailing

G Ω (μ, V , b, ĉ) = p + log det Σ + log(1 -λ).
As a consequence, to minimize G Ω , λ should be the largest eigenvalue. The end of the proof follows the same lines as the one of Proposition 1.

Proof of Proposition 3. Let P be the projection matrix on S d :

P = d j=1 qj qt j ,
and, for all i = 1, . . . , n consider the projected predictor defined by Xi = P X i . Introducing Γ = P ΓP the empirical "between slices" matrix associated to X1 , .

  This yields N samples, with common predictors but different noises, as well as N estimators b(r) , r = 1, . . . , N of β. Two criteria are used to assess the estimates accuracy: The Proximity Criterion (PC) and the Stability Criterion (SC) defined as (s) ) t b(r) ) 2 .

Figure 1 :

 1 Figure 1: Influence of the regularization parameter. The link function is g 1 and the condition number is fixed (θ = 2). Horizontally: log(τ ), vertically: (a) PC and (b) SC. Continuous line: Ω 0 (SIR), -• -: Ω 1 (ridge), dashed line: Ω 3 (Tikhonov).

  SC as a function of τ

Figure 2 :

 2 Figure 2: Influence of the regularization parameter. The link function is g 1 , the "cut-off" dimension and the condition number are fixed (d = 20 and θ = 2). Horizontally: log(τ ), vertically: (a) PC and (b) SC. Continuous line: Ω 2 (PCA+SIR), -• -: Ω 4 (PCA+ridge), dashed line: Ω 5 (PCA+Tikhonov).

Figure 3 :

 3 Figure 3: Influence of the regularization parameter. The link function is g 2 and the condition number is fixed (θ = 2). Horizontally: log(τ ), vertically: (a) PC and (b) SC. Continuous line: Ω 0 (SIR), -• -: Ω 1 (ridge), dashed line: Ω 3 (Tikhonov).

Figure 4 :

 4 Figure 4: Influence of the regularization parameter. The link function is g 2 , the "cut-off" dimension and the condition number are fixed (d = 20 and θ = 2). Horizontally: log(τ ), vertically: (a) PC and (b) SC. Continuous line: Ω 2 (PCA+SIR), -• -: Ω 4 (PCA+ridge), dashed line: Ω 5 (PCA+Tikhonov).

22

 22

Figure 5 :

 5 Figure5: Comparison of pairs (β t X i , bt X i ), i = 1, . . . , n obtained with SIR (•) and GR-SIR, PCA+Tikhonov prior (×). The link function is g 1 and the optimal regularization parameter is used. The "cut-off" dimension and the condition number are fixed (d = 20 and θ = 2).

Figure 6 :Figure 7 :

 67 Figure6: Comparison of pairs (β t X i , bt X i ), i = 1, . . . , n obtained with SIR (•) and GR-SIR, PCA+Tikhonov prior (×). The link function is g 2 and the optimal regularization parameter is used. The "cut-off" dimension and the condition number are fixed (d = 20 and θ = 2).

Figure 8 :Figure 9 :

 89 Figure8: Sensitivity of GRSIR with respect to the condition number of the covariance matrix. The link function is g 1 , the "cut-off" dimension is fixed to d = 20 and the optimal regularization parameter is used for each value of θ. Horizontally: θ, vertically: PC. Continuous line: Ω 2 (PCA+SIR), -• -: Ω 4 (PCA+ridge), dashed line: Ω 5 (PCA+Tikhonov).

Figure 10 :

 10 Figure10: Estimated functional relationship using piecewise linear regression between reduced spectra βt X on the first GRSIR (PCA+ridge prior) direction and Y , the grain size of CO 2 . Horizontally: reduced spectra from the learning database on the first GRSIR direction. Vertically: Grain size of CO 2 .

Figure 11 :

 11 Figure 11: Grain size of CO 2 ice estimated by piecewise linear regression from the hyperspectral image observed on Mars during orbit 61 when the dimension of the predictor space is reduced by SIR.

Figure 12 :

 12 Figure 12: Grain size of CO 2 ice estimated by piecewise linear regression from the hyperspectral image observed on Mars during orbit 61 when the dimension of the predictor space is reduced by GRSIR (PCA+ridge prior).

Acknowledgments

This research is partially supported by IAP research network grant nr. P6/03 of the Belgian government (Belgian Science Policy). We are grateful to Sylvain Douté for his important contribution to the data and to the referees whose remarks led to an important improvement of the presentation.

Appendix: Proofs

Proof of Lemma 1. Let us remark that

where we have defined for i = 1, . . . , n,

Since Z 2,. and Z 3,. are centered, replacing the previous expansion in [START_REF] Cook | Fisher lecture: Dimension reduction in regression[END_REF] yields

where

and the conclusion follows.

Proof of Lemma 2. The joint distribution of (X, b) given Y denoted by p(X, b|Y ) is calculated as the product p(X|Y, b)p(b|Y ) where p(X|Y, b) is given by (4) and p(b|Y ) is given by [START_REF] Cook | Regression graphics. Ideas for studying regressions through graphics[END_REF]. The estimators are obtained by minimizing

Proof of Proposition 2. The proof is similar to the one of Proposition 1. First, remark that equation ( 11) still holds and thus μ = X -s t ĉ V b. Let us recall the following definitions

Annulling the gradients of G Ω (µ, V, b, c) yields the system of equations

Assuming W is regular, equation ( 28) entails ĉ = W -1 M b/((1+η( b))( bt V b)) and replacing in [START_REF] Schmidt | Wavanglet: An efficient supervised classifier for hyperspectral images[END_REF] yields

Now, multiplying (29) on the left and on the right by V , it follows that

which means that Σb is proportional to V b. As a consequence, one also has

and replacing in (


Peggy Cénac
email: peggy.cenac@u-bourgogne.fr

Stéphane Loisel
email: stephane.loisel@univ-lyon1.fr

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

Clémentine Prieur
email: clementine.prieur@imag.fr

P C Énac

Risk indicators with several lines of business: comparison, asymptotic behavior and applications to optimal reserve allocation

come   L'archive ouverte pluridisciplinaire

Introduction

The current change of regulation leads the insurance industry to address new questions regarding solvency. In Europe, insurance groups will soon have to comply with new rules, namely Solvency II. In comparison to the previous regulation system, Solvency II aims at defining solvency requirements that are better adjusted to the underlying risks. Solvency requirements may either be computed thanks to a standard formula, or with internal models that companies are encouraged to develop. While a bottom-up approach is used in the standard formula to aggregate risks (one first studies each small risk separately and then aggregates them thanks to a kind of correlation matrix), a top-down approach may be used in some internal or partial internal models to allocate economic capital: once the main risk drivers for the overall company have been identified and the global solvency capital requirement has been computed, it is necessary to split this overall buffer capital into marginal solvency capitals for each line of business, in order to penalize as equitably as possible lines of business and customers according to the share of the overall risk they represent after diversification. Once this is done, one may also want to allocate some free additional surplus in some different zones in order to avoid as far as possible that some lines of business become insolvent too often. Capital fungibility between lines of business or between entities of a large insurance group that lie in different countries is indeed limited by different entityspecific or country-specific solvency constraints. For example, when AIG experienced problems in the USA, it was impossible for the group to transfer some funds from some European branches to the ones in distress. Some European branches were able to continue business almost as usual and were not impacted, while some other lines would have to get up to 37.8 billion of fresh liquidity.

One possible way to define optimality of the global reserve allocation is to minimize the expected sum of penalties that each line of business would have to pay due to its temporary potential insolvency, following ideas presented in [START_REF] Loisel | Differentiation of some functionals of risk processes and optimal reserve allocation[END_REF]. If one neglects discounting factors, a first approximation of this penalty is given by the time-integrated expected negative part of the surplus process. Closed-form formulas were available in the classical risk model for exponentially distributed claim amounts, which led to a semiexplicit optimal reserve allocation. However, this approach does not take into account the dependence structure. Also, one could object that if the company is ruined at the group level (i.e. the sum of the surpluses is negative), then the allocation does not change anything. Following [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF], we define and study the expected orange area: the reserve of line k at time j is denoted by R k j (defined as u k plus the aggregate premium minus the aggregate claim of the k th branch over the first j periods). We look for the optimal allocation (u 1 , . . . , u d ) -where u k stands for the part of capital allocated to the k th branch -that minimizes

I 1 = p∧τ j=1 d k=1 E |R k j |1 1 {R k j <0} 1 1 {R 1 j +•••+R d j ≥0}
under the constraint that

u 1 + • • • + u d = u, where τ = inf j ∈ N * , R 1 j + • • • + R d j < 0 .
The I 1 indicator is the stopped version of the indicator considered in [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF]:

I 2 = p j=1 d k=1 E |R k j |1 1 {R k j <0} 1 1 {R 1 j +•••+R d j ≥0} .
Remark that in the case p = 1, I 1 = I 2 . Below, I 1 is referred to as the stopped orange area and I 2 is referred to as the orange area. In the present paper, we also consider the alternative risk indicator, called expected violet area:

J = p j=1 d k=1 E |R k j |1 1 {R k j <0} 1 1 {R 1 j +•••+R d j ≤0} .
The three studied cases are represented in the figure below. In the case when p = 1, these three indicators may be considered as particular cases of Dhaene et al.'s [START_REF] Dhaene | Optimal Capital Allocation Principles[END_REF] general framework which resumes as the minimization of

d j=1 E ξ j D X j -u j v j
where X j stands for the loss of the j th branch, (v j ) for a weight sequence and D is a deviation function.

Our proposal corresponds to take ξ j = 1 1 {S≤u} or ξ j = 1

Local insolvency

Global insolvency u1 R 1 1 u2 R 2 1 R 1 1 + R 2 1 u = u1 + u2 d k=1 R k p R 1 p R 2 p

Stopped orange area

Figure 1. Representation of the three indicators in [START_REF] Dhaene | Optimal Capital Allocation Principles[END_REF]).

The number p of periods could also correspond to one year in the Solvency II Solvency Capital Requirement computation problem, or to 3 to 5 years in the Own Risk and Solvency Assessment (ORSA) framework. ORSA is defined as the process undertaken by an insurer or an insurance group to assess the adequacy of its risk management and current and future solvency position; in that setting, time is often discrete to limit the complexity and one is interested in ruin at yearly inventories. The Article 45 of the Solvency II directive defines the ORSA as a part of the risk-management system of an insurance or reinsurance undertaking to conduct its own risk and solvency assessment (see [START_REF]EIOPA Consultation Paper On the Proposal for Guidelines on Own Risk and Solvency Assessment[END_REF] for example). The National (U.S) Association of Insurance Commissioner (NAIC) had also initiated reflections and proposed principles on ORSA, a review may be found in [START_REF]Risk Metrics for Decision Making and ORSA, Edited by the SOA[END_REF].

Here we consider very basic insurance discrete time models for different lines, and we try to study the optimal reserve allocation. We shall see that considering the orange area (I 1 or I 2 indicator) is a completely different approach from the economic capital allocation problem, where one tries to penalize contributions to the overall risk and to diversification effect.

As noticed in [START_REF] Dhaene | Optimal Capital Allocation Principles[END_REF], different capital allocations must in some sense correspond to different questions that can be asked within the context of risk management. A review on various ways to allocate capital may be found in [START_REF] Decupère | Agrégation des risques et allocation de capital sous Solvabilité II, Mémoire d'actuariat[END_REF]. Also J.-P. Laurent [START_REF] Laurent | Analytics of risk allocations within financial institutions and at a systemic level[END_REF] gives a global view of capital allocation and exhibits connections between several theories (including Euler, Shapley-Aumann and Pareto optimal allocations). Let us remark that even if our purpose implies a risk driven allocation of capital, we do not refer explicitly to risk measure theory (see [START_REF] Artzner | Heath Coherent Measures of Risk[END_REF][START_REF] Goovaerts | Decision principles derived from risk measures[END_REF][START_REF] Goovaerts | Actuarial risk measures for financial derivative pricing[END_REF]). An example is the TVaR allocation principle (which appears also as a particular case of [START_REF] Dhaene | Optimal Capital Allocation Principles[END_REF], see equation (54)), according to which each line of business k has to pay for cost of capital corresponding to average marginal loss for line k given that things go wrong for the global company (the VaR of the sum is overshot). This allocation principle is similar in spirit to the minimization of the indicator J (violet area), while the orange area is placed in a different context and relevant when the global company is safe. We are more interested in the stability and liquidity issues with limited capital fungibility. When the global company is ruined, this allocation has no importance (even if in practice, one may break the mutual solidarity principle and decide to run separately healthy lines). Finally, note that we do not take into account the penalties in the overall reserve process; in theory the line of business or the global company would have to pay the penalty upfront, and this could accelerate ruin.

The present article is an attempt to explore the quantitative behavior of the allocation procedure, based on the minimization of a risk indicator. While Dhaene et al. (see Theorem 3 in [START_REF] Dhaene | Optimal Capital Allocation Principles[END_REF]) show that the optimization problem leads to compute quantile functions and comonotonic sums related to each risk and to the specific form of the indicator (see Theorem 3), we propose explicit computations in dimension 2 and provide a simulation study in higher dimension and for multi-periodic models. The simulation study is made possible thanks to the algorithmic procedure proposed in [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF]. Let us emphasize that in [START_REF] Dhaene | Optimal Capital Allocation Principles[END_REF], the question of the effective computation of the solution to the optimization problem is not addressed.

Our paper is organized as follows. In Section 1, we tackle the case when there are two lines of business (d = 2) and one period (p = 1). We derive semiexplicit formulas and/or asymptotics in the cases where aggregate claim amounts are independent exponentials, independent Generalized Pareto Distribution (GPD), correlated GPD, for the minimization of the orange area.

In Section 2, we present the same computations for the minimization of the violet area, and compare both minimization strategies. Section 3 consists in a simulation study for the case treated in Section 1 and Section 2, using the algorithm developed in [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF]. This section helps for benchmarking the parameters of the algorithm. In Section 4, models in higher dimension (d = 10), still for one period, are considered. In Section 5, we investigate time horizons p > 1 and d = 2, . . . , 5 lines of business. For the models considered in these two last sections, semi-explicit formulas are generally not available.

Time horizon p = 1, d = 2 lines

In this section, we consider the one-period problem with two lines of business. Let

R k 1 = u k + c k -X k
, where for line k = 1, 2, u k is the initial reserve, c k is the premium and X k is the aggregate claim amount during the first period. We start with a basic example with a discrete distribution. It illustrates the fact that the use of the orange area may lead to inconsistent allocation, for example, if the total capital u is too small. Then we consider more meaningful models for the claims, for which we may derive explicit solutions.

1.1. A basic example. The following example is very basic but a caricature. It was chosen to highlight the fact that the expected orange area must not be used to penalize lines of business according to their risk contributions. Imagine that c 1 = c 2 = 10, that X 1 = 1000Y 1 and X 2 = 20Y 2 , where Y 1 and Y 2 are independent Bernoulli random variables with parameter 1/2, and that u = 10. The expected orange area is equal to 10-u 2 4 for 0 ≤ u 2 ≤ u, because the aggregate reserves are negative if X 1 = 1000. Consequently, the optimal reserve allocation would be to allocate all the capital to the less risky line of business (line 2). This decision corresponds to minimizing penalties in case of temporary insolvency of one line, but is in complete disagreement with any reasonable economic capital allocation principle. Let us remark that in this example, the violet area equals 1 4 (1980 -u 1 ), so that it is minimized by allocating all the capital to the most risky branch, as expected. We remark that in this example, u and the premiums c 1 and c 2 are not consistent with the claims X 1 and X 2 . Below a , we have taken u = 960, c 1 = 50 and c 2 = 1. We have considered two cases: firstly Y 1 and Y 2 are independent Bernoulli variables with parameter 1 2 . Secondly, they are still independent Bernoulli variables but Y 1 has parameter 1 2 and Y 2 has parameter

1 3 . 1.1.1. Y 1 ; B( 1 2 ), Y 2 ; B( 1 2
). A direct computation gives:

I 1 = I 2 = J = 1 4 (950 -u 1 )1 1 {u 1 <950} + (19 -u 2 )1 1 {u 2 <19} .
The violet and the orange area are the same. Under the constraint u 1 +u 2 = u, it is minimized for 10 ≤ u 2 ≤ 19 (in this case, there is not a unique minimum). [START_REF] Artzner | Heath Coherent Measures of Risk[END_REF] ). A direct computation gives:

1.1.2. Y 1 ; B( 1 2 ), Y 2 ; B( 1
I 1 = I 2 = 1 6 2 * (950 -u 1 )1 1 {u 1 <950} + (19 -u 2 )1 1 {u 2 <19} , J = 1 6 (950 -u 1 )1 1 {u 1 <950}} + (19 -u 2 )1 1 {u 2 <19} .
There is a unique minimum for I 1 = I 2 , which is reached for u 2 = 10 and u 1 = 950.

In the case of the violet area J, the minimum is reached for u 2 ∈ [10; 19], the minimum is not unique.

Let us now assume that the insurance company has previously correctly penalized lines of business according to their risks, and then wants to allocate some additional safety capital u in order to avoid temporary insolvencies of some lines. The orange area has been designed to answer this question.

1.2. General remarks. We now consider more general cases, still for the period p = 1 and focus on the orange area. For the sake of simplicity, in the remainder of this section, we will consider u k instead of u k + c k . The premium is thus included in the allocation here. Remark that for p > 1 this would be a non sense to include the premium in the allocation.

a The premiums have been taken equal to 10% of the expected value of the claims.

The gradient with respect to (u 1 , u 2 ) of the average orange area I 2 is given by

P(R 1 1 < 0, R 1 1 + R 2 1 ≥ 0), P(R 2 1 < 0, R 1 1 + R 2 1 ≥ 0)
. We are looking for non degenerate solutions (i.e. each line receives a positive initial reserve u k in the optimal strategy), an easy calculation with Lagrange's multipliers with the constraint u 1 + u 2 = u leads to an optimal allocation when

P R 1 1 < 0, R 1 1 + R 2 1 ≥ 0 = P R 2 1 < 0, R 1 1 + R 2 1 ≥ 0 , which can be rewritten as (1.1) P X 1 > u 1 , X 1 + X 2 ≤ u = P X 2 > u 2 , X 1 + X 2 ≤ u .
We shall consider several probabilistic models for the claims X 1 and X 2 : independent exponentials (Section 1.3), independent GPD (Section 1.4), as well as conditionally independent exponentials (Section 1.5) which are correlated Pareto distributed (GPD).

1.3. Independent exponentials. Assume X 1 and X 2 are independent exponential random variables with respective parameters µ 1 and µ 2 . In the particular case where µ 1 = µ 2 , the optimal allocation is

u 1 = u 2 = u/2.
From now on, assume that µ 2 > µ 1 .

Optimal allocation.

Let us denote α and β the real numbers such that µ 2 = αµ 1 and u 1 = βu, α > 1 and 0 ≤ β ≤ 1. Equation (1.1) leads to the optimal allocation

(1.2) (α -1) (h (β) -h (α(1 -β))) -αh(1) + (1 + α)h (β + α -αβ) = h(α)
where h is the function defined by h(x) = exp -uµ 1 x .

Proposition 1.1. There exists a unique β = ψ(α, u, µ 1 ) ∈ [0, 1] satisfying Equation (1.2). Moreover, 0 < ψ(α, u, µ 1 ) ≤ α α+1 and α → ψ(α, u, µ 1 ) is non decreasing.

Proof. The proof is a straightforward application of the convexity of h and the implicit function theorem.

Remark 1.2. The fact that ψ(α, u, µ 1 ) (and thus u 1 ) increases with α is consistent with the fact that if α increases then the first branch becomes riskier.

1.3.2. Asymptotic as the capital u goes to infinity. Studying the asymptotic as u → ∞ is useful in order to make comparison between models easier.

Proposition 1.3. For α > 1, lim u→∞ ψ(α, u, µ 1 ) = α α + 1 . Proof. Multiplying (1.2) by h(-β), we get (α-1)-αh(1-β)+(1+α)h (α(1 -β))-(α-1)h (α(1 -β) -β)-h(α-β) = 0.
Remark that lim u→∞ ψ(α, u, µ 1 ) = 1 since the left hand side of the previous equation has to go to 0 when u tends to ∞. Then, one has

lim u→∞ -αh 1 -ψ(α, u, µ 1 ) +(1+α)h α 1 -ψ(α, u, µ 1 ) -h α -ψ(α, u, µ 1 ) = 0,
and consequently h α 1 -ψ(α, u, µ 1 ) -ψ(α, u, µ 1 ) -→

f X 1 (x) = a b 1 + x b -a-1 , f X 2 (x) = αa b 1 + αx b -a-1
and

F X 1 (x) = 1 + x b -a , F X 2 (x) = 1 + αx b -a
.

Equation (1.1) yields in this case

F X 1 (u 1 ) - a b u u 1 1 + x b -a-1 1 + α(u -x) b -a dx = F X 2 (u 2 ) - a b u u 2 1 + αx b -a-1 1 + (u -x) b -a dx. (1.3)
We do not know how to derive an explicit optimal allocation from this equation. Nevertheless, an asymptotic as u → ∞ is reachable. We begin with a general result for independent variables. Theorem 1.4. Let X 1 and X 2 be two independent, continuous and non negative random variables such that b F X 1 (x) = Θ(F X 2 (x)) as x → ∞. The solution (u 1 , u 2 ) to Equation (1.1), under the constraint

u 1 + u 2 = u is equivalent, as u → ∞ to F X 1 (u 1 ) = F X 2 (u 2 ).
Proof. First, we start with the proof of the existence of two real numbers κ 1 and κ 2 such that 0 < κ 1 < κ 2 < 1 and for u large enough, (1.4)

κ 1 ≤ u 1 u ≤ κ 2 . Let β(u) = u 1 u . If (1.4) was not satisfied then taking if necessary a sequence u k → ∞, we could assume that either β → 0 or β → 1. If β → 0 then F X 1 (u 1 ) → 0 and F X 2 (u 2 ) → 1. Remark that P(X 1 > u 1 , X 1 + X 2 ≤ u) = F X 1 (u 1 ) -P(X 1 > u 1 , X 1 + X 2 > u),
so that (1.1) rewrites:

F X 1 (u 1 )-P(X 1 > u 1 , X 1 +X 2 > u) = F X 2 (u 2 )-P(X 2 > u 2 , X 1 +X 2 > u).
Both terms P(X 1 > u 1 , X 1 + X 2 > u) and P(X 2 > u 2 , X 1 + X 2 > u) go to 0 as u → ∞, so that we can not have F X 1 (u 1 ) → 0 and F X 2 (u 2 ) → 1.

b Recall that f (x) is said to be Θ(g(x)) as x → ∞ if the ratio f (x) g(x) is bounded (from below and from above) as x → ∞. and thus β can not go to 0. The same reasoning proves that β cannot go to 1 and thus (1.4) is satisfied. Now,

P(X 1 > u 1 , X 1 + X 2 > u) = P(X 1 > u 1 , X 2 > √ u, X 1 + X 2 > u) (1)
+ P(X 1 > u 1 , X 2 < √ u, X 1 + X 2 > u) (2)
.

The term (1) is less than

P(X 1 > u 1 )P(X 2 > √ u) = o(F X 1 (u 1 )). The term (2) is less than F X 2 ( √ u)F X 1 (u - √ u)
. Thus, (1.4) leads to

F X 1 (u - √ u) ≤ F X 1 (u 1 )
for u large enough, so that

P(X 1 > u 1 , X 1 + X 2 > u) = o F X 1 (u 1 )
and

P(X 2 > u 2 , X 1 + X 2 > u) = o F X 2 (u 2 ) .
The fact that F X 1 (x) = Θ(F X 2 (x)) yields the desired result.

Proposition 1.5. The unique solution to (1.3) satisfies

lim u→∞ β(u) = α α + 1 ,
where β(u) = u 1 u . Proof. The proof is a straightforward consequence of Theorem 1.4 together with the resolution of F X 1 (u 1 ) = F X 2 (u 2 ). 1.5. Correlated Pareto. In this section, we consider a model inspired from [START_REF] Albrecher | Explicit ruin formulas for models with dependence among risks[END_REF]. Let X be a random variable with exponential distribution E(Θ), where the parameter Θ is Γ(a, b) distributed. Recall that the resulting mixed survival function of X is given by

1 -F X (x) = ∞ 0 e -θx f Θ (θ)dθ = 1 + x b -a , x ≥ 0.
In other words, X has a Generalized Pareto distribution (GPD).

1.5.1. Explicit solution for fixed u. Let X 1 and X 2 be two GPD according to the above model. Conditionally to µ 1 , we assume that X 1 ∼ E(µ 1 ), X 2 ∼ E(α µ 1 ) with 1 < α, X 1 and X 2 are independent. We assume moreover that µ 1 ∼ Γ(a, b), a, b > 0. Following [START_REF] Albrecher | Explicit ruin formulas for models with dependence among risks[END_REF], we have that the dependence structure of the X i 's is given by a survival Clayton copula. Conditionally to µ 1 equation (1.1) leads to:

(α-1)h(β)-αh(1)+h(α-αβ+β) = (α-1)h(α(1-β))-αh(α-αβ+β)+h(α).

We now integrate with respect to µ 1 and get the same equation replacing the function h by s(x) = (1 + x u b ) -a . Hence this equation can be rewritten as f (α, β, u) = 0 with

f (α, β, u) = (α -1) (s(β) -s (α(1 -β)))
-αs(1) + (α + 1)s (α -αβ + β) -s(α). (1.5) Proposition 1.6. There exists a unique β = Φ(α, u, a, b) ∈ [0, 1] satisfying equation (1.5). Moreover, 0 < Φ(α, u, a, b) ≤ α α+1 . Proof. We remark that f (α, 0, u) = -f (α, 1, u) = (α -1) + s(α) -αs(1) ≥ 0 as by convexity of s, one has s(α) -s(1) α -1 ≥ s(1) -s(0) .

Moreover one has

b a(α -1)u ∂f ∂β = -g(β) + (α + 1)g(α -αβ + β) -αg(α(1 -β)) = -g(β) + g(α -αβ + β) + αg(α -αβ + β) -αg(α(1 -β)) (1.6) ≤ 0 with g(x) = (1 + x u b ) -(a+1)
, as g is decreasing. Hence there exists a unique β 0 ∈ [0, 1] such that f (α, β 0 , u) = 0. We also have

(1.7) f (α, α α + 1 , u) = s(γ) -s(α) + αs(γ) -αs(1)
with γ = 2α α+1 . We remark that 1 < γ < α, hence using the convexity of s we have s(γ) -s(α) ≤ (s(1) -s(γ)) α -γ γ -1 which implies that f (α, α α+1 , u) ≤ 0 hence β 0 ≤ α α+1 . 1.5.2. Asymptotic as the capital u goes to infinity. One gets the following result, considering the asymptotic as u goes to ∞ : Proposition 1.7. For α > 1, lim u→∞ Φ(α, u, µ 1 ) = Φ 0 is the solution to:

(1.8) (α-1) Φ -a 0 -(α(1 -Φ 0 )) -a -α+(α+1)(α-αΦ 0 +Φ 0 ) -a -α -a = 0. We have Φ 0 < α α+1 . Proof. Equation (1.5) is equivalent, when u → ∞ to (1.9) (α -1) β -a -(α(1 -β)) -a -α + (α + 1) (α -αβ + β) -a -α -a = 0.
One wants to prove that the solution of (1.9) above satisfies β < α α+1 when α > 1. Indeed, define

g(β) = (α -1)β -a -α + (α + 1)(α -αβ + β) -a -(α -1)α -a (1 -β) -a -α -a .
One has g(0 + ) = +∞ and g(1 -) = -∞. Straightforward computations prove that g is decreasing in β and that g( α α+1 ) < 0. This proves that the solution Φ 0 of (1.9) satisfies Φ 0 < α α+1 .

Comparison with an indicator designed for crisis

Recall that the violet area is given by the indicator

J(u 1 , . . . , u d ) = p j=1 d k=1 E |R k j |1 1 {R k j <0} 1 1 {R 1 j +•••+R d j ≤0}
under the constraint

u 1 + • • • + u d = u.
As mentioned in the introduction, this corresponds to a more classical economic capital allocation principle, closely related, in spirit, to the TVaR based allocation.

In this section, we consider the capital allocation that minimizes this violet area. As the computations are similar to the ones corresponding to the orange area, we do not go into the details. Firstly, we focus on the independent exponential distributions case and on the conditionally exponential independent distributions one. Then, we give asymptotic results for the independent GPD case. For this last purpose, we study the asymptotic behavior of probabilities for subexponential distributions. An easy calculation leads to the convexity of the function J since

(u 1 , . . . , u d ) → R k j 1 1 {R k j <0} 1 1 { d k=1 R k j ≤0}
is a convex function on the set

{(v 1 , . . . , v d ) ∈ (R + ) d , v 1 + . . . + v d = u}
and the sum and expectation of convex functions is convex. Moreover, since x → |x| is strictly convex so does the function J. Thus there exists a unique minimum u * to the optimal allocation problem with the indicator J. It is reached for

u 1 + u 2 = u with (2.1) P(X 1 ≥ u 1 , X 1 + X 2 ≥ u) = P(X 2 ≥ u 2 , X 1 + X 2 ≥ u).
We have the following implicit solutions to the optimal allocation issue (recall that u 1 = βu):

• Independent exponential distributions:

(2.2) (α + 1)e -uµ 1 (α-αβ+β) -αe -uµ 1 -e -αuµ 1 = 0.

• Conditionally exponential independent distributions:

(2.3) (α + 1) 1 + u(α -αβ + β) b -a -α 1 + u b -a -1 + uα b -a = 0.
Going through the limit u → ∞ in equation (2.2) leads to β → 1 while in equation ( 2.3) it leads to

α + α -a -(α + 1)(α -αβ + β) -a = 0.
These asymptotic behaviors are illustrated on Figure 2 below. The bold lines are for the independent exponential model, with parameter µ 1 = 1 20 , the simple lines are for the conditionally independent model with parameter a = 1 and b = 20, in both cases, α = 5. The violet lines are for the J indicator and the orange ones for the I 2 indicator.

2.1. Some general results for independent and subexponential distributions. In this section, we give some results for the asymptotic behavior of the allocation with respect to the violet area, in the case where the two distributions are subexponential and independent. For results on subexponential distributions, we refer to Asmussen [START_REF] Asmussen | Ruin probabilities[END_REF]. Recall that a distribution is subexponential if it is concentrated on [0, ∞[ and its distribution function F satisfies:

(2.4) F * 2 (x) F (x) x→∞ -→ 2,
where F * 2 is the convolution square, that is, the distribution function of two independent variables with distribution function F . The class of subexponential distributions is denoted by S. In [START_REF] Asmussen | Ruin probabilities[END_REF] (Chapter IX, Proposition 1.4), We derive the following result from the properties of the subexponential class.

Theorem 2.2. Let X and Y be two independent random variables, concentrated on R + . Let F X be the distribution function of

X, u, v ∈ R+, v ≤ u. Assume that (1) v u → β with 0 < β < 1, (2) F X ∈ S, (3) for any t > 0, F X (tu) F X (u) u→∞ = O(1). Then (2.5) lim u→∞ P(X ≥ v, X + Y ≥ u) F X (u) = 1.
Proof. We closely follow ideas from [START_REF] Asmussen | Ruin probabilities[END_REF], Chapter IX. The fact that lim inf

u→∞ P(X ≥ v, X + Y ≥ u) F X (u) ≥ 1 follows from P(X ≥ v, X + Y ≥ u) ≥ P(X ≥ v, max(X, Y ) ≥ u)
and

P(X ≥ v, max(X, Y ) ≥ u) = P(X ≥ u) + P(X ≥ v, Y ≥ u) -P(X ≥ u, Y ≥ u) = F X (u) + F X (v)F Y (u) -F X (u)F Y (u).
Now, we shall prove that lim sup

u→∞ P(X ≥ v, X + Y ≥ u) F X (u) ≤ 1.
Let w ∈ R + be fixed. We have

P(X ≥ v, X+Y ≥ u) = P(X ≥ v, Y ≤ w, X+Y ≥ u)+P(X ≥ v, Y ≥ w, X+Y ≥ u).
For u large enough,

P(X ≥ v, Y ≤ w, X + Y ≥ u) = w 0 F Y (dy)F X (max(v, u -y)) = w 0 F Y (dy)F X (u -y).
Now, Proposition 2.1 gives

lim u→∞ P(X ≥ v, Y ≤ w, X + Y ≥ u) F X (u) = F Y (w).
On the other hand,

P(X ≥ v, Y ≥ w, X + Y ≥ u) = u-v w F Y (dy)F X (u -y) + ∞ u-v F Y (dy)F X (v) ≤ F X (u -w)F Y (w) + F X (v)F Y (u -v). So that lim sup u→∞ P(X ≥ v, Y ≥ w, X + Y ≥ u) F X (u) ≤ F Y (w).
We deduce the following result for independent GPD.

Corollary 2.3. If X 1 and X 2 are like in Section 1.4, then the unique solution to (2.1) satisfies

β(u) u→∞ -→ 1 where β(u) = u 1 u .
Proof. We apply Theorem 2.2 to X 1 and X 2 independent GPD with respective parameters ( 1 a , b a ) and ( 1 a , b αa ), a, b > 0, α > 1. We consider the quantities

ℓ 1 = P(X 1 ≥ u 1 , X 1 + X 2 ≥ u) P(X 1 ≥ u 1 ) and ℓ 2 = P(X 2 ≥ u 2 , X 1 + X 2 ≥ u) P(X 1 ≥ u 1 ) . If u 1 + u 2 = u and (2.1) is satisfied then ℓ 1 = ℓ 2 . If β(u) = u 1
u is bounded away from 0 and 1, then up to considering a sequence of real numbers u going to infinity, we may assume that

u 1 u → β with 0 < β < 1.
In that case, applying Theorem 2.2 gives

lim u→∞ ℓ 1 = β a and lim u→∞ ℓ 2 = β α a which is contradictory with the fact that ℓ 1 = ℓ 2 .
So that, up to considering a sequence of real numbers u going to infinity, we have that either u 1 u → 0 or u 1 u → 1. Let us assume that u 1 u → 0. Following the lines of the proof of Theorem 2.2, we prove that

lim u→∞ P(X 1 ≥ u 1 , X 1 + X 2 ≥ u) F X 1 (u) = 1 + 1 α a and lim u→∞ P(X 2 ≥ u 2 , X 1 + X 2 ≥ u) F X 1 (u) = 1
α a . This is in contradiction with the fact that the two above expressions are equal. We conclude that β(u) → 1 as u → ∞.

Summary of comparisons.

Below is summarized the asymptotic behavior for the three models: (M1) independent exponential, (M2) independent GPD, (M3) conditionally independent exponential and the two indicators (I) orange area, (J) violet area. The asymptotic of β is expressed with respect to α and to the parameters of the model.

M1 M2 M3 I α α + 1 α α + 1 (α -1)β -a -α + (α + 1)(α -αβ + β) -a -(α - 1)(α(1 -β)) -a -α -a = 0 J 1 1 α + α -a -(α + 1)(α -αβ + β) -a = 0
The algorithm given in [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF] can be easily adapted to the estimation of allocation minimizing the violet area. The keystone for the proof of convergence of this algorithm relies in the strong convexity of the auxiliary function used in the mirror algorithm, in the properties of the Legendre transform and in martingales techniques. The violet area as well as the stopped orange area satisfy the hypothesis of convergence of the algorithm (in [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF], it has been verified for the orange area). The stochastic algorithm is developed in order to estimate the optimal solution (u 1 , . . . , u d ) under the constraint

u 1 + • • • + u d = u.
This algorithm is distribution free and is efficient even when d is quite large (d = 30 or d = 50). This is a Kiefer-Wolfowitz version of the stochastic mirror algorithm.

2.3. Rapid description of the algorithm. The procedure may be summarized in the following picture:

0 1 2 ζ n-1 ζ n α Slope: 1/γ n Gradient descent x n-1 x n W β
The algorithm constructs two random sequences:

⊲ (χ n ) in C = {x ∈ R d , x i ≥ 0, x 1 + • • • + x d = u}, ⊲ (ξ n ) in the dual space E * of R d Algorithm 2.4.
⊲ Initialization:

ξ 0 = 0 ∈ E * , χ 0 ∈ C ⊲ Update: for n = 1 to N do • ξ n = ξ n-1 -γ n ψ(χ n-1 ) • χ n = ∇W βn (ξ n ) ⊲ Output: S N = N n=1 γ n χ n-1 N n=1 γ n
where ψ is the discrete gradient of the vector I ∈ R d whose kth coordinate is

I k =   p j=1 g k (R k j )1 1 {R k j <0} 1 1 { d i=1 R i j ≥0}
  and respectively of the vector J ∈ R d whose kth coordinate is

J k =   p j=1 g k (R k j )1 1 {R k j <0} 1 1 { d i=1 R i j ≤0}   ,
W stands for an auxiliary function whose gradient maps E * into C (see [START_REF] Cénac | Some multivariate risk indicators; minimization by using a Kiefer-Wolfowitz approach to the mirror stochastic algorithm[END_REF] for details as well as the proof of the consistency of the estimator S N of the argmin of I 1 , I 2 or J).

Simulations for

p = 1, d = 2
We shall perform simulations for p = 1 and d = 2 for the models described above (exponential independent = IE, conditionally exponential= CIE, independent GPD = IGPD). Recall that the IGPD and CIE models have the same margins. For the IE and CIE models, we compare the simulated results with the theoretical ones. The simulations are done for the non crisis indicator (orange area) as well as for the crisis indicator (violet area). Then, we have performed simulations in higher dimension (d > 2). for the IGPD model. In order to get the estimation of the minimum, we have performed 10 times the stochastic algorithm on data of length 15 000. The mean and the standard deviation over the 10 estimations are given below. For α = 5 and α = 10, we have taken u = 50, we compare with the theoretical value using the relative mean squared error (rmse) and the mean squared error (mse) for the IE and CIE models. Recall that if u i k is the estimated value of u i on the k th sample, then

rmse(u i ) = 1 k k j=1 u i k -u i u i 2 , and
mse(u i ) = 1 k k j=1 u i k -u i 2 .
3. This table shows that, for the three models, the more the first branch is risky, the more we allocate capital to it. Also, we recall that, asymptotically, the capital allocated to the riskiest branch for the IGPD model is higher than for the CIE model (see Figure 2) . In the table above, the values of u are (relatively) small, we conclude that the asymptotic has not been reached. The same observations as in dimension 2 can be done. Moreover unsurprisingly, the allocation is identical for the branches 2 to 7 and for the branches 8 to 10.

4.2.

A block of correlated GPD and a block of independent GDP, in dimension 10. In this section, we consider a block of 5 conditionally exponential variables and an independent block of 5 independent GPD variables, with the same margins as the 5 correlated GPD variables of the first block. We refer this model to the mixed model (MM).

• Conditional exponential: X 1 1 ; E(Θ) and for i = 2, . . . , 5,

X i 1 ; E(α 1 Θ), where Θ ; Γ(α, b). • Independent GPD: X 1 1 ; GP D( 1 a , b a
) and for i = 6, . . . , 10,

X i 1 ; GP D( 1 a , b α 1 a ).
The MM model has to be compared with the IGPD and CIE models with the same margins. We have chosen α 1 = 5 and u = 80. We have performed our stochastic algorithm 10 times on data sets of length 20 000. This example illustrates the fact that for the expected orange area, lines correlated to the most dangerous one receive less capital that lines with the same marginal distributions but independent from the riskiest line. This is in accordance with the basic example presented in the introduction. If line 1 is insolvent, then it is likely that the group is ruined (the aggregate reserves are likely to be negative). Consequently, a greater part of the insolvencies of lines 2 to 5 than the ones of lines 6 to 10 will not contribute to the orange area because the sum of reserves is negative. Once again, the orange area should not be used for economic capital allocation, as it encourages to take additional risks that are correlated to the main source of risk. Logically, the effect is reversed if one uses the expected violet area. In that case, a greater part of insolvencies of lines 6 to 10 will not contribute, and consequently the required capital for lines 2 to 5 will be higher than for lines 6 to 10, as in the TVaR allocation principle. Below are the results for the same models as above, in the case of the violet area (indicator J). We have performed 10 simulations of length 23 000. 5.2. The stopped orange area, the I 1 indicator. If there is only one period (p = 1), there is no difference between the I 2 and the I 1 indicators. In order to compare the optimization with the I 1 indicator with respect to the I 2 , we have performed simulation with p = 5 and d = 5 on the same models as above (independent exponential, conditionally independent exponential, independent GPD, mixed GPD). We remark that the only notable difference is observed for the mixed model.

5.3.

The violet area, the J indicator. As already noticed, for the violet area, more data is needed because for less of the aggregate loss go over u.

The simulation below are done on data of length 16 000, 10 times.

u→∞ 1 , which leads to lim u→∞ α 1 - 1 . 4 .

 1114 ψ(α, u, µ 1 ) -ψ(α, u, µ 1 ) = 0 and Proposition 1.3 follows. Independent Pareto. We consider two independent Generalized Pareto Distribution (GPD) denoted by X 1 and X 2 with respective parameters 1 a , b a and 1 a , b αa , a, b > 0, α ≥ 1. The density and survival functions are given by:

Figure 2 .

 2 Figure 2. Asymptotic behavior for the independent exponential and conditionally independent exponential models

3. 1 .

 1 A comparison of our three models. We have chosen µ 1 = 1 20 for the IE model, a = 3 and b = 60 for the CIE model and ξ = 1 a , σ = b a

  2. The indicator I 2 : the orange area.

						α = 5
			IE model			CIE model			IGPD model
		u 1	u 2	β	u 1	u 2	β	u 1	u 2	β
	mean	38.37	11.63 0.767	36.8	13.2	0.735 35.83 14.17 0.717
	sd dev	0.085	0.085	. . .	0.115	0.115	. . .	0.133 0.133	. . .
	th. √ rmse 0.0032 0.0105 38.46 11.54 0.769 . . . √ mse 0.121 0.121 . . .	36.84 0.0032 0.009 13.16 0.737 . . . 0.119 0.119 . . .	non available non available non available
						α = 10
	mean	42.64	7.36	0.852	40.95	9.05	0.819 40.05 9.48 0.801
	sd dev	0.161	0.161	. . .	0.138	0.138	. . .	0.164 0.164	. . .
	th. √ rmse 0.0085 0.0516 42.96 7.04 √ mse 0.364 0.363	0.859 . . . . . .	41.22 0.0073 0.0342 8.78 0.301 0.3	0.824 . . . . . .	non available non available non available

Acknowledgments. We are grateful to Dominik Kortschak for fruitful discussions on subexponential distributions. This work has been partially supported by the BNP Paribas Cardif Chair "Management de la modélisation". The views expressed in this document are the authors own and do not necessarily reflect those endorsed by BNP Paribas Cardif.

The fact that the rmse and mse are larger than in the case of the orange area is due to the fact that less data satisfy the condition d k=1 R k j ≤ 0 and thus the estimation is less accurate.

Large number of lines of business

We consider 10 lines, with different situations. We begin with some generalization of the previous models and then propose a model for which we have a block of 5 correlated lines of business, and a block of lines which are mutually independent. 4.1. Generalizations of the previous models in dimension 10.

• Independent exponential: X 1 1 ; E(µ 1 ), for i = 2, . . . , 7, X i 1 ; E(α 1 µ 1 ) and for i = 8, . . . , 10, X i 1 ; E(α 2 µ 1 ). • Conditional exponential: X 1 1 ; E(Θ), for i = 2, . . . , 7, X i 1 ; E(α 1 Θ) and for i = 8, . . . , 10, X i 1 ; E(α 2 Θ), where Θ ; Γ(α, b).

). The last two models have the same margins. We have chosen α 1 = 5 and α 2 = 8 and u = 80. We have performed our stochastic algorithm 10 times on data sets of length 20 000 for the orange area. As expected, the optimal allocation with respect to the violet area leads to penalize all the branches correlated to the more risky one; the independent branches are less penalized.

Time horizon p > 1, d lines of business

We finish our simulation series with several multi-periodic simulations. We shall consider p = 2, p = 3, p = 4 and d = 3 lines of business. We have performed simulation only for models with independence in time. As already mentioned, for multi-periodic cases, the premium c has to be taken into account. We fix it for each branch as 5% of the expectation of the branch. We have performed 10 simulations of length 15 000, for our three models, with α = 5, u = 30, a = 3, b = 60. The simulations are done for the orange area (indicator I 2 ) then for the orange area with the stopping time (indicator I 1 ) and finally for the violet area (indicator J). It is expected that the behavior of the stopped and the non stopped orange area are quite similar because in both cases, the indicator does not take into the event which lead to global insolvency. These multi-periodic simulations show the same kind of behavior than for p = 1. Nevertheless, we remark that for the IGPD model, with p = 4, the more risky branch becomes less allocated. This is a quite surprising result which may be explained by the fact that u = 30 in this context is quite small (recall that the premiums are added at each time n = 1, . . . , p). For larger values of u, this phenomenon does not appear anymore as it can be seen below (u = 80). We have also performed simulations for the mixed model (the first 3 lines are conditionally exponential and the last two lines are independent GPD with the parameters are the same as above).

Conclusion

In this article, we have presented some results for the allocation of capital with respect to the minimization of some risk indicators. In dimension d = 2 for one period (p = 1), explicit formulas and some asymptotics are given for specific models, a simulation study is done, it allows to benchmark the parameters of the optimization algorithm. Then, simulations are done in higher dimension and for several periods (p > 1). It emphasizes the fact that minimization driven allocations require to chose the right risk indicator. As an example, we would recommend to use the orange area for free capital allocation, while the violet area is more consistent for economic capital allocation.