Dominique Guegan
email: dguegan@univ-paris1.fr

Fatima Jouad
email: fatima.jouad@gmail.com

Dominique Guégan

Documents de

The advent of the Internal Model Approval Process within Solvency II and the desirability of many insurance companies to gain approval has increased the importance of some topics such as risk aggregation in determining overall economic capital level. The most widely used approach for aggregating risks is the variance-covariance matrix approach. Although being a relatively well-known concept that is computationally convenient, linear correlations fail to model every particularity of the dependence pattern between risks. In this paper we apply different pair-copula models for aggregating market risks that represent usually an important part of an insurer risk profile. We then calculate the economic capital needed to withstand unexpected future losses and the associated diversification benefits. The economic capital will be determined by computing both 99.5 th VaR and 99.5 th ES following the requirements of Solvency II and SST.

Introduction

Insurance companies face a multitude of risks that could cause financial losses. Economic capital is the amount of capital that is needed to cover losses at a certain risk tolerance level. It captures a wide range of risks such as financial, insurance and operational risks and expresses all of this as a single number. Recently, within the regulatory framework of Solvency II, EIOPA addressed risk aggregation as one of the main challenges in risk management. Risk aggregation is about studying aggregate risks of a given portfolio. If we consider a set of risks (e.g. returns or losses in a portfolio), denoted by a random vector X = (X 1 , . . . , X d ), the aim is to build an aggregation function of those risks, denoted g(X), where g : R d → R and to measure the risk of the quantity F (g(X)), for some risk measure F over a given period. For the purpose of economic capital allocation, one can employ different approaches to aggregate risks. The most commonly used method is the variance-covariance matrix approach. While being easily implemented and convenient to work with, this approach can lead to a misestimation of the aggregated economic capital. Hence, a copula-based approach is chosen to achieve this purpose. However, one major shortcoming of copulas is their use in high dimensions. Indeed, elliptical copulas can be easily extended to higher dimension, but they are unable to model financial asymmetries (see [START_REF] Patton | Copula-based models for financial time series[END_REF], while Archimedean copulas are not satisfactory to describe multivariate dependence in dimensions higher than two [START_REF] Joe | Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters[END_REF]. In this paper, using the method of pair-copula models introduced first by Joe (1997), we analyze the impact of aggregating risks for the purpose of computing the economic capital through the calculation of two risk measures, the Value-at-Risk (99.5 th VaR) -required by Solvency II -and the Expected Shortfall (99 th ES) -required by the Swiss Solvency Test (SST). Furthermore, we compare the economic capital resulting from the different copula models with the most widely used approach in the insurance industry: the variance-covariance matrix and the simple summation approach. The two last approaches are introduced in the last section of this paper. The diversification benefits derived from the different methodologies are then discussed. Papers in the insurance area, where the application deals with risk aggregation are quite few. Most of them compare the copula methodology with linear aggregation approaches, thereby documenting how the latter may underestimate total risk (Schuermann and Rosenberg 2006, Tang and[START_REF] Tang | Economic Capital and the Aggregation of Risks using Copulas[END_REF]. Among the other approaches reviewed in the risk aggregation area, the methodology we have chosen has not been applied yet in the insurance area even if this approach is not new. Indeed, pair-copulas are used to achieve our purpose. Risk aggregation provides necessary information that enables effective group-wide risk management, as well as a wide variety of other key business decisions and business processes. However, the financial crisis that began in 2007 highlighted at least some degree of failure of risk aggregation methods. Many of the firms acknowledge that model risk in this area may be higher than previously recognized. Despite that recognition, there has been surprisingly little movement by most of these firms to reassess or revise risk aggregation practices in significant ways. Actually, in many Solvency II internal models, dependence is modeled with correlation matrix or Gaussian copula. As the dimension is large in general, the correlation matrix or the Gaussian copula remain the easiest tool (it is sometimes replaced by a Student copula). Nonetheless, other alternatives may be provided such as the use of pair-copula models [START_REF] Joe | Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters[END_REF]Bedford andCooke, 2001-2002; Guégan and Maugis, 2010 and references therein). This represents a radically new way of constructing complex multivariate highly dependent models, which parallels classical hierarchical modeling (Green et al., 2003and Guégan and Hassani, 2010[START_REF] Guégan | Multivariate VaRs for operational risk capital computation : a vine structure approach[END_REF].

A key contribution of our paper is to provide existing analytical techniques around the dependence modeling that are supervisory expectations in the context of Solvency II. First, we assess several types of vine and nested copula models. This is accomplished through the application of canonical maximum likelihood (CML) estimation methodology. We use empirical marginal distributions and calibrate the copula parameters by maximum pseudo-likelihood estimation (MPLE). We also conduct benchmarking of models by comparing two risk measures and the associated diversification benefits across different frameworks for aggregating risks. Indeed, we compare various copula models, as well as other methods for the construction of a joint distribution of losses, such as the variance-covariance matrix approach and the simple summation.The empirical exercise uses financial data collected from Bloomberg to proxy for various market risk types. The empirical analysis focuses on five of the largest market sub-risks available on a daily frequency (Equity, Interest rates, Spread, Foreign Exchange and Implied Volatility). The rationale for concentrating on these risks is motivated in large part by the fact that they usually represent more than 75% of the market risk portfolio and about 60% of the total risk portfolio. While the focus herein is upon the insurance sector, our methodology could just as easily be extended to any other kind of financial institution such as a bank. This study is part of a literature that performs a comprehensive analysis around how to combine a set of risk factors influencing the total risk of large insurance companies. The five main market 2 Documents de Travail du Centre d'Economie de la Sorbonne -2012.31 sub-risks are considered from January 1, 1999 to March 31, 2009.. Each data set is filtered by applying AR-likeGARCH process using skewed elliptic and GED innovations. We then determine the dependence structure between these five proxies applying the pair-copula methodology. Two kinds of dependence prevail: the Student-t copula and the Student C-vine. We eventually provide the aggregated economic capital through 99.5 th VaR and 99 th ES and determine the corresponding amount of diversification benefits for each approach. The rest of the paper is structured as follows. Section 2 provides the methodology used to model the dependence structure. Section 3 deals with the modeling and the filtering process of market risks needed to build and calibrate the pair-copula models produced in Section 4 where the numerical results of the estimation are presented. Section 5 provides a comparison of the economic capital and the diversification benefits resulting from the assuming different copula models to the variance-covariance matrix and simple summation approach. We conclude this paper in Section 6 with few remarks on the findings and an outlook for future research on the subject of aggregating risks for the purpose of setting capital requirements.

Methodology used to Model Dependency

A common approach for deriving economic capital is to first assess the individual risk components and then consider possible techniques to aggregate these components to derive an overall capital number. Insurance companies differ in their approaches to economic capital risk aggregation, some techniques being more sophisticated than others. The most commonly used methodologies are simple summation, variance-covariance matrix and copula-based approaches. The method we have chosen to calculate the aggregate economic capital is the pair-copula approach [START_REF] Joe | Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters[END_REF], which is more flexible than the use of the variance-covariance matrix. While the latter only requires measurement of the sub-risks' economic capital estimates and the correlation between the sub-risks, copula methods of aggregation depend on the whole distribution of the sub-risks. However, the full marginal risk distribution is needed for each risk, rather than the relevant capital number in the case of the variance-covariance matrix approach. In the simple case of two risks X and Y , a copula C : [0, 1]2 → [0, 1] is a joint cumulative distribution function (cdf) of two random variables which are uniformly distributed on [0,1]. A two-dimensional cdf H with continuous margins F 1 and F 2 can be written as H(x 1 , x 2 ) = C(F 1 (x 1 ), F 2 (x 2 )). In other words, a copula is the description of a multivariate dependence structure. It is the link between the one-dimensional marginal distribution and the n-dimensional joint distribution in the general case. This approach can be extended to aggregate more than two marginal distributions into a joint multivariate distribution. In the sequel, a short reminder about the pair-copula decomposition is given. For more details about the models and its estimation inference, we refer to [START_REF] Joe | Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters[END_REF], Bedford andCooke (2001, 2002), Kurowicka and Cooke (2006) and [START_REF] Guégan | New prospects on Vines[END_REF].The decomposition of a multivariate distribution in a cascade of pair-copulas was originally proposed by [START_REF] Joe | Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters[END_REF], and later discussed in details by the cited authors. The method of extending two-dimension copulas to higher dimensions is not very flexible and additional assumptions are often needed. In this section we use the results of [START_REF] Joe | Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters[END_REF] and follow his exposition very closely. The basic idea behind the pair-copula construction is to decompose an arbitrary distribution function into simple bivariate building blocks that are two-dimensional copulas and stitch them together appropriately. Besides, this method is recursive by nature. Consider the random vector X = (X 1 , . . . , X n ) whose joint density f (x 1 , . . . , x n ) can be factorized as hereafter:

f (x 1 , . . . , x n ) = f (x n ) • f (x n-1 |x n ) • f (x n-2 |x n-1 , x n ) • • • • • f (x 1 |x 2 , . . . , x n ) (1)
From Sklar's theorem we know that there exists a copula C such that:

F (x 1 , . . . , x n ) = C(F 1 (x 1 ), . . . , F n (x n )) (2)
Differentiating this formula by applying the chain rule and working with densities, one gets:

f (x 1 , . . . , x n ) = c 12...n (F 1 (x 1 ), . . . , F n (x n )) • f 1 (x 1 ) • • • • • f n (x n ) (3)
for some (uniquely identified) n-variate copula density c 1...n (.). For the base case in two dimensions we can easily see that the density function f (x 1 , x 2 ) is given by

f (x 1 , x 2 ) = c 12 (F 1 (x 1 ), F 2 (x 2 )) • f 1 (x 1 ) • f 2 (x 2 ) (4)
This follows immediately by taking partial derivatives with respect to both arguments in

F (x 1 , x 2 ) = C(F 1 (x 1 ), F 2 (x 2 )
), where C is the copula associated with F via Sklar's Theorem. Besides, notice that from the above formula the conditional density of X 2 given X 1 can be determined:

f 2|1 (x 2 |x 1 ) = f (x 1 , x 2 ) f 1 (x 1 ) = c 12 (F 1 (x 1 ), F 2 (x 2 )) • f 2 (x 2 ) (5)
In Section 4, we use the previous formula to represent this conditional distribution function when x and ν are uniform, i.e.f (x

) = f (ν) = 1, F (x) = x and F (ν) = ν. That is, F (x|ν, θ) = δC xν (x, ν, θ) δν (6)
where the second parameter always corresponds to the conditioning variable and θ denotes the set of parameters for the copula of the joint distribution function of x and ν. For sampling from pair-copula distributions, one needs the inverse of F (x|ν) with respect to the first variable x. In Table 1 we give some examples of that function that have been used in our case study.

Copulas We present now other ways to compute practically a high dimensional distribution: Fully nested Archimedean Copulas (FNAC) and vines (more particularly, C-vines).

F (u 1 |u 2 , θ) function Bivariate Gaussian F (u 1 |u 2 , ρ 12 ) = φ( φ -1 (u1)-ρ12φ -1 (u2) √ 1-ρ 2 12 ) Bivariate Student F (u 1 |u 2 , ρ 12 , ν 12 ) = t ν12+1 ( t -1 ν 12 (u1)-ρ12•t -1 ν 12 (u2)
(nu 12 +(t -1 ν 12 2 ))(1-ρ 2 12 ) ν 12 +1 Bivariate Clayton F (u 1 |u 2 , ρ 12 ) = u -δ12-1 2 (u -δ12 1 + u -δ12 2 -1) -1-1 δ 12 Bivariate Gumbel F (u 1 |u 2 , ρ 12 ) = C 12 (u 1 , u 2 ) • 1 u2 • (-log(u 2 )) δ12-1 • ((-log(u 1 )) δ12 + (-log(u 2 )) δ12 ) 1 δ 12 -1

Nested Copulas

The idea behind FNAC models is simply to add a dimension step by step. For the seek of illustration, take for instance the nodes u 1 and u 2 coupled through copula C 1 , node u 3 coupled with C 1 (u 1 , u 2 ) through copula C 2 , node u 4 coupled with C 2 (u 3 , C 1 (u 1 , u 2 )) through copula C 3 , and finally node u 5 coupled with C 3 (u 4 , C 2 (u 3 , C 1 (u 1 , u 2 ))) through copula C 4 . Hence, the copula for the 5-dimensional case requires four bivariate copulas C 1 , C 2 , C 3 and C 4 , with corresponding generators φ 1 , φ 2 , φ 3 and φ 4 :

C(u 1 , . . . , u 5 ) = C 4 (u 5 , C 3 (u 4 , C 2 (u 3 , C 1 (u 2 , u 1 )))) = φ -1 4 (φ 4 (φ -1 3 (φ 3 (φ -1 2 (φ 2 (φ -1 1 (φ 1 (u 1 ) + φ 1 (u 2 )) + φ 2 (u 3 ))) + φ 3 (u 4 ))) + φ 4 (u 5 )))
In this structure, all margins are themselves Archimedean copulas. It allows for the specification of d-1 copulas and corresponding distributional parameters, while the remaining (d-1)(d-2)/2 copulas and parameters are implicitly given through the construction. More specifically, the two pairs (u 1 , u 3 ) and (u 2 , u 3 ) both have copula C 2 with dependence parameter θ 2 . Moreover, the three pairs (u 1 , u 4 ), (u 2 , u 4 ) and (u 3 , u 4 ) all have copula C 31 with dependence parameter θ 3 . Hence, when adding variable x to the structure, we specify the relationships between x pairs of variables.

The FNAC is a construction where some technical conditions need to be satisfied for being a proper d-dimensional copula. The consequence of these conditions for the FNAC, if all the generators are of the same type, is that the degree of dependence, as expressed by the copula parameter, must decrease with the level of nesting in order for the resulting d-dimensional distribution to be a proper copula. If the generators belong to different families, the parameter restrictions are even stronger. There are only few generators that can be combined and this method often fails in practice. For nested copula models, the parameters may be estimated by maximum likelihood estimation.

Vines

There exists many ways of building a pair-copula model. This leads to a large number of possible pair-copulas constructions, for instance in the case of five variables there are 240 different constructions. In order to organize them one can use a graphical model denoted as the regular vine introduced by Bedford andCooke (2001-2002). In the sequel, we are going to focus on one special case of these regular vines: C-vines (Kurowicka and Cooke, 2004). Fitting a canonical vine to our data set is relevant in our case since a particular variable is known to be a key variable that governs interactions in the data set. In such a situation we decide to locate this variable at the root of the canonical vine. In a case of a d-dimensional C-vine model, the pairs at tree 1 are 1, i, for i = 2, . . . , d, and for tree T (2 ≤ T ≤ d), the conditional pairs are T, i|1, . . . , T -1 for i = T + 1, . . . , d. That is, for the C-vine, conditional copulas are specified for variables T and i = 1, . . . , T -1.

In other words, each tree T j has a unique node that is connected to n -j edges. We have then n possibilities to choose the key variable in the first tree, n -1 possibilities in the second tree,... and 3 possibilities in tree n -2. Thus, there are n! 2 different C-vines. The general expression for the five-dimensional C-vine model is given hereunder:

f (x 1 , x 2 , x 3 , x 4 , x 5 ) = f (x 1 ) • f (x 2 ) • f (x 3 ) • f (x 4 ) • f (x 5 )• c 12 (F (x 1 ), F (x 2 )) • c 13 (F (x 1 ), F (x 3 )) • c 14 (F (x 1 ), F (x 4 )) • c 15 (F (x 1 ), F (x 5 ))• c 23|1 (F (x 2 |x 1 ), F (x 3 |x 1 )) • c 24|1 (F (x 2 |x 1 ) • F (x 4 |x 1 )) • c 25|1 (F (x 2 |x 1 ), F (x 5 |x 1 ))• c 34|12 (F (x 3 |x 1 , x 2 ), F (x 4 |x 1 , x 2 )) • c 35|12 (F (x 3 |x 1 , x 2 ), F (x 5 |x 1 , x 2 ))• c 45|123 (F (x 4 |x 1 , x 2 x 3 ), F (x 5 |x 1 , x 2 , x 3 ))

Estimation Methodology

For evaluating the log-likelihoods in pair-copula models, two assumptions are needed. The observations of each variable are independent over time and the variables are uniform in [0, 1]. Then, the following steps should be taken when fitting a pair-copula decomposition.

• A specific factorization is defined. We determine which factors are most dependent by using the Kendall's tau or Spearman's rho coefficients. These ones should be coupled. In case of Student-t copulas these can also be the factors where the bivariate Student-t copula estimation reveals the smallest number of degrees of freedom. In case of a C-vine this can be performed in each tree, or one uses the order of covariates in the first tree to determine the order of nodes (which will be used here).

• Pair-copula types are chosen. The so far simplest and most common way to find out which copula type to use is the graphical analysis.

• The copula parameters are estimated via maximum likelihood estimation, see [START_REF] Guégan | New prospects on Vines[END_REF] for more details about the inference procedure. To do so, we first estimate the parameters of the copulas in tree 1 from the original data. Then, we compute the observations (i.e. conditional distribution functions) for tree 2 using the copula parameters from tree 1 and the marginal conditional distributions. Besides, we iterate steps 1 and 2 over all trees. In each tree, we determine which copula types to use.

• We maximize the full log-likelihood function using the parameters obtained from the stepwise procedure as starting values.

Dynamics of Market risks

In this section, we summarize the filtering process of the five market risks needed to perform the pair-copula modeling.

Let consider the daily Euro Stoxx index for Equity risk, denoted (X 1 ) t , the daily Euro 10 year-IR swap rates for interest rates risk, denoted (X 2 ) t , the daily Euro spread figures for spread risk, denoted (X 3 ) t , the daily Euro Implied Volatility rates for Implied volatility risk, denoted (X 4 ) t and the daily USD/EUR Foreign Exchange rates for FX risk, denoted (X 5 ) t . The data set starts from January 1, 1999 to March 31, 2009. These samples cover from 10 years and encompass several world crises such as the 9/11 terrorist attack of 2001, the mortgage crisis of 2007 and the huge credit and market turmoil of 2008. This data set have been collected from Bloomberg.

To achieve stationarity, the log returns of the time series, denoted (R i ) t for i = 1, . . . , 5 are computed: R t = ln((X t )/(X t-1 )). The Augmented Dickey Fuller (ADF) test performed for all data sets confirm the stationarity of the series. Table 2 summarizes the descriptive statistics of the five market risk returns (R 1 ) t , (R 2 ) t , (R 3 ) t , (R 4 ) t and (R 5 ) t , the ADF test and the normality tests. From these findings, all data sets appear to have heavy tails and the normality assumption is inconsistent with all time series. When performing a model selection in the scope of the generalized hyperbolic distribution class based on the Akaike information criterion, we obtain as the best model the symmetric variance gamma distribution (defined as the normal variance-mean mixture where the mixing density is the gamma distribution) for Equity, IR, Implied Volatility and FX risk data while the best model obtained for Spread risk is a symmetric generalized hyperbolic distribution function. The estimated parameters and the Kolmogorov-Smirnov statistics, given in the appendix, confirm theses findings. Hence, as normality assumption is clearly not consistent with all time series, we use Quasi Maximum Likelihood (QML) for estimating the parameters of the AR-likeGARCH models in the following section.

Descriptive Statistics

Statistics (R 1 ) t (R 2 ) t (R 3 ) t (R 4 ) t (R 5 ) t Mean -1.

Modeling and Filtering Process

When plotting the trajectory for the series (R 1 ) t ,(R 2 ) t ,(R 3 ) t ,(R 4 ) t and (R 5 ) t on the full period from 1999 to 2009 (sample size n=2673), the plots show that the log-returns were more volatile toward the end of the study period. These results might be triggered by the instability in the markets due to the market turmoil starting from 2007. One can expect that the conditional variance of all these series vary over time and hence use GARCH models. Some disadvantage of GARCH processes in modeling financial returns is the symmetry of conditional variance of t with respect to positive and negative values of t-1 , t-2 ,. . . . In practice, one observes a leverage effect, meaning asymmetric consequences of positive and negative innovations (the conditional variance tends to decrease if noise is positive implying bigger returns). An Exponential GARCH process (Nelson, 1991) does not have this disadvantage. This is also the case for Asymmetric

Power ARCH process introduced by Ding, Granger and Engle (1999). Therefore, in order to filter our market risk data, we apply then AR(k)-likeGARCH(p,q) process:

R t = φ 1 R t-1 + • • • + φ k R t-k + t , t = 1, . . . , T. (7)
and

t = σ t ν t (8)
where (ν t ) is a strict white noise, and (σ t ) satisfies the following recurrence equations:

• GARCH(p,q) (Bollerslev, 1986)

σ 2 t = ω + p j=1 α j 2 t-j + q j=1 β j σ 2 t-j (9)
where ω denotes the variance intercept, α j the ARCH(p) parameters, β i the GARCH(q) parameters and σ 2 the conditional variance,

• EGARCH(p,q) (Nelson, 1991)

log e (σ 2 t ) = ω + p j=1 (α j z t-j + γ j (|z t-j | -E|z t-j |)) + q j=1 β j log e (σ 2 t-j ) (10)
where the coefficients α j captures the sign effect, γ j the size effect and z t are standardized innovations.

• APARCH(p,q) (Ding, Granger and Engle, 1993)

σ 2 t = ω + p j=1 α j (| t-j | -γ j t-j ) δ + q j=1 β j σ δ t-j (11)
where the power parameter δ ∈ R + , being a Box-Cox transformation of σ t , and γ j the leverage(p) parameters. Various models arise from this model such as GARCH (Bollerslev, 1986), AVGARCH (Taylor, 1986 and Schwert, 1990), GJRGARCH (Glosten et al., 1993), TGARCH (Zakoian, 1994), Non linear ARCH (Higgins and Bera, 1992) and the Log ARCH (Geweke, 1986 and Pantula, 1986).

Thus, V ar( t | t-1 , t-2 , . . . ) = σ 2 t , E( t ) = 0, Cov( t , s ) = 0, t = s.
The distribution of the residuals ν t can be Gaussian, Student, GED or the generalized hyperbolic (GH) distributions (Barndorff-Nielsen,1977). Different AR models are fitted to our time series by applying QML estimation and AIC is used for selecting the number of AR parameters. it comes out that an AR(3) process, an AR(4) process and an AR(6) process can describe respectively the observed returns (R 2 ) t , (R 3 ) t and (R 4 ) t , and (R 1 ) t while (R 5 ) t can be modeled by a simple white noise. In Table 3 we provide the estimated parameters of the AR part with the corresponding standard error. For further analysis 4 GH distributions can model asymmetry and heavytailedness ant its density is written as:

gh(x|λ, α, β, δ, µ) = a(λ, α, β, δ, µ)(δ 2 + (x -µ) 2 ) (λ-1/2)/2 K λ-1/2 (α δ 2 + (x -µ) 2 )exp((β(x -µ)) with a(λ, α, β, δ, µ) = (α 2 -β 2 ) λ/2 /( √ 2πα λ-1/2 δ λ K λ (δ α 2 -β 2 )
) where K λ is the modified Bessel function of the third kind, δ ≥ 0 and 0 ≤ |β| < α. It is symmetric if β = 0. The normal inverse Gaussian is the special case λ = -1/2. The hyperbolic distribution is the special case λ = 1 and many distributions are limiting cases: normal, student-t and variance-gamma -see Prause (1999) or Eberlein and Von Hammerstein (2003). of residuals ˆ t behavior, volatility models listed previously with skewed Normal, skewed Student and skewed GED have been fitted using QML method. The estimated parameters, standard errors and information criteria (AIC, BIC and HQIC) are given in the sequel (Tables 4 to 8). Note that if more than one model have comparable values of these criteria, then the model with lower number of parameters or lower p-values of significance parameters' tests is chosen. As a result, the best fitting models for each market risk are given hereunder.

• The daily EuroStoxx returns will be modeled using an AR(6) -EGARCH(2, 1) process with a skewed student-t distribution.

• The daily 10y-IR returns will be modeled applying an AR(3) -EGARCH(1, 1) process with a skewed GED distribution.

• The daily Euro Spread returns will be modeled employing an AR(4) -EGARCH(1, 1) process with a skewed GED distribution.

• The daily Euro Implied Volatility returns will be modeled using an AR(4)-EGARCH(1, 2) process with a skewed student-t distribution.

• The daily USD/EUR FX returns can be modeled applying an EGARCH(2, 1) process with a skewed GED distribution.

These modelings provide the residuals (ν t ) t for each market risk that we use in the sequel in order to characterize the dependence structure between these five market risks through the differnet copula approaches. In order to account for possible time-dependence and volatility on one hand and to fulfill the assumptions of maximum likelihood used for selecting and calibrating copula models on the other hand, the financial time series for Equity, Interest Rate, Spread, Implied volatility and FX risks have been filtered by applying AR-likeGARCH models. Given the resulting residuals, respectively denoted ((ν 1 ) t , . . . , (ν 5 ) t ), the joint density of several copula models are estimated by using a semi-parametric method called Canonical Maximum Likelihood estimation (CML). Indeed, the residuals (ν t ) t are transformed via the empirical cumulative distribution function into uniform random variables U t = F N ((ν t )) on the unit hypercube. Note that the filtering process do not alter the dependence structure when using copulas. The resulting "pseudo observations" are used to estimate the unknown parameters of the copula models by using MPLE. Goodness-of-fit measures are the log-likelihood and BIC criterion (to account for the number of parameters). Given their properties, the parametric copulas used in our study are the following:

• Copulas with no tail dependence such as the normal copula,

• Copulas with upper tail dependence only such as the Gumbel copula,

• Copulas with lower tail dependence only such as the Clayton copula,

• Copulas with symmetric upper and lower tail dependence such as the t copula.

In the sequel, we are going to treat two ways of modeling the dependence structure: the variance-covariance matrix approach which is the most currently used approach in the industry and copulas. Among the copula methodology, we consider 5-variates copulas, FNAC and C-vines. The maximized log-likelihood function and BIC criterion for each copula model are given in Tables 12, 13 and 16. Note that in order to have a positive dependence overall and an appropriate copula fit, the negative residuals for Spread and Implied Volatility data are taken. Figures 9 and 10 display the plot of the residuals (ν i ) t for i = 1, . . . , 5 and the plot of the normalized ranks in the unit hypercube.

Variance-Covariance Matrix Approach

The variance-covariance matrix approach is a convenient and commonly used analytical technique that allows managers to combine marginal distributions of losses or distinct tail losses into a single aggregate loss distribution. The sole requirement is to characterize the level of interdependence of standalone losses, which is typically accomplished with a matrix of linear correlations. The lower the correlations on the non-diagonal elements of the matrix, the greater the level of diversification that can be realized with incremental (long) exposure to a risk component. The variance-covariance matrix approach is well known from multivariate normal distributions where the value-at-risk (VaR) of a risk category X := X 1 + • • • + X n is calculated based on the VaRs of the single risks and their correlation matrix Σ. In order to introduce this approach, we start with the formula for the standard deviation. This formula holds true for all random variables with finite variances. The standard deviation σ of X is calculated as

σ 2 = (σ 1 , . . . , σ n )Σ(σ 1 , . . . , σ n ) T ,
where σ i denotes the standard deviation of the single risk X i . When (X 1 , . . . , X n ) are multivariate normally distributed and centered, a similar formula can be used for calculating the value-at-risk of X:

V aR(α) = n i=1 V aR 2 i (α) + 2 • n i=1 n j=1 ρ i,j V aR i (α)V aR j (α) (12)
where ρ i,j is the correlation between risks (here the five assets) i and j and α is the confidence level of the aggregated VaR. More generally, this formula holds true for all centered elliptical distributions. Typical examples of elliptical distributions are the Normal and the Student-t distributions. The latter are quite flexible in what concerns the modeling of heavy tails and the corresponding extreme events. This aggregation approach is based on the assumption that risks are measured as the difference between value-at-risk and the mean.

The main advantages of this approach are that it is simpler relative to other methods, can be evaluated through a simple formula and does not require fundamental information about the lower-level risks. However, while the variance-covariance matrix approach is simple and highly tractable approach to risk aggregation, the cost to the unwary user is that it effectively fills in unspecified details about the nature of the loss distributions, which may or may not be accurate or intended (e.g. assumption of elliptical distributions). The variance-covariance matrix method imposes thus a simple dependence structure on what is believed to be a more complex web of dependences.

In order to derive the economical capital through VaR and ES (the risk measures we consider in this paper) from the variance-covariance matrix approach, we used an empirical correlation matrix inferred from the historical data.

Economic capital and Diversification benefits

In this section, the total market economic capital is given through the calculation of the VaR and the ES from the previous copulas modelings. The risk measures are computed based on a specified level of probability. For our purposes, we have chosen q = 99.5% for the VaR and q = 99% for the ES. These confidence levels are chosen to be respectively consistent with Solvency II and SST requirements.

In order to achieve our purpose, we have created a virtual portfolio consisting of 61 underlying investments representing the five risks with respective weights ω 1 , . . . , ω 61 so that the change in value of the portfolio over the given holding period (the so-called P&L or profits and losses) can be written as X = Σ i ω i X i , where X i denotes the change in value of the i th investment.

Measuring the risk of this portfolio essentially consists in determining its distribution function in order to derive then the functionals describing this distribution function such as its mean, the 99 th quantile and the 99.5 th quantile. Knowing that, let us consider this virtual portfolio and a fixed time horizon, and denote by F L (l) = P (L ≤ l) the distribution of the corresponding loss distribution. Given some confidence level α ∈ (0, 1), the VaR of the portfolio at the confidence level α is given by the smallest number l such that the probability that the loss L exceeds l is no larger than (1 -α). Formally,

V aR α = inf (l ∈ R : P (L > l) ≤ 1 -α) = inf (l ∈ R : F L (l) ≥ α). (13)
For an integrable loss L with continuous distribution function F L and any α ∈ (0, 1) we have ES α = E(L|L ≥ V aR α ). Note that this risk measure is coherent and reflects the severity of the losses which is not the case for the VaR measure (Artzner et al, 1997). Now, let denote by µ the mean of the loss distribution. The statistics V aR α -µ and ES α -µ are used for economic capital purposes instead of ordinary VaR and ES measures.

In our study, the VaR measure is calculated as the point on the ranked P&L distribution that corresponds to the particular level of q. That is, for each output distribution, of the 110000 simulated P&L, the 99.5 th VaR corresponds to the 550 th worst value if the distribution is ranked in decreasing magnitude. To calculate the ES, we simply take the arithmetic average, or expected value, of the values subsequent to the value corresponding to the VaR measure. Therefore, 99 th ES is calculated as the average of the 1 st to 1100 th worst value of the ranked P&L distribution. Continuing with the procedure outlined earlier, the simulated P&L for each market sub-risk is then aggregated to produce a distribution of the aggregate P&L under each copula model. The single 99.5 th VaR and 99 t hES for the five market sub-risks are displayed in Table 17.

We then conduct benchmarking of these models by comparing the two risk measures and the associated diversification benefits across the different risk aggregation approaches. Actually, we compare the estimated copula models with the variance-covariance matrix approach and the simple summation. We also examine the resulting impact in terms of economic capitals when using these different approaches models. We then present the resulting economic capitals hereafter in monetary units. To make the assessment, Table 18 exhibits the economic capitals for each linear approach and for each pair-copula model and their corresponding diversification benefits in percentage.

Concluding Remarks

From Table 18, we observe that the simple summation approach is the most conservative method since it does not allow for diversification, which is equivalent to say that the risks are perfectly positively correlated. Then comes the C-vine copula approach followed then by the variancecovariance matrix method. We also note that the use of 5-variates copulas leads to lower aggregated economic capital than when using the variance-covariance matrix approach. When focusing on the copula methodology, we see that there is a consistent effect due to the choice of the structure on the total economic capital regardless of the risk measure used. The most conservative copula approach is the C-vine approach. This can be explained by the fact that we have a more appropriate structure for gathering the variables and more parameters to describe the dependence pattern than when applying the other copula approach.

The choice of copula has also a profound effect on the resulting economic capital and the diversification benefits. As we expected, due to its lack of account of the lower tail dependence of losses, the Gumbel copula results in the lowest capital requirement while the Clayton and the Student copulas result in the highest capital requirement. From these findings, we can conclude that among the copula approaches, the Clayton C-vine copula results in the highest capital requirement (which makes sens because it allows for lower tail dependence), followed by the Student C-vine. This difference is significant for both risk measures with the extreme being the 20 Documents de Travail du Centre d'Economie de la Sorbonne -2012.31 ES (99%) case. We can also underline that the observed discrepancy of the economic capital between copulas highlights the importance of correctly modeling the dependence structure, in particular the structure and tail dependence between losses across risks. We make the following key observations: In Table 18, the value of the diversification benefits according under each copula is calculated and this provides an initial crude measure of the magnitude of the diversification benefits resulting from each copula. In decreasing order, the ranking of the copulas is 5-variates Gumbel (41.7%), 5-variates Gaussian (41.2%) and 5-variates Student (39.6%). Therefore, modeling dependencies between risks from different risks with 5-variates copulas allowing for higher tail dependence will give by far the highest level of diversification benefits while using a C-vine that allows for lower tail dependence will result in the lowest level of diversification benefits.

Then, the choice of the structure and the choice of copula have a paramount effect on the economic capital as well as on the diversification benefits for an insurer. This effect is driven mainly by the amount of tail dependence that the copula allows for losses between risks. The more lower tail dependence allowed by a copula, the higher the economic capital is if losses are aggregated under that copula. Hence, even if the best fitting copula models are the Student-t C-vine and the 5-variate Student-t copula, the copulas that exhibit lower tail dependence will yield the highest EC figures. This result makes sense due to the nature of these copulas. For all copula assumptions except the C-vine model, we see that they underestimate the economic capital for the virtual portfolio. This means that if insurers are willing to use a 5-variates copula based model rather than a vine method to aggregate their risks for economic capital purposes, they may not cover their risks to withstand unexpected losses. Therefore, in making that assertion, we must be aware that this result is sensitive to the choice of the copula. Consequently this leads to the need for more accurate and flexible internal models to be developed for capital determination purposes where copulas are incorporated for the aggregation of losses from different risks. However, as the actual economic capital are found to be extremely sensitive to the choice of copulas, and further this sensitivity varies depending on the portfolio composition of the particular insurer, it is imperative that in constructing the internal models, an appropriate copula assumption and an appropriate structure are made.

Conclusion

To ensure solvency, insurers are required both for regulatory purposes and as a going business concern to hold capital to back their insurance liabilities. In aggregating losses from different risks for the purpose of capital determination, insurers have traditionally either ignored the dependence structure between risks or used simple linear correlations to model such dependence.

In this paper, we aggregate each risk losses using the variance-covariance matrix approach, simple summation, 5-variates copulas and two variants of pair-copula models. We then assess the capital requirements in each case using the value-at-risk and tail conditional expectation risk measures at 99.5% and at 99%, respectively. Further, we calculate the diversification benefits from holding capital in aggregate risks. We analyze this diversification benefit for the different copula assumptions. The following are the key findings of this paper. First, the choice of the structure for gathering variables has a huge effect on both the economic capital and diversification benefit for an insurer. In our case study, the more sophisticated or appropriate structure a model allows (vines), the higher the required capital is. The opposite relationship between the choice of the structure and the diversification benefit exists. Therefore, because of the potential for massive modeling errors, it is imperative for insurers to select a structure that is most reflective of their own unique situation to avoid the risk of mis-calculating their capital requirement.

Second, the choice of the copula model impacts also both the economic capital and diversification benefits for an insurer. Indeed, the more tail dependence a copula allows (here Clayton and Student), the higher is the required capital. Lastly, the adequacy of the variance-covariance matrix approach is assessed against the capital requirements implied by the copula models.The C-vine Student and the C-vine Clayton appear to yield more conservative economic capital than when using the variance-covariance matrix approach. Again, the higher the copula's allowance for tail dependence and the more appropriate structure a copula model allows, the higher the resulting capital requirement will be under that model.

The reader of this paper should be aware that there are simplifying assumptions inherent in our analysis. In calculating the capital requirements and diversification benefits, we have limited our focus on the five main market sub-risks of an insurer while ignoring other sources of market risk such as for instance hedge funds, private equity and real estate risks. However, mainly due to data limitations, amongst other reasons, we have excluded the former risks from our scope. Due to these assumptions, the single most significant limitation to the results of this paper is that they merely serve to compare the effect of the dependence structure on the capital requirements and does not quantify the required capital that a particular insurer should hold for the total market risk. The numerical results in monetary terms presented in this paper must also be taken with caution as this is only representative of a virtual portfolio held by an insurer and again, we emphasize that it only accounts for the market risk component of the total capital requirements. However, the modeling procedure of aggregating risks using copulas as demonstrated in this paper can be readily adapted. Given that our findings in this paper indicate a pronounced effect of the copula structure on the capital requirements for insurers, we suggest further research into using other copulas to model the dependence structure of insurers portfolios. In this paper, we have explored using pair-copulas to model the dependence structure but these are only a few of the vast pool of copula structures that one can draw from. For example, rather than modeling in a static way the dependence between risk sources, one may investigate dynamic dependencies between these risks using dynamic pair-copula models. It would be also interesting to conduce the same analysis considering parametric margins instead of empirical ones. All these topics will be discussed in future papers.

8 |Figure 1 :

 81 Figure 1: Plots of the filtered data (left) and normalized ranks (right).

Figure 2 :

 2 Figure 2: Diagnostic plots of standardized residuals from the fitted AR6-EGARCH(2,1) model with skewed student distribution function for Equity risk data

Figure 3 :

 3 Figure 3: Diagnostic plots of standardized residuals from the fitted AR3-EGARCH(1,1) model with skewed GED distribution function for IR risk data

Figure 4 :

 4 Figure 4: Diagnostic plots of standardized residuals from the fitted AR4-GARCH(1,2) model with skewed Student distribution for Implied Volatility risk data

Figure 5 :

 5 Figure 5: Diagnostic plots of standardized residuals from the fitted GARCH(2,1) model with skewed GED distribution for FX risk data

Table 1 :

 1 Examples of the F (u 1 |u 2 , θ) function for different bivariate copulas.

Table 3 :

 3 5 AIC = -2ln(L) + 2k where k is the number of unknown parameters, n is the sample size and L denotes the maximum likelihood function 8 Documents de Travail du Centre d'Economie de la Sorbonne -2012.31 Parameter estimates of AR models for (R 1 ) t , (R 2 ) t , (R 3 ) t and (R 4 ) t with their corresponding standard deviations.

	AR Par.	Φ 1	Φ 2	Φ 3	Φ 4	Φ 5	Φ 6	AIC
	Estimates -0.038	-0.037	-0.068	0.052	-0.077	-0.049	-14677
	(R 1 ) t	(0.019)	(0.019)	(0.019)	(0.019)	(0.019)	(0.019)	-
	Estimates -0.0170	0.0461	-0.0588	-	-	-	-17186
	(R 2 ) t	(0.019)	(0.019)	(0.019)	-	-	-	-
	Estimates -0.2435	-0.0372	0.0460	0.0535	-	-	-12784
	(R 3 ) t	(0.019)	(0.020)	(0.020)	(0.019)	-	-	-
	Estimates 0.1785	0.0493	0.0634	-0.0796	-	-	-14825
	(R 4 ) t	(0.019)	(0.020)	(0.020)	(0.019)	-	-	-

Table 4

 4

	: Estimates, standard errors and information criteria for various GARCH(2,1), EGARCH(2,1) and APARCH(2,1) models fitted	to AR(6) residuals for Euro Stoxx data over the study period.

10 Documents de Travail du Centre d'Economie de la Sorbonne -2012.31 11 Documents de Travail du Centre d'Economie de la Sorbonne -2012.31

Table 6 :

 6 Estimates, standard errors and information criterion for various GARCH(1,1), EGARCH(1,1) and APARCH(1,1) models fitted to AR(4) residuals for Euro Spread data over the study period.

	APARCH	Skew GED	5.874e-6	(1.0e-6)	3.358e-1	(0.032)	3.980e-1	(0.078)	3.107e-1	(0.070)	1.663e-7	(0.034)	1.624	(0.048)	1.051	(0.024)	1.426	(0.051)	-7.3569	-7.3392	-7.3505
	APARCH	Skew t	6.297e-6	(1.0e-6)	3.819e-1	(0.036)	3.484e-1	(0.075)	3.192e-1	(0.067)	3.596e-9	(0.032)	1.653	(0.045)	1.048	(0.029)	6.965	(0.791)	-7.3680	-7.3540	-7.3617
	EGARCH APARCH	Skew GED Skew Normal	-0.182 3.244e-6	(0.055) (1.0e-6)	2.202e-3 2.671e-1	(0.004) (0.022)	0.679 4.482e-1	(0.068) (0.076)	0.303 3.091e-1	(0.069) (0.072)	0.558 4.937e-1	(0.043) (0.030)	-1.694	(0.084)	1.060 1.037	(0.023) (0.024)	1.401 -	(0.052) -	-7.3489 -7.3239	-7.3335 -7.3085	-7.3433 -7.3183
	EGARCH	Skew t	-0.241	(6.623e-2)	2.358e-3	(7.013e-3)	0.694	(6.417e-2)	0.282	(6.456e-2)	6.146e-1	(4.590e-2)	-		1.054	(2.747e-2)	6.427	(0.752)	-7.3589	-7.3435	-7.3533
	GARCH GARCH EGARCH	Skew t Skew GED Skew Normal	1.048e-6 8.668e-7 -0.101	(1.0e-6) (1.0e-6) (3.560e-2)	3.426e-1 3.017e-1 1.212e-3	(3.673e-2) (3.554e-2) (1.117e-3)	3.496e-1 4.001e-1 0.657	(7.983e-2) (9.084e-2) (6.369e-2)	3.068e-1 2.971e-1 0.331	(6.653e-2) (7.530e-2) (6.384e-2)	--0.481	--(3.348e-2)	---	-	1.045 1.047 1.047	(2.704e-2) (1.950e-2) (2.434e-2)	7.653 1.448 -	(9.320e-1) (5.675e-2) -	-7.3680 -7.3559 -7.3137	-7.3536 -7.3426 -7.3005	-7.3620 -7.3511 -7.3089
	GARCH	Skew Normal	6.574e-7	(2.0e-6)	;2.526e-1	(2.664e-2)	4.450e-1	(8.402e-2)	3.015e-1	(7.099e-2)	-	-	-	-	1.036	(2.346e-2)	-	-	-7.3240	-7.3129	-7.320
	Models		ω	(σ)	α 1	(σ)	β 1	(σ)	β 2	(σ)	γ 1	(σ)	δ	(σ)	Skew	(σ)	Shape	(σ)	AIC	BIC	HQIC

Table 7 :

 7 Estimates, standard errors and information criterion for various GARCH(1,2), EGARCH(1,2) and APARCH(1,2) models fitted to AR(4) residuals for Euro Implied volatility data over the study period.

	APARCH	Skew GED	2.379e-9	(1.0e-6)	1.146e-2	(0.015)	2.030e-2	(0.015)	9.496e-1	(0.014)	2.374e-2	(0.011)	2.371e-2	(0.030)	2.752	(0.085)	9.899e-1	(0.023)	1.575	(0.052)	-7.3625	-7.3426	-7.3553
	APARCH	Skew t	1.187e-8	(1.0e-6)	1.230e-2	(0.018)	2.233e-2	(0.014)	9.353e-1	(0.011)	2.383e-2	(0.070)	2.366e-2	(0.060)	2.725	(0.047)	9.984e-1	(0.024)	9.983	(1.322)	-7.3536	-7.3337	-7.3464
	EGARCH APARCH	Skew GED Skew Normal	-0.032 7.394e-10	(0.026) (1.0e-60)	0.012 1.332e-2	(0.011) (0.013)	-0.012 2.052e-2	(0.011) (0.007)	0.997 9.475e-1	(0.003) (0.015)	-0.108 2.367e-2	(0.054) (0.091)	0.177 2.369e-2	(0.054) (0.074)	-2.760	(0.077)	0.994 9.888e-1	(0.022) (0.022)	1.555 -	(0.065) -	-7.3658 -7.3246	-7.3482 -7.3091	-7.3595 -7.3190
	EGARCH	Skew t	-0.035	(0.021)	0.006	(0.006)	-0.007	(0.006)	0.997	(0.002)	-0.120	(0.046)	0.189	(0.047)	-		0.992	(0.025)	11.926	(2.698)	-7.3611	-7.3435	-7.3547
	GARCH EGARCH	Skew GED Skew Normal	7.374e-8 -0.030	(1.0e-6) (0.021)	2.949e-5 0.022	(4.2e-5) (0.018)	3.053e-2 -0.021	(1.496e-3) (0.018)	9.684e-1 0.997	(1.675e-3) (0.002)	--0.095	-(0.044)	-0.163	-(0.044)	--		9.930e-1 0.992	(2.270e-2) (0.024)	1.564 -	(6.207e-2) -	-7.3697 -7.3531	-7.3564 -7.3376	-7.3649 -7.3475
	GARCH	Skew t	7.743e-8	(1.0e-6)	8.549e-7	(0.019)	3.133e-2	(0.019)	9.677e-1	(0.002)	-	-	-	-	-		9.908e-1	(0.025)	12.621	(2.851)	-7.3648	-7.3516	-7.3601
	GARCH	Skew Normal	7.685e-8	(1.0e-6)	2.515e-4	(3.92e-4)	2.960e-2	(2.181e-3)	9.691e-1	(2.281e-3)	-	-	-	-	-	-	9.908e-1	(2.428e-2)	-	-	-7.3574	-7.3464	-7.3534
	Models		ω	(σ)	α 1	(σ)	α 2	(σ)	β 1	(σ)	1	(σ)	γ 2	(σ)	δ	(σ)	Skew	(σ)	Shape	(σ)	AIC	BIC	HQIC

Table 8 :

 8 Estimates, standard errors and information criterion for various GARCH(2,1), EGARCH(2,1) and APARCH(2,1) models fitted to FX USD/EUR returns data over the study period.

Table 11

 11

	displays the correlation matrix between the

Table 13 :

 13 Values of the maximized log-likelihood function and BIC criteria for the elliptic and Archimedean C-Vine fitted to the transformed data set.

Table 14 :

 14 Single VaR(99.5 th ) and ES(99 th ) for the five market sub-risks (in EUR m).

	Market risks	VaR(99.5 th ) ES(99 th )
	Equity		4840	4918
	Interest Rates	872	902
	Spread		2341	2453
	FX		383	393
	Implied Volatility	1897	1916
	Methodologies	(99.5 th )V aR (99 th )ES Diversification Benefits
	Variance-covariance approach	6800	6956	34.2%
	Simple summation	10334	10582	0%
	Gaussian	6063	6220	41.2%
	Student	6188	6387	40.0%
	Gumbel	6026	6173	41.7%
	Clayton	6276	6444	39.2%
	C-Vine Student	7303	7502	29.2%
	C-Vine Clayton	7467	7674	27.6%
	C-Vine Gumbel	7159	7355	30.6%

Table 15 :

 15 Single VaR(99.5 th ), ES(99 th ) (in EUR m) and diversification benefits (in %) for both linear and copula approaches.

  Estimated parameters and Kolmogorov-Smirnov statistics of VG distribution for (R 1 ) t , (R 2 ) t , (R 4 ) t and (R 5 ) t and of GHYP distribution for (R 3 ) t .

						S test (p-value)
	(R 1 ) t	4.595e-11 8.025e-1 1.554e-2 0	-	0.110 (0.197)
	(R 2 ) t	5.183e-11	0.5	9.501e-3 0	-	0.121 (0.117)
	(R 3 ) t	5.902e-13 -7.485e-2 2.969e-1 0 1.772e-12	0.114 (0.192)
	(R 4 ) t	2.026e-11	0.5	8.003e-3 0	-	0.122 (0.112)
	(R 5 ) t	-6.190e-5	1.903	6.511e-3 0	-	0.124 (0.103)
	Figure 6:

The risks are represented by the daily Euro Stoxx index for Equity risk, the daily Euro 10 year-IR swap rates for interest rates risk, the daily Euro spread figures for spread risk, the daily Euro Implied Volatility rates for Implied volatility risk and the daily USD/EUR Foreign Exchange rates for FX risk

Documents de Travail du Centre d'Economie de la Sorbonne -2012.31

Documents de Travail du Centre d'Economie de la Sorbonne -2012.31

Graphs are provided in the appendix at the end of the paper.

* Research supported by Fonds AXA pour la Recherche. 1 advisory body to the European Parliament and the Council of the European Union.

insurance policyholders, pension scheme members and beneficiaries.

Table 9: Correlation matrix for market risks.

be intuitively expected. Note that the simple summation method is a particular case of the variance-covariance matrix approach when all correlation coefficients are set to 100%. This approach does not allow for diversification and is known to yield usually very conservative economic capital.

4.2 A more general approach: Copulas

5-variate copulas

We first fit 5-variate copulas to the transformed data set of market risks, among the Gaussian, the Student, the Gumbel and the Clayton. We restrict to these copulas since the two first ones are widely used in the industry and because all verify the constraints allowing to perform the calculation. We then estimate the parameters of these copulas and compute the corresponding maximized log-likelihood function and BIC criterion. These values are provided in Table 12. From these findings, the best fitting copula is given by the Student-t with ν * = 9 degrees of freedom. Note also that, in that case, the elliptic copulas show a better fit than the Archimedean copulas because the former present more parameters to describe the dependence structure than the latter ones.

FNAC models

In a second time, we apply the FNAC models. Once again, the Gumbel and the Clayton copulas are used since these are the only Archimedean copulas that one can employ to model FNAC. Note that in this approach, the dependence structure between our five market risks is described through four parameters instead of only one as in the previous methodology. This really improves the fit of the Archimedean copulas. Regarding the structure, we first arrange the variables according to their dependence (Interest Rate, Spread, Equity, FX and Implied Volatility risks) then in terms of economic capital (Equity, Spread, Implied Volatility, Interest rates and FX). However, both structures do not satisfy the constraint of decreasing dependence. Therefore, this approach will not be used for the purpose of economic capital calculation.

C-Vines

Kendall In order to use the vine methodology, we first compute the samples measures of Kendall's tau and Spearman's rho. Both measures give the same results to determine the C-vine structure we to use to describe the dependence structure between our five market risks. The values of these samples measures are given in Tables 14 and15. The variable that is the most correlated with the other variables is taken as the main node of the tree. In our case, Interest rate (IR) risk is the key variable of the first tree. Then comes as the main node Spread risk conditionally on IR risk in the second tree, IV risk conditionally on IR and Spread risks in the third tree and finally Equity risk conditionally on IR, Spread and IV risks in the last tree.

In other words, we have in the first layer, the dependence between the variable "IR" and all the other variables in this system is modeled with bivariate copulas. The second layer consists in modeling the dependence of variables Spread with variables IV , Equity and F X, conditionally on the variable IR. In the last layer, one uses a bivariate copula to model the dependence between the variables Equity and F X, conditionally on the variables IR, IV and Spread. Finally, we estimate this dependence structure. In Table 16, we provide the values of maximized loglikelihood function and BIC criterion for the different C-Vines copula models.

Based on these figures, it appears that the Student C-vine copula model is the model that best matches our market data set. It is important to keep in mind that during all this exercise, we estimate non-parametrically the margins by using the empirical distributions.

From our findings, one can say that the multivariate Student-t and the Student-t C-vine model appear to be the best approaches of fitting copulas. Then, using pair-copula construction models really improves the fit of the Gumbel and Clayton copulas. Note that the C-vines show a better fit than the 5-variates copula models. Eventually, all the above copula models are used to estimate the economic capital that is needed to cover losses at a certain risk tolerance level through the calculation of the VaR(99.5%) and the ES(99%).


Allen L Husker

Vladimir Kostoglodov

Victor M Cruz-Atienza

Denis Legrand

Nikolai M Shapiro

Juan S Payero

Michel Campillo

Eduardo Huesca-Pérez

Temporal variations of non-volcanic tremor (NVT) locations in the Mexican subduction zone: Finding the NVT

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

Introduction

Non-volcanic tremor (NVT) in the Japan, Cascadia, and Costa Rican subduction zones typically occurs when the slab is at the depth of 30 -50 km [e.g. [START_REF] Obara | Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[END_REF][START_REF] Brown | Deep low-frequency earthquakes in tremor localize to the plate interface in multiple subduction zones[END_REF][START_REF] Ghosh | Tremor patches at Cascadia revealed by array analysis[END_REF][START_REF] Kao | Northern Cascadia episodic tremor and slip: A decade of tremor observations from 1997 to 2007[END_REF]. Recent precise detections of Low Frequency Earthquakes (LFE), which make up NVT, locate the LFE's on the interface between the subducting oceanic plate and overriding mantle wedge [START_REF] Brown | Deep low-frequency earthquakes in tremor localize to the plate interface in multiple subduction zones[END_REF][START_REF] Kato | Variations of fluid pressure within the subducting oceanic crust and slow earthquakes[END_REF]. In Guerrero, Mexico, the slab descends from the Middle American Trench to ~40 km depth at about 150 km from the trench. It then remains horizontal near 40 km depth for ~140 km [START_REF] Pérez-Campos | Horizontal subduction and truncation of the Cocos Plate beneath central Mexico[END_REF][START_REF] Kim | Geometry and seismic properties of the subducting Cocos plate in central Mexico[END_REF] (Figure 1). At ~290 km from the trench the slab plunges deeper into the mantle [START_REF] Pérez-Campos | Horizontal subduction and truncation of the Cocos Plate beneath central Mexico[END_REF]Husker and Davis, 2009;[START_REF] Kim | Geometry and seismic properties of the subducting Cocos plate in central Mexico[END_REF]. NVT in Mexico occur in a 75 km long zone above the 140 km horizontal slab, far from the mantle wedge (Figure 1) [START_REF] Payero | Nonvolcanic tremor observed in the Mexican subduction zone[END_REF]. As a result of the slab geometry, the commonly used terms 'updip' and 'downdip' can be misleading in the sub-horizontal section of the slab where there is no dip, but we use the terms for clarity when describing location along the slab. SSE has been observed in two types: (1) short-term (small repeat interval), small magnitude SSE and (2) long-term (long repeat interval), large magnitude SSE [e.g. [START_REF] Hirose | Repeating short-and long-term slow slip events with deep tremor activity around the Bungo channel region, southwest Japan[END_REF]Obara and Hirose, 2006;[START_REF] Obara | Phenomenology of deep slow earthquake family in southwest Japan: Spatiotemporal characteristics and segmentation[END_REF][START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF]. The short-term SSE and NVT occur simultaneously in space and time. On the other hand, the long-term SSE occurs updip from the source region of NVT and short-term SSE. The episodic tremor and slip (ETS) observed in Cascadia is then the short-term SSE, while the 2006 Mexican SSE observed to be updip from the NVT belongs to the long-term SSE [START_REF] Rogers | Episodic Tremor and Slip on the Cascadia Subduction Zone: The Chatter of Silent Slip[END_REF][START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF][START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF].

Previous studies have linked evidence of high pore fluid pressure to subduction zones where NVT occurs [e.g. [START_REF] Obara | Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[END_REF][START_REF] Shelly | Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip[END_REF][START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF]. In Mexico, studies have presented a number of possible sources of fluid in the crust above the flat slab. These include the release of trapped water from the slab after passing overpressurized impermeable boundaries [START_REF] Song | Subducting Slab Ultra-Slow Velocity Layer Coincident with Silent Earthquakes in Southern Mexico[END_REF][START_REF] Kim | Geometry and seismic properties of the subducting Cocos plate in central Mexico[END_REF], water releasing metamorphic phase changes [START_REF] Manea | Thermal structure, coupling and metamorphism in the Mexican subduction zone beneath Guerrero[END_REF][START_REF] Manea | Flat-Slab Thermal Structure and Evolution Beneath Central Mexico, Pure and App[END_REF], and water squeezed from the slab due to slab bending [START_REF] Jödicke | Fluid release from the subducted Cocos plate and partial melting of the crust deduced from magnetotelluric studies in southern Mexico: Implications for the generation of volcanism and subduction dynamics[END_REF]. Evidence from a magnetotelluric study [START_REF] Jödicke | Fluid release from the subducted Cocos plate and partial melting of the crust deduced from magnetotelluric studies in southern Mexico: Implications for the generation of volcanism and subduction dynamics[END_REF]] and a seismic velocity tomography [Huesca-Pérez and Husker, submitted] suggests that water is present throughout the entire lower crust above the flat slab and that the amount of fluid increases in the downdip direction. The evidence from the magnetotelluric study and tomography do not distinguish which of the fluid sources is most probable and it may be a combination of the three.

The seismic network that measured the NVT in Mexico during the 2006 SSE was a very dense (~6 km spacing) temporary line of broadband seismometers installed perpendicular to the trench from 2005 -2007 called the MesoAmerican Subduction Experiment (MASE) [http://www.tectonics.caltech.edu/mase/] (Figure 1). The original NVT study [START_REF] Payero | Nonvolcanic tremor observed in the Mexican subduction zone[END_REF] of data from MASE used cross-correlations of envelopes of NVT [START_REF] Obara | Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[END_REF] in order to locate the tremor. The method suffers from imprecise locations because correlating emergent windows over many minutes or hours causes inexact, rolling cross-correlation peaks. In addition, the locations were limited because the MASE station arrangement was approximately a one-dimensional line. In order to bypass these limitations, [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF] plotted the NVT energy (square of velocity) of the north component recorded at each station during the two and a half year duration of MASE to map the region of NVT along the MASE line.

In this study, we determine epicenters of all the individual NVT's during that period using an energy location method. We then analyze the epicentral locations of individual NVT's over time. When there are a large number of individual NVT's over an interval of a few days, this is considered an episode of NVT as defined by other authors [e.g. [START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF][START_REF] Obara | Phenomenology of deep slow earthquake family in southwest Japan: Spatiotemporal characteristics and segmentation[END_REF]. [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF] exclusively analyzed episodes of activity, although they were described as bursts since the individual NVT's that make up an episode were not seen in that study. We then analyze the energy of the NVT in different regions and its association with long term and short term SSE's in Mexico.

NVT Location Method

The location technique employed in this study is similar to amplitude location methods in volcano studies [e.g. [START_REF] Legrand | Moment tensor analysis of near-field broadband waveforms observed at Aso Volcano, Japan[END_REF][START_REF] Battaglia | Location of seismic events and eruptive fissures on the Piton de la Fournaise volcano using seismic amplitudes[END_REF][START_REF] Taisne | Imaging the dynamics of magma propagation using radiated seismic intensity[END_REF]. The difference between the amplitude location method with volcanic tremor and with NVT is that volcanic tremor is largely isotropic with P waves dominating while NVT appears to have a double couple source with S waves dominating [START_REF] Obara | Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[END_REF][START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF]. In our data, the energy of the north component is systematically greater than the east component, which is greater than the vertical component with the location of the energy maximums and the energy falloff corresponding between all components (Figure 2). It is very difficult to generate radiation patterns with corresponding shapes between all components from an ideal single double couple source, as maximums in the radiation pattern for S waves are nodes for P waves and vice versa. There is evidence that NVT is made up of many small low frequency earthquakes (LFE's) [START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF]. Multiple small horizontal approximately east-west striking double-couple sources with rakes pointing northward generated over a small area could create the observed energy pattern. This NVT mechanism model is supported by the general north-south Cocos slab and North American plate convergence in this region [START_REF] Demets | Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions[END_REF]. Scattering plays a role in diminishing the observed radiation pattern as well, but it obviously has not completely eliminated the radiation pattern because the energy level is different in each of the three components. Assuming this dislocation model for the observed data, the energy maximum at the surface is the location of the epicenter. There may be small deviations based on small differences in the focal mechanism, particularly the dip, but this effect must be small or the observed energy shapes would substantially differ between components.

Onsets and durations of NVT were automatically determined using the energy of median filtered envelopes of seismograms [START_REF] Husker | Automatic non-volcanic tremor detection in the Mexican subduction zone[END_REF]]. An individual NVT was established when the energy averaged across all stations surpassed an empirically determined threshold. The energy was then summed into two-week intervals and presented over time (2005 -2007 the duration of MASE) [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF]. The background of Figure 3 was created in the same way. This gave the broad details of NVT location along the MASE profile during this period and showed a significant increase of NVT during the 2006 SSE [START_REF] Payero | Nonvolcanic tremor observed in the Mexican subduction zone[END_REF][START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF][START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. However, the summing of NVT energy in previous studies only showed the dominant NVT trend and not the complete spatio-temporal variations of smaller individual NVT's.

In order to determine the temporal and spatial variations of the NVT, each tremor was located using the energy falloff from the source with distance due to geometrical spreading and intrinsic attenuation (Equation 1). The goal of the inversion technique to locate the NVTs was to find the epicenter of each event along the MASE array. The onedimensional shape of the seismological network does not allow for solving for a focal mechanism, the area, or the depth of the NVT, as there are too many degrees of freedom that trade-off with one another. We only determine the epicenters using the following simple equation:

, ( 1
)
where E is the energy (summed squared velocity) of the NVT measured at each station, C is a constant which is related to the energy at the source, r is the distance between source and receiver, is the peak frequency, v is the average seismic velocity of the medium, and Q is the quality factor. The latitude and longitude of the epicenter and C were the parameters solved to minimize the least-square error between the observed and theoretical energy profiles along the array [e.g. [START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF]. Only the energy, E, from the northern component was used in the inversion as it had the greatest signal to noise ratio. Site effects were removed before inverting for the location using coda-determined site factors [START_REF] Husker | Automatic non-volcanic tremor detection in the Mexican subduction zone[END_REF]. ycu"ugv"vq"dg"307"J¦"*z"4ヾ+"ukpeg"vjg"ugkuoqitcou"ygtg"hknvgtgf"3"-2 Hz due to the range of detectable NVT frequencies (~1 -10 Hz) and ambient noise (> 2 Hz).

NVT's have a more or less fixed spectrum, unlike earthquakes, so a small frequency band is representative of the overall energy [START_REF] Aguiar | Moment release rate of Cascadia tremor constrained by GPS[END_REF][START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF]. The quality factor Q was set to 276 determined for this frequency range from attenuation studies in Mexico [START_REF] Garcia | Inslab Earthquakes of Central Mexico: Q, Source Spectra, and Stress Drop[END_REF]. The depth was conditionally set to 40 km as many recent studies have found all or most NVT at or just above the plate interface [START_REF] Shelly | Non-volcanic tremor and low-frequency earthquake swarms[END_REF][START_REF] Brown | Deep low-frequency earthquakes in tremor localize to the plate interface in multiple subduction zones[END_REF]. However these studies use limited time windows that may create limitations of allowed moveouts forcing only detections of deeper sources. Kao et al.

[2009] found sources throughout the crust in Cascadia with concentrations of NVT 10 -20 km above the plate interface. Therefore, a range of depths was implemented in our inversion, but the results changed very little because it is just a search for the maximum of the energy parabola as already mentioned. Initial values for latitude, longitude, and C were chosen from the station that recorded the maximum energy. Location results with standard deviations greater than 1% of the final value of the latitude or the longitude were thrown out. Figure 4 depicts two examples of the inversion results.

NVT Locations over time

The NVT epicenter location results are shown in Figure 1 and largely agree with previous location studies [START_REF] Payero | Nonvolcanic tremor observed in the Mexican subduction zone[END_REF][START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF]. Figure 3 the trench, which we label the Sweet Spot. The term "Sweet Spot" is most commonly used in baseball and acoustics to refer to a limited zone with the correct conditions to create a particular effect (e.g. http://en.wikipedia.org/wiki/Sweet_spot), and it is used here with the same understanding. There is inter-episode NVT on the downdip side of the NVT activity in Cascadia as well, but it is not always in the same location along strike of the trench [START_REF] Wech | Automated detection and location of Cascadia tremor[END_REF]. In southwest Japan, the studies there focus on NVT concentrations along strike, rather than perpendicular to the trench [e.g. [START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF].

There are Sweet Spots observed there as well where certain points along strike exhibit continuous inter-episode NVT, but it is not clear if they are on the downdip side of the NVT zone as in Cascadia and Mexico [START_REF] Obara | Phenomenology of deep slow earthquake family in southwest Japan: Spatiotemporal characteristics and segmentation[END_REF].

In Lower energy NVT epicenters were also detected further updip with their own separate NVT episodes. Figure 3 depicts clusters of many NVT's during the 2006 SSE that are updip from the high-energy episodes. Their summed NVT energy is much lower than the downdip epicenters, and therefore they were not observed by [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF] and do not appear as episodes in Figure 3. In order to observe the summed energy of the updip tremors exclusively, only the energy profiles from NVT epicenters < 190 km from the trench are summed in the background of Figure 5. The updip, low energy NVT epicenters (green dots, Figure 5), which have total energy about four times lower than the dowdip NVTs, start to occur about 1 month after the beginning of the SSE, and the strongest updip NVT's do not occur until at least 3 months after the initiation of the SSE.

NVT durations have been found to vary from seconds to days, but the amplitude of tremor remains nearly constant [e.g. [START_REF] Aguiar | Moment release rate of Cascadia tremor constrained by GPS[END_REF]. Thus, the energy of NVT is nearly linearly correlated with the duration of the individual NVT. Although we find that the energy (C) and the duration ( v) are correlated (Figure 6a), NVT in the high-energy zone (> 190 km from the trench) also exhibits greater amplitudes than the NVT in the low energy zone (< 190 km from the trench) (Figure 6b). The duration-normalized energy (C/ v, also called the energy release rate) of NVT in Figure 6 would be roughly constant over time if C were always proportional to v. Instead, the updip NVT's are limited to ~0.4 of the largest downdip NVT's (Figure 6b) implying limited amplitudes as well. In addition, apart from a few events in June-July 2005 (which may acutally be a small episode) the NVT activity between episodes during the whole time span is upper-bounded by the ~0.4 level below which the updip, low-energy NVTs are confined (Figure 6b). This shows that during the episodes within the high-energy zone, the NVT energy release rate increases, contrary to what has been observed in Cascadia, where the average rate does not fluctuate [START_REF] Aguiar | Moment release rate of Cascadia tremor constrained by GPS[END_REF]. [START_REF] Wech | A continuum of stress, strength and slip in the Cascadia subduction zone[END_REF] also found that there is greater NVT energy downdip in Cascadia, but this is purely due to duration of the NVT episodes and the quantity of NVT's in those episodes.

The majority of the 2006 SSE dislocation extended to approximately 165 km from the trench [START_REF] Larson | The 2006 aseismic slow slip event in Guerrero, Mexico: New results from GPS[END_REF][START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF][START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. A more recent inversion found small amounts of slip as far as 215 km from the trench, but the maximum was always updip from the NVT Sweet Spot [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF] (Figure 7).

NVT did not occur within the peak of the slipping front of the SSE and primarily only overlapped the downdip edge of it (Figure 7) [START_REF] Larson | The 2006 aseismic slow slip event in Guerrero, Mexico: New results from GPS[END_REF][START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF][START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF][START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. This confirms, along the whole NVT region, the spatio-temporal separation between NVT and SSE previously observation by [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF] over the NVT high-energy segment. However, the updip, low energy NVT occurs concurrently with the 2006 SSE and appears to be a direct effect of the SSE although the low energy NVT is also limited spatially (> 150 km from the trench). The updip limit of the NVT is very close to the corner of the slab (Figures 1,7).

Discussion

NVT episodes provide evidence for short term SSE

The observation of spatially separated SSE and NVT may occur because there are two types of slow slip events as seen in the southwest Japan subduction zone [START_REF] Obara | Phenomenology of deep slow earthquake family in southwest Japan: Spatiotemporal characteristics and segmentation[END_REF][START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF]. There are short-term (small repeat interval), small magnitude SSE that occur downdip from the long-term (long repeat interval), large magnitude SSE [START_REF] Hirose | Repeating short-and long-term slow slip events with deep tremor activity around the Bungo channel region, southwest Japan[END_REF][START_REF] Obara | Phenomenology of deep slow earthquake family in southwest Japan: Spatiotemporal characteristics and segmentation[END_REF][START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF]. In the case of Japan, both types of SSE are accompanied by NVT, however more NVT occur during the short term SSE. Thus, Obara

[2011] categorized the Cascadia episodic tremor and slip (ETS) as short term SSE. Due to the previously observed separation of NVT's from the SSE's in Mexico [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF], those SSE's were classified as long term SSE's [START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF].

Varios studies have suggested that with the absence of detailed geodetic measurements, NVT episodes could be used to identify short term SSE [START_REF] Aguiar | Moment release rate of Cascadia tremor constrained by GPS[END_REF][START_REF] Wech | Automated detection and location of Cascadia tremor[END_REF][START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF]. The observation of NVT episodes downdip from the long term SSE zone in Mexico provides evidence for short term SSE. In addition, possible short term SSE's have been detected within the NVT region [START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. Detections were so small within the GPS signal that it was impossible to determine the precise location, but the strongest signals were all located well within the NVT region [START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. With each detectable short term SSE a high energy NVT episode occurred [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF][START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. High-energy NVT episodes with no clear short term SSE were also observed, but short term SSE detection is problematic because the corresponding displacements recorded by GPS are almost at the noise level of the position time series [START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. Thus, the model seen in Japan of cyclical short term SSE's with strong NVT episodes and long term SSE's with weak NVT episodes follows the evidence seen in Mexico [START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF]. The increased rate of the NVT episodes during the 2006 long term SSE suggests that the short term SSE's are triggered by the long term SSE.

Conditions for NVT

The long term SSE region in Mexico aligns with an ultra slow velocity layer (USL)

found between the slab and overriding plate that is probably a remnant mantle wedge from before the slab was horizontal [START_REF] Pérez-Campos | Horizontal subduction and truncation of the Cocos Plate beneath central Mexico[END_REF][START_REF] Song | Subducting Slab Ultra-Slow Velocity Layer Coincident with Silent Earthquakes in Southern Mexico[END_REF]. The high energy NVT region, which is now considered due to a short term SSE, aligns with a moderate USL [START_REF] Song | Subducting Slab Ultra-Slow Velocity Layer Coincident with Silent Earthquakes in Southern Mexico[END_REF]. The USL's are evidence of high pore fluid pressure that allows the two types of SSE's to occur. These observations coincide with those by [START_REF] Kato | Variations of fluid pressure within the subducting oceanic crust and slow earthquakes[END_REF] where the long term SSE's occur in a highly pressurized region, while short term SSE's and their associated NVT's occur in a modestly pressurized region. NVT's in Mexico do occur directly because of the long term SSE, but they do not extend throughout the entire long term SSE region (Figures 1, 3, 5, 7), and there is no evidence in Mexico of how the short term SSE region aligns with the down dip NVT. Thus, the limited NVT zone within the long term SSE is evidence that NVT might require modestly high pore fluid pressure.

This may suggest that the Sweet Spot is the zone with the proper pressure and/or temperature for NVT to persistently occur and that NVT that occur outside of the Sweet Spot may require an additional input. Evidence of a pressure change was observed during the 2006 SSE when the average seismic velocity within NVT region in the crust was observed to decrease, indicating increased dilation [START_REF] Rivet | Seismic evidence of nonlinear crustal deformation during a large slow slip event in Mexico[END_REF]. As mentioned, during the 2006 SSE, NVT was observed farther from the Sweet Spot than at any other time during this study coinciding with the observed dilation increase.

Conclusions

This study found that persistent background NVT occurs approximately 215 km from the trench in a flat section of the slab in Guerrero, Mexico. High-energy NVT tremor episodes occur every few months with durations of a few weeks. [START_REF] Hirose | Repeating short-and long-term slow slip events with deep tremor activity around the Bungo channel region, southwest Japan[END_REF][START_REF] Obara | Phenomenology of deep slow earthquake family in southwest Japan: Spatiotemporal characteristics and segmentation[END_REF][START_REF] Obara | Characteristics and interactions between non-volcanic tremor and related slow earthquakes in the Nankai subduction zone, southwest Japan[END_REF]. Figure 3. The NVT energy and epicenter locations. In the top panel the background color is the NVT seismic energy measured at the surface summed in two-week intervals as described by [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF]. This figure differs slightly from Figure 2B of [START_REF] Kostoglodov | Slow Slip Event and Nonvolcanic Tremor in the Mexican Subduction Zone[END_REF] because energy profiles that did not meet the error criteria described in the localization method were considered spurious, non-NVT events and not used in the sum. The dots are the epicenters of the NVT from the inversion. The dots sometimes do not align exactly with the background energy profile due to how the summed energy is binned into 15-day intervals. The green dots are the trenchward, updip, lowenergy epicenters. The orange dots are the downdip epicenters with mostly high energies as noted by the background colors. The gray dots are the inter-episode epicenters. The y-axis is the time and the x-axis is distance from the trench. The cross-section of the continental crust and slab in the lower panel (profile A -A' from Figure 1) aligns with the upper panel in order to demonstrate the details from the electrical conductivity measured in the crust [START_REF] Jödicke | Fluid release from the subducted Cocos plate and partial melting of the crust deduced from magnetotelluric studies in southern Mexico: Implications for the generation of volcanism and subduction dynamics[END_REF]. The blue color in the cross-section represents high resistivity and low conductivity, and red is the opposite. Temperature contours in degrees Celsius from [START_REF] Manea | Thermal structure, coupling and metamorphism in the Mexican subduction zone beneath Guerrero[END_REF] are also shown in the cross-section. storms as seen in these two events. Although noisy storm data from the coastal stations was included in the inversion, the effect to the epicentral locations was quite small (ø"3"mo+0"

When the coastal stations were removed from the inversion the locations were 182.7 km ± 3.2 km and 214.9 km ± 12.2 km respectively. 1) related to the NVT duration, v, for the downdip high-energy NVT (orange dots) and the updip low-energy NVT (green dots). The gray dots are the inter-episode NVT's. These are the same events and color scheme shown in Figures 1 and3. C is directly related to the seismic energy of the NVT at the source. shows that C is roughly proportional to v as has been observed in other studies [START_REF] Aguiar | Moment release rate of Cascadia tremor constrained by GPS[END_REF], however the relation is very disperse. shows the value of v over time.

  and the NVT.mp4 animation in the auxiliary material show how the NVT epicenter locations change over time. Figure 3 depicts almost continuously occurring NVT near 215 km from

  addition to the continuous NVT activity at the Sweet Spot, there are NVT episodes occurring every few months. The episodes, which last a few days to weeks, are composed of many NVT's. The background of Figure3is the sum of all the individual NVT's over two week intervals. Thus, higher sums (E > 40 X 10 -12 m 2 /s; red background) can be due to more individual NVT's and/or higher energy tremors. The coincidence of large numbers of NVT's and the high summed energy (red background) are considered high-energy episodes. It is difficult to determine if there is periodicity of the episodes as only 3 intervals between episodes occur outside of the 2006 SSE. Throughout the 2006 SSE, the frequency and duration of the NVT episodes increased. During the episodes, NVT epicenters extend updip from the Sweet Spot, near 215 km from the trench (Figure3). The NVT epicenters in the episodes are observed to migrate or progress updip at a rate of ~1 -20 km / day similar to migration observed in Cascadia[START_REF] Wech | A continuum of stress, strength and slip in the Cascadia subduction zone[END_REF]. During the 2006 SSE it was not possible to observe migration due to much higher activity throughout the NVT region.

Figure Captions Figure 1 .

 Captions1 Figure Captions

Figure 2 .

 2 Figure 2. Example tremor seismogram. This is from an NVT episode recorded at the station

Figure 4 .

 4 Figure 4. Example results from the inversion. Energy in the figure was normalized to the

Figure 5 .

 5 Figure 5. The seismic energy (background colors) measured at the surface from only the

Figure 6 .

 6 Figure 6. Normalized value of C from the inversion (Equation1) related to the NVT

Figure 7 .

 7 Figure 7. NVT epicenters from this article and the 2006 SSE slipping front from Radiguet

  The recurrence rate of those episodes increased during the 2006 SSE. During that same SSE, low energy NVT were observed closer to the trench than at any other time, between 150 and 190 km from the trench. This NVT activity was about four times less energetic than further down-dip and had a lower energy-release-rate. The high-energy NVT zone, between 190 to 220 km from the trench, coincides with a region of modestly high pore fluid pressure in between the slab and overriding plate. The spot 215 km from the trench has the necessary conditions to generate NVT persistently and so we refer to it as the NVT Sweet Spot. All NVT generated out of the Sweet Spot require an additional input to generate NVT. We suggest that the additional stress and increased dilation for fluid content necessary to generate NVT outside of the Sweet Spot are the SSE's. Assuming that each high-energy NVT episode therefore indicates a short term SSE, the increase in the number of NVT episodes during the long term SSE suggests that short term SSE's can be triggered by SSE's. These observations and interpretation of long and short term SSE coincide with observations in Japan where SSE's are found to occur closer to the trench than short term SSE's, the two interact, and NVT's coincide with short term SSE's [Hirose

Acknowledgments

We would like to thank Guillermo González for preparing the NVT.mp4 video in the supplementary material. This study was supported by PAPIIT IN110611, CONACYT