Malik Kamel

Bensafta
email: bensafta.malik@yahoo.fr

Kamel Malik Bensafta

« A Regional Analysis of Markets Uncertainty Spillover »

Keywords: JEL classification: C5, E5, F3, G1 MGARCH, Structural breaks, subprime crisis, volatility transmission MGARCH, Ruptures Structurelles, Crise des Subprimes, transmission de volatilité. Classification JEL : C5, E5, F3, G1

This study aims to describe the transmission of uncertainty between the stock markets of four aggregate regions: North America, Europe non Euro-zone, Asia and the euro area. We use a non-linear VAR model with innovations following a Multivariate GARCH with variance regime change. The interest of the model with regime change is to correct the estimation bias caused by the overestimation of the shocks persistence. We apply the non-linear VAR model with regime change in daily MSCI data aggregated from four regions over the period from June 2005 to October 2013. This period included the crisis episodes in 2007 and 2011. Our results indicate the importance of taking into account changes in variance in measuring the persistence of volatility shocks. They also show the high exposure of European and Asian markets to the uncertainties of North American markets. The transmission in time of crisis is higher compared to the quiet period. This result confirms the contagious nature of the crises of 2007 and 2011 and supports the thesis of the contingency theory to crisis.

Introduction

A number of empirical studies look at the effects of contagion that were observed following the American stock exchange crash of October 1987 and the frequent crises of the emergent countries during 1990's (King and Wadhawani, 1990and Edwards, 1998[START_REF] Forbes | No contagion, only interdependence: measuring stock market comovements[END_REF][START_REF] Mcaleer | Testing for Contagion in Asian Exchange Rates[END_REF]. Studying spillover mechanisms and market correlations leads to two conclusions: first, market comovement is important in any strategy of international diversification [START_REF] King | Volatility and Links between National Stock Markets[END_REF]. This international diversification is based on a weak cross-correlation of geographically spaced markets. Secondly, in the current context of financial globalization, mean and variance spillover reinforce market interdependencies and reduce dramatically any benefit of diversification strategies.

International transmission mechanisms were established between stock markets and well demonstrated by studying volatility [START_REF] Ng | Volatility Spill-over Effects from Japan and the U.S. to the Pacific-Basin[END_REF][START_REF] Granger | A Bivariate Causality Between Stock Prices and Exchange Rates: Evidence from Recent Asia Flu[END_REF]. The volatility transmission translates the exogenous part of market turbulence linked to other market-uncertainties [START_REF] Engle | Meteor Showers or Heat Waves? Heteroskedastic Intra-Daily Volatility in the Foreign Exchange Market[END_REF] 2 The GARCH models are widely used for modelling the volatility of financial series. Under GARCH process, shocks to volatility persist according to ARMA process of squared innovations. Empirical findings show strong persistence of high frequency financial series and this is usually near unity . According to this definition, the dependence in variance is a sign of market imperfections and allows risks and returns predictability. It appears that those markets are increasingly dependent in variance [START_REF] Hamao | Correlations in Price Changes and Volatility across International Stock Markets[END_REF][START_REF] Koutmos | A symmetric Volatility Transmission in International Stock Markets[END_REF], since there is more information in market volatility than market prices [START_REF] Kyle | Continuous Auctions and Insider Trading[END_REF].

3

In this paper, the model of market volatility is a multivariate GARCH process. Short-run Mean and variance spillovers are based on Granger Causality and variance Causality. First, we look for bias estimates in volatility spillover estimates in the standard BEKK model [START_REF] Engle | Multivariate Simultaneous Generalized ARCH[END_REF] compared to the BEKK model with a Structural Break in Variance subsequently called BEKK-BSV [START_REF] Bensafta | Chocs, Chocs de Volatilité et Contagion entre les Marchés Boursiers : Application d"un modèle ICSS-MGARCH[END_REF]. In the last section, we provide ways to extend this work. . However, [START_REF] Lamoureux | Persistence in Variance, Structural Change, and the GARCH Model[END_REF] show that misspecification in conditional variance processes explain higher persistence measurement. They found that time varying coefficients may exhibit persistence and they proposed time variation of unconditional variance. When taking into account structural change in unconditional variance, they obtained mode reduces persistence value. Theoretically it is hard to detect such structural change but there exist various methods to detect structural change empirically such as regime switching models [START_REF] Susmel | Switching volatility in private international equity markets[END_REF]. In this paper we consider a method based on the CUSUM test, namely the Iterated Cumulative Sum of Squared ICSS algorithm developed by Inclan and Tiao (1994) and Sanso et al. (2004). 2 Engle et al (1990) described volatility spillover as "meteor showers".

3 Shock persistence is the measurement of cumulative effects of shocks on volatility. For a GJR-GARCH (p, q) process, persistence is equal to 𝜑 𝑥 = ∑ �𝛼 𝑖𝑥 + 1 2 𝛾 𝑖𝑥 � 𝑝 𝑖=1

+ ∑ 𝛽 𝑗𝑥 𝑞 𝑗=1

. The measurement of unconditional variance is (ℎ 𝑡 ) = 𝜔 1-𝜑

. A higher unconditional variance leads to highly persistence estimates.

The Econometric model

Volatility is the fluctuation of the value of an unspecified asset around a central tendency; it can also be measured as variation asset return over time. Negative shocks have more impact on volatility than positive ones of the same magnitude say asymmetric effect [START_REF] Engle | Autoregressive Conditional Heteroskedasticity with estimates of the variance of United Kingdom inflation[END_REF]NG, 1993, Koutmos and[START_REF] Koutmos | A symmetric Volatility Transmission in International Stock Markets[END_REF]. These interactions give useful information about market dynamics, such as the value of financial assets, market indexes, and the foreign exchange rate. Since [START_REF] Engle | Autoregressive Conditional Heteroskedasticity with estimates of the variance of United Kingdom inflation[END_REF] first developed the ARCH process, generalized by [START_REF] Bollerslev | Generalized Autoregressive Conditional Heteroskedasticity[END_REF], GARCH models have proved their ability to capture properties of high frequency financial series such as volatility clustering and thick tails. Volatility clustering indicates a persistence of volatility shocks. For reasons of parsimony and simplicity, the GARCH(1,1) model could suffice. Non-linear univariate GARCH models meet the objective of estimating volatility and the management characteristics of financial series such as volatility clustering. However, univariate models contain some important biases: omitted variable bias, overestimating the persistence of volatility shocks, and quadratic specification of the conditional variance masks the asymmetric effects of shocks. Moreover, they do not show volatility transmission. Multivariate GARCH models such as BEKK-GARCH and Vech meet the latter objective, but still overestimate persistence4 . This bias is particularly responsible for the overvaluation of the conditional correlations between markets. Because of this bias, these standard models underestimate the proportion of exogenous volatility related to external events. [START_REF] Ewing | Re-examining the Asymmetric Predictability of Conditional Variances: The Role of Sudden Changes in Variance[END_REF], and [START_REF] Bensafta | Non-Stationary Variance and Volatility Causality[END_REF]Semedo (2009, 2011) used a MGARCH-BEKK model with structural breaks in variances to study the cross-market volatility transmission and contagion. [START_REF] Ng | Volatility Spill-over Effects from Japan and the U.S. to the Pacific-Basin[END_REF] and [START_REF] Bensafta | Non-Stationary Variance and Volatility Causality[END_REF]Semedo's model (2009, 2011), confirms the leading place of the U.S. market worldwide. The introduction of structural break in variance corrects the overestimation of the persistence of volatility shocks [START_REF] Lamoureux | Persistence in Variance, Structural Change, and the GARCH Model[END_REF]. This correction shows that crises are not always contagious, but markets should be cautious of crises originating in the USA, supporting the results of "No Contagion only Interdependence" from [START_REF] Forbes | No contagion, only interdependence: measuring stock market comovements[END_REF].

Econometric model of market volatility

In this paper, we used a quadrivariate MGARCH model with structural changes in the variance. We used a quadrivariate version to introduce regional as well as overall effects of the American aggregate stock markets. Intuitively, the introduction of the structural change in variance should improve the standard BEKK models in terms of log-likelihood, while being closer to the real model.

The market volatility model used in this study is a MGARCH. First of all, we produced a BEKK-diagonal asymmetric model [START_REF] Engle | Multivariate Simultaneous Generalized ARCH[END_REF]. This standard model was then compared with our model: a BEKK model increased by structural change in variance.

Let 𝑟 𝑡 = ( 𝑟 1𝑡 … 𝑟 𝑀𝑡) 𝑡𝑟 denote the vector of logarithmic returns of MSCI indices, and 𝑢 𝑡 = ( 𝑢 1𝑡 … 𝑢 𝑀𝑡) 𝑡𝑟 the vector of dynamic VAR(n1) residual, such that:

Because the data are daily, we included the possibility of holiday's effects and weekends effects. We also included central bank disclosures, the conduct of monetary policy, and surprise effects. In this multivariate model, mean transmission is measured by the VAR coefficients of conditional mean equations.

Φ(𝐿)(𝑟

𝑡 -𝜇) = 𝑢 𝑡 (1)
Where Ф (L) is the function with lags in the VAR (n 1 )5 Φ(𝐿)(𝑟 𝑡 -𝜇) = 𝐶(𝐿)𝒢(𝑋 𝑡 ) + 𝑢 𝑡 process. Mean equations are augmented with exogenous regressors such as monetary policy variables. This general case is:

(2)

Where 𝒢(𝑋 𝑡 ) is exogenous variables vector and 𝐶(𝐿) the exogenous coefficient matrix. The mean crossmarket transmission is described by Φ ij k coefficients of the VAR process. Suppose that 𝑢 𝑡 is a vector of non-autocorrelated VAR residual, and:

𝑢 𝑡 = 𝐻 𝑡 -1/2 𝜀 𝑡 (3)
Where ε t is an N-dimension vector of white noise elements, such that ε t ↝ i. [START_REF] Tse | A Test for Constant Correlations in a Multivariate GARCH Model[END_REF]Tsui, 2002, and[START_REF] Engle | Autoregressive Conditional Heteroskedasticity with estimates of the variance of United Kingdom inflation[END_REF]Shepphard, 2001). [START_REF] Bauwens | Multivariate GARCH Models: A Survey[END_REF] provided an extensive literature review of the MGARCH model (see appendix A. for details). The BEKK-diagonalasymmetrical model with transmission in variance and structural break is defined as follows:

(4) 𝐶 𝐵 , 𝐴 𝑖 , 𝐵 𝑗 𝑎𝑛𝑑 𝐺 𝑗 are coefficients matrices in conditional variance-covariance equations of standard model. 𝑢 𝑡 is innovation vector, 𝐷 𝑀𝑡-𝑖 is M dimension vector of dummy variables 7 𝐶 𝐵 is a constant coefficient matrix in conditional variance-covariance equations. . In the diagonal, asymmetric model with structural changes in variance and variance spillover (BEKK-BSV model (𝑑 𝑖𝑡 ) where:

�𝑑 𝑖,𝑡 � = � 1 𝑤ℎ𝑒𝑛 𝑢 𝑖𝑡 < 0 0, 𝑒𝑙𝑠𝑒 6
thereafter), the (𝐶 𝐵 ) 𝑡𝑟 (𝐶 𝐵 ) take into account structural changes in variance. Each diagonal element of 𝐶 𝐵 is defined as follows:

�𝐶 𝐵𝑗𝑗 � = 𝜔 0𝑗 + ∑ 𝜔 𝑗𝑖 𝑆 𝑗𝑡 𝑖 𝑁𝑅 𝑗 𝑖=1
(5)

𝑆 𝑗𝑡

𝑖 are dummy variables for variances regimes detected by ICSS-heteroskedastic algorithm (see appendix A for details). * is an element-by-element matrix product, 𝑢 𝑡 the innovations vector and p and q the GARCH process order.

𝕋 𝑖 andℤ 𝑖 are coefficient matrix of volatility transmission and p and q are GARCH process orders. 𝕋 𝑖 is a coefficient matrix of elements �𝑡𝑐𝑣 𝑖 � for shock to volatility transmission between markets and ℤ 𝑖 a coefficient matrix of element �𝑡𝑣 𝑖 � for volatility transmission between markets. The volatility transmission regressors are defined as follows:

𝕋 𝑖 ′ 𝑢 𝑡-𝑖 𝑢 𝑡-𝑖 ′ 𝕋 𝑖 = � �∑ �𝑡𝑐𝑣 𝑗,1 𝑖 �𝑢 𝑗,𝑡-𝑖 𝑀 𝑗=2 � 2 ⋯ ⋯ ⋮ ⋱ ⋮ �∑ �𝑡𝑐𝑣 𝑗,1 𝑖 �𝑢 𝑗,𝑡-𝑖 𝑀 𝑗=2 ��∑ �𝑡𝑐𝑣 𝑗,𝑀 𝑖 �𝑢 𝑗,𝑡-𝑖 𝑀-1 𝑗=1 � ⋯ �∑ �𝑡𝑐𝑣 𝑗,𝑀 𝑖 �𝑢 𝑗,𝑡-𝑖 𝑀-1 𝑗=1 � 2 � and ℤ 𝑖 ′ 𝐻 𝑡-𝑖 ℤ 𝑖 = � ∑ �𝑡𝑣 𝑗,1 𝑖 � 2 ℎ 𝑗𝑗,𝑡-𝑖 𝑀 𝑗=2 ⋯ ⋯ ⋮ ⋱ ⋮ ⋮ ⋯ ∑ �𝑡𝑣 𝑗,𝑀 𝑖 � 2 ℎ 𝑗𝑗,𝑡-𝑖 𝑀-1 𝑗=1 � Volatility spillover named Volatility Causality is measured by the sum ∑ �𝑡𝑐𝑣 𝑖 � 2 + ∑ �𝑡𝑣 𝑖 � 2 𝑞 𝑖=1 𝑝 𝑖=1
. The last regressor ∑ (𝛿 𝑖 ) ′ 𝐷 𝑖 𝛿 𝑖 5 𝑖=1 permits « day of the week effects » and « holiday effects » in variance. 𝐷 𝑖 is a diagonal matrix whose element �𝐷 𝑗𝑗𝑡 � = 𝐷 𝑖𝑡 are :

𝐷 𝑖𝑡 = {𝐷 1𝑡 ; 𝐷 2𝑡 ; 𝐷 3𝑡 ; 𝐷 4𝑡 ; 𝐷 5𝑡 } (6)
Where 𝐷 1𝑡 = 1 (resp 𝐷 2𝑡 = 1, 𝐷 3𝑡 = 1, 𝐷 4𝑡 = 1 𝑎𝑛𝑑 𝐷 5𝑡 = 1) for Monday (resp, Tuesday, Wednesday, Friday and holiday day) and 𝐷 1𝑡 = 0 (resp 𝐷 2𝑡 = 0, 𝐷 3𝑡 = 0, 𝐷 4𝑡 = 0 𝑎𝑛𝑑 𝐷 5𝑡 = 0)otherwise. This effect may be present in daily data frequencies [START_REF] Solnik | Day-of-the-Week Effect on the Paris Bourse[END_REF][START_REF] Barone | The Italian Stock Market: Efficiency and Calendar Anomalies[END_REF][START_REF] Agrawal | Anomalies or Illusions? Evidence from Stock Markets in Eighteen Countries[END_REF]. 𝛿 𝑖 is a vector of coefficients to be estimated.

Transmission in means and variances: Estimation and tests

Estimation

The standard model and the BEKK-BSV model are estimated in two steps: first, we obtain the VAR residual, and secondly we estimate the conditional variance-covariance matrix parameters. This two-step estimation is possible because of the block-diagonal character of the variance-covariance matrix. Errors are conditionally normal and the likelihood function of all distribution is the sum of log-likelihood of each element. Let f t denote the joint conditional distribution and 𝐿 𝑁 𝑏 its log-likelihood function:

𝑙𝑜𝑔(𝐿 𝑁 𝑏 ) = � 𝑙𝑜𝑔�𝑓(𝑟 𝑡 𝛺 𝑡-1 ⁄ ; 𝜗)� 𝑁 𝐵 𝑡=1 𝑙𝑜𝑔�𝑓 𝑡 (𝑟 𝑡 𝛺 𝑡-1 ⁄ ; 𝜗)� = -(𝑀 2 ⁄ )𝑙𝑜𝑔2𝜋 -0.5𝑙𝑜𝑔 �𝐷𝑒𝑡�𝐻 𝑡,𝜗 �� -(𝑢 𝑡 ) 𝑡𝑟 �𝐻 𝑡,𝜗 � -1 𝑢 𝑡 2 ⁄ (7) 𝑙𝑜𝑔(𝐿 𝑁 𝑏 ) = -0.5 � �𝑀 𝑙𝑜𝑔(2𝜋) + 𝑙𝑜𝑔 �𝐷𝑒𝑡�𝐻 𝑡,𝜗 �� + (𝑢 𝑡 ) 𝑡𝑟 �𝐻 𝑡,𝜗 � -1 𝑢 𝑡 � 𝑁 𝑏 𝑡=1
Where 𝑁 𝑏 is the number of observations, M the number of markets, 𝜃 = 𝑣𝑒𝑐ℎ(𝐴, 𝐵, 𝐶 𝐵 , 𝐺, 𝐽, 𝛿, 𝜁 , 𝜉 and 𝜚) the vector parameters to be estimated, and 𝑢 𝑡 the normally distributed vector of innovation. Products 𝑢 𝑖𝑡 𝑢 𝑗𝑡 are second order correlated and the joint distribution of 𝑢 1𝑡 ; 𝑢 2𝑡 ; … ; 𝑢 𝑀𝑡 may not be normal. For this reason, θ is estimated by the Quasi-Maximum likelihood method (QML) of [START_REF] Bollerslev | Quasi Maximum Likelihood Estimation and Inference in Dynamic Models with Time-varying Covariances[END_REF]. Optimization is obtained with BHHH algorithm [START_REF] Berndt | Estimation and Inference in Nonlinear Structural Models[END_REF] which is very suitable for non-linear maximizations [START_REF] Engle | Multivariate Simultaneous Generalized ARCH[END_REF].

Tests

Mean spillover can be assessed using the Granger Causality test [START_REF] Lütkepohl | New Introduction to Multiple Time Series Analysis[END_REF]. These regimes clearly indicate a structural break in the variance during long-run bull and bear market volatility. These distinctions produce a more accurate market volatility model and permit a better comprehension of mean and volatility transmission and interdependences between stocks markets.

Empirical results

All

Mean and Prices Causality

Prices and asset returns spillover between stock exchanges markets. It is also shown that market prices are often transmitted unilaterally from US market to stock markets around the world. We produce measures of the mean transmission coefficient by the mean of linear VAR model.

Coefficients estimates show that recent information had more impact on returns than older information.

The bigger impact came from the American market. One day lagged American returns explain nearly 62% of European (Non EMU) prices, 60% of EMU prices and 48% of Asian markets prices. Markets returns are also explained by own lagged returns. Prices are also transmitted from European (Non EMU) markets to Asian and EMU markets (table III, panel A).

Looking for Granger Causality between markets confirms the largest impact of American markets on all other markets. Granger Causality test of non-transmission is accepted from Asian market to all other markets (table III, panel B).

Long-Run Mean Spillover (LRMS hereafter) is measured by

Φ = �I 4 -∑ Φ i 5 i=1 � -1
. One can see that American market LRMS estimates are more important than others LRMS. At long-run more than 50% of European (Non EMU) and Asian prices are exogenous and determined by American prices. It is about 49% in the case of EMUs markets (table III, panel C).

These findings are confirmed with Cumulative Impulse-Response Functions (CIRF hereafter). Figure 2 plots CIRF for linear VAR. CIRF of the American aggregate market are much important than others. There is no accumulate responses of American market to others market. American prices seem to be exogenous (figure2, panel A); whether Asian and European (Non EMU) prices are highly affected by American prices (figure 2, panel B and C).

EMU aggregate markets are highly responsive to American market and quite responsive to European market. The long-run effect of American market is more important than the endogenous prices of EMU theme self (figure 2, panel D).

Concerning the "Days of the Week" effects and according to our results, there are no such effects.

Concerning the causality during crisis, one can see a higher prices transmission during subprime crisis compared to transmission during the entire period (table IV, panels A and B). American aggregate prices are still affecting all others market (table IV, panel C). According to Wald test of coefficient equality, mean transmission coefficients are significantly different during the crisis than those of the entire period (table IV, panel D). First, all asymmetric coefficients are significant and confirm the asymmetric behavior of market's volatility according to positive and negative shocks. The American market is the most asymmetric one

Conditional Variances-Covariance's Estimates

(table V panel A and table VI panel A).

Second, in the case of the standard BEKK model, volatility persistence is close to unity (>.94) particularly for American markets. Results show lesser persistence estimates with BEKK-BSV. Half-life of volatility is tow times shorter than that estimated by the standard model and varies from 1 day for the Asian market to 8 days for the American market. Volatility is much more persistent in the American markets than in the Asian markets. Hence, the standard model without multiple regimes variances tends to overestimate the persistence (Lamoureux and Lastrapes, 1994). This overestimation leads to a bias in market spillover estimates and market interdependence estimates (table V There are three main results:

First of all, transmission between markets during crisis period is higher than transmission outside the crisis. The transmission from the American market to European markets almost doubled in intensity during the crisis: from 5% to 15% for European market, from 13% to 29% for Asian market and from 4% to 13% for EMU market (table VI pane, B).

Second, there are also new transmissions during subprime crisis: for example, during crisis period, 26% of European (Non EMU) market is transmitted to the Asian markets (table VI pane, B). These changes during the crisis undoubtedly indicate shift contagion [START_REF] Forbes | No contagion, only interdependence: measuring stock market comovements[END_REF].

Concerning the stability of transmission during crisis, Wald test of no-strict contagion during subprime crisis indicate that volatility transmission from American market to European markets was significantly different from the one during calm period. This results support the shift contagion theory. There is also contagion from European non-EMU market to EMU market (table VI panel C).

4.

In this empirical paper, we demonstrate that linear VAR is misspecified because of the significance of second order autocorrelations. This is a serious bias in mean spillover estimates if we consider as here the links between the American, the European and the Asian markets.

Concluding remark

This paper highlights the importance of structural breaks in variance. The standard BEKK model tends to overestimate the persistence of volatility shocks, while also overestimating the endogenous part of volatility in market uncertainty. By introducing multiple regimes of variance, we reduced the persistence and improved the model in terms of maximum log-likelihood. The Granger Causality test revealed transmission from the American to European markets in means and variances. Moreover, it appears that it takes a week for information to filter through to all the markets. We have shown that the standard BEKK model underestimates the transmission variance, including the impact of the American market on European markets. The model with structural breaks in variance confirms unilateral transmission from the American market to European and Asian markets.

This study of the transmission of uncertainty Markets has confirmed the key role of the market in North America region. It appears that almost half the prices of European markets and Asian markets are determined by the U.S. market. This price transmission is valid in quiet times and in times of crisis.

Having regard to the transmission of volatility, it also appears that the uncertainty of the U.S. market is transmitted unilaterally to other regions. This transmission is very important in the case of the Asian markets. The market impact is also important in European markets outside the euro area and in the markets of the euro area. In times of crisis , such as the subprime crisis between July 2007 and July 2009 , the transmission of American uncertainty to other regions has increased significantly. This episode of crisis that shook the world of finance seems to be contagious. Our results showed the appearance of new transmission channels during the crisis and support the contingent crises, for which the transmission mechanism during the crisis ( or just after ) is fundamentally different from that which prevailed before the crisis theory. The latter causes a structural change such that the shocks are propagated through a channel that does not exist in the periods of financial stability. The regional analysis of aggregated markets showed the vulnerability of markets in the euro area face shocks to U.S. and European shocks outside the euro zone.

Figure 2 : Cumulative Impulse-Response-Function CIRF.

CIRF are from linear VAR(5) for 10 periods. Legend: y1 (America), y2 (Europe exclude EMU), y3 (Asia) and y4 (EMU). = 33). ***, **, * and # Significant coefficients at 1%, 5%, 10% and 20%. . A, B and C Subprime crisis period. D Debt crisis.

a Under null hypothesis of no contagion, Wald statistic follow 𝜒(2).

  Cholesky One S.D. Innovations ± 2 S.E.

Figure 3 :

 3 Figure 3: American, European, Asian and EMU market's Volatility estimates

  tests of variance effects (volatility transmission, variance regime are performed by the mean of likelihood Ratio (LR) test and/or Wald test. Mean transmission from market « i » to market « j » is measured by Φ 𝑖𝑗 The data cover America, aggregate Europe, aggregate EMU and aggregate Asia markets. MSCI standard capitalization indices are used since they are better adapted than usual market indices. Indeed they include the mid-cap and the large-cap companies' capitalization. Daily observations from June, 2005 to October 4, 2013 are used.

	3.1.	Data descriptive

1 … Φ 𝑖𝑗 𝑛1 VAR coefficients. No transmission in mean is a Granger-Causality test called NMT. Variance transmission is linked to �𝑡𝑐𝑣 𝑖𝑗 𝑘 � and �𝑡𝑣 𝑖𝑗 𝑘 � coefficients. There are no variance transmission (NVT thereafter) from market « i » to market « j » when all �𝑡𝑐𝑣 𝑖𝑗 𝑘 � and �𝑡𝑣 𝑖𝑗 𝑘 � coefficients are zero. NVT is a Wald-test of No-Causality in variance. Finally, we test the equality of volatility transmission during calm period and crisis period. This is Wald test of No strict contagion (see Appendix. C). Next we discuss the main empirical results.

Table I :

 I Descriptive statistics of markets returns. statistics show the usual characteristics of financial data: asymmetry, excess kurtosis and nonnormality (TableI). MSCI indices returns are weak. The EMU markets are more volatile than others according to standard deviation. The asymmetry is most pronounced in Asian and American markets than the EMU and European markets. Excess kurtosis shows that extreme values are more frequent than predicted by normality. The Jarque-Bera statistic confirms the non-normality of data.The Variance regimes detected by the ICSS-H algorithm are given in table II and represented in figure1. The number of regimes is 11 for Europe, 10 for America and 8 for the EMU and Asian markets. Start and end dates of regimes are not equal although there are some coincidences during the multiple stage of subprime crisis period between July 2007 till December 2009:

	3.2.	Conditional variance model
		Mean Max	Min Std. Dev. Skewness Kurtosis	JB	N. Obs
	AMERICA EUROPE Non-EMU 0.017 10.885 -9.907 0.018 10.626 -9.708 EMU 0.012 10.869 -10.443	1.289 1.420 1.653	-0.316 -0.110 -0.086	12.371 10876.920 10.650 7221.861 8.164 3291.764	*** *** ***	2959 2959 2959
	ASIA	0.017 9.691 -8.848	1.258	-0.330	8.129 3297.742	2959

*** *** Significant at 1%. JB Jarque-Bera statistic.

Descriptive

Table II :

 II Number and date of variance regimes. Variance regimes dates obtained by ICSS-H algorithm. Start: the beginning date of the regime. End: end date of the scheme. Variance: the measure of non-conditional variance during the regime. A, B and C Subprime crisis period. D Debt crisis Return to a quite calm period from July 2009 to June 2011. • 07/2011-12/2011: European debt crisis. • Return to a quite calm period from July 2011 to nowadays.

			America markets					Europe Non Emu
	Regime	Start	End	Variance		Regime	Start	End	Variance
	1 03/06/2005 20/07/2007	0.649		1 03/06/2005 20/07/2007	0.832
	2 23/07/2007 11/09/2008	1.263	A	2	23/07/2007 10/09/2008	1.487	A
	3	12/09/2008 28/11/2008	4.679	B	3	11/09/2008 05/12/2008	4.648	B
	4	01/12/2008 15/05/2009	2.301	C	4	08/12/2008 14/07/2009	2.155	C
	5	18/05/2009 06/11/2009	1.248		5 15/07/2009 16/12/2009	1.443
	6 09/11/2009 23/04/2010	0.823		6 17/12/2009 23/04/2010	1.011
	7 26/04/2010 31/08/2010	1.496		7 26/04/2010 10/08/2010	1.844
	8 01/09/2010 29/07/2011	0.795		8 11/08/2010 01/08/2011	1.080
	9 01/08/2011 19/12/2011	2.039	D	9	02/08/2011 14/12/2011	2.227	D
	10	20/12/2011 01/10/2013	0.755		10 15/12/2011 02/08/2012	1.237
						11 03/08/2012 01/10/2013	0.836
			Emu					Asia
	Regime	Start	End	Variance		Regime	Start	End	Variance
	1 03/06/2005 20/07/2007	0.911		1 03/06/2005 10/05/2006	0.909
	2 23/07/2007 02/09/2008	1.459	A	2	11/05/2006 07/08/2006	1.623
	3 03/09/2008 05/12/2008	4.528	B	3	08/08/2006 25/07/2007	0.829
	4 08/12/2008 14/07/2009	2.384	C	4	26/07/2007 04/09/2008	1.507	A
	5	15/07/2009 06/07/2011	1.627		5 05/09/2008 18/05/2009	2.606	B, C
	6	07/07/2011 29/11/2011	2.953	D	6	19/05/2009 28/07/2011	1.010
	7 30/11/2011 02/08/2012	1.793		7 29/07/2011 20/12/2011	1.608	D
	8	03/08/2012 01/10/2013	1.178		8 21/12/2011 01/10/2013	0.878

• 07/2007 -09/2008: the first phase of the sub-prime crisis. • 09/2008 -12/2008: the most turbulent regime and the most volatile during the second phase of the economic and financial crisis. • 12/2008 -07/2009: a third phase of the crisis which is less violent than the two previous phases. •

  panel A and table VI panel A). Concerning days of the week and holiday's effects we did not found any significant results9 The LR test confirms that the model with structural breaks is better than the standard model (table VI panel B). Figure2shows the volatilities estimated by BEKK-BSV models for the four aggregate markets between 2005 and 2013. Figures indicate clearly high volatility regimes during Sub-prime crisis period between July 2007 and April 2009 and during economic crisis period from 2010 to 2011

	Third, the omegas estimates in the conditional variance equation are higher during the crisis periods:
	regimes 2, 3 and 4 for the American market, European (Non EMU) markets and EMU markets. Those
	regimes coincide with the multiple steps of the sub-prime crisis between July 2007 and July 2009 (table
	VI panel A) 8
	3.5.	Volatility spillover and Variance Causality
			.
		10	.
	3.6.	Volatility spillover during crisis period
	Are the transmission phenomena different during crises? We identify crisis period from 23 July 2007 to
	14 July 2009. This sub-period of turbulence is determined by successive more volatile regimes obtained
	by ICSS-H algorithm 11 .

. A very high volatility regime in the American market from Volatility spillover is given for both models (tables V and VI). The comparison between the standard model and the BEKK-BSV model shows important differences: the standard model underestimates the intensity of European market exposure to American market turbulence. Results confirm the unilateral spillover from American markets to al others regions (table VI panel B).

Overall, the share of the volatility of the American market in the European market volatility is about 5%, it is about 13% to Asian market and 4.5% to EMU market. Volatility spillover estimates from American market to European market obtained with standard model is two times less than the one obtained by structural break model. Feedback from European and Asian markets are low (less than 1%) and non significant. There is no sensitivity of the European and EMU markets to the turbulence of the Asian market and vice-versa. Within Europe, there is no volatility spillover between EMU and non-EMU market.

8 High volatility regimes are indicated by bolt script.

9

These results are not presented in the paper but can be provide upon request. 10 Unconditional variance measurements confirm the presence of high and low variance regimes (table 4). The measurement of unconditional variance is E(ℎ 𝑡 ) = 𝜔 1-𝜑 . Estimates are available upon request. 11 See Table 1 for regimes detected by ICC-H algorithm.

Table III :

 III Price return transmission, Mean Causality and Long-Run Mean Spillover.

		Panel A: Mean Causality (Entire period 2006-2013)			Mean Causality (Entire period 2006-2013)
		America	Europe	Asia		EMU		America	Europe	Asia	EMU
	America (-1)	-0.103 ***	0.626 ***	0.484 ***	0.609 *** Asia (-1)	0.003	0.017	-0.332 ***	0.025
	(-2)	-0.011	0.252 ***	0.190 ***	0.231 ***	(-2)	-0.004	-0.006	-0.146 *** -0.007
	(-3)	0.056	0.150 ***	0.131 ***	0.138 ***	(-3)	-0.002	-0.023	-0.079 *** -0.023
	(-4)	-0.013	0.061	0.028		0.047	(-4)	-0.007	-0.006	-0.056 **	-0.009
	(-5)	-0.016	0.095 ***	-0.033		0.048	(-5)	0.015	0.003	-0.020	-0.001
	Europe (-1)	-0.104 *	-0.371 ***	0.098 **		-0.395 *** EMU (-1)	0.114 **	-0.053	0.080 **	-0.022
	(-2)	-0.141 **	-0.193 ***	-0.006		-0.158 **	(-2)	0.068	-0.019	0.055	-0.036
	(-3)	-0.044	-0.116 *	0.006		-0.039	(-3)	0.012	-0.017	0.021	-0.075
	(-4)	0.009	0.016	0.082 *		0.061	(-4)	0.025	-0.020	-0.029	-0.046
	(-5)	-0.027	-0.122 **	-0.018		-0.104	(-5)	-0.027	0.001	0.096 **	0.011
		Panel B: Granger Causality tests (entire period)	b			Panel C: Long-Run Mean Spillover LRMS	a
		America	Europe	Asia		EMU		America	Europe	Asia	EMU
	America	16.737 ***	408.012 ***	461.665 ***	268.993 *** America	0.858	0.535	0.541	0.491
	Europe	7.842	47.046 ***	7.668		33.951 *** Europe	-0.204	0.452	-0.114	-0.432
	Asia	0.522	0.917	207.021 ***	0.911	Asia	0.003	-0.003	0.613	-0.003
	EMU	7.811	1.434	13.676 **		2.594	EMU	0.160	0.046	0.216	0.976
	(..) Delay. ***, ** and * Significant coefficients at 1%, 5% and 10%. hypothesis of non Granger Causality test statistic follow 𝜒 (5) 2 . C Under null hypothesis of Equal transmission test statistic follow 𝜒 (5) a Long-run mean spillover is measured as Φ = (I m -Φ 1 ⋯ -Φ 5 ) -1 . 2 .	b Under null

Table IV :

 IV Price return transmission, Mean Causality and Long-Run Mean Spillover during Crisis periods.

			Panel A:USA Subprime Crisis (2007-2009)		Panel B: EMU Debt Crisis (2010-2011)
			America	Europe	Asia	EMU	America	Europe	Asia	EMU
	America	(-1)	-0.136 ***	0.677 ***	0.451 ***	0.631 ***	-0.022 --	0.547 ***	0.481 ***	0.565 ***
		(-2)	-0.013	0.286 ***	0.165 ***	0.218 ***	0.150 ***	0.339 ***	0.224 ***	0.341 ***
		(-3)	0.150 **	0.218 ***	0.178 ***	0.172 ***	-0.095 **	0.020 --	0.107 ***	0.030 --
		(-4)	-0.054	0.039	0.039	-0.004	0.089 *	0.040 --	0.077 **	0.056 --
		(-5)	0.014	0.144 ***	-0.006	0.080	-0.143 ***	-0.123 ***	-0.089 ***	-0.229 ***
	Europe	(-1)	-0.162	-0.320 ***	0.117 *	-0.223 **	0.046 --	-0.355 ***	0.019 --	-0.545 ***
		(-2)	-0.205 **	-0.096	-0.021	-0.048 --	0.017 --	-0.177 **	0.159 ***	-0.154 --
		(-3)	-0.177 *	-0.203 **	-0.079	-0.103 --	-0.086 --	0.011 --	0.107 **	0.050 --
		(-4)	-0.114	0.003	0.096	0.052 --	0.059 --	-0.017 --	-0.004 --	0.016 --
		(-5)	-0.035	-0.213 **	-0.167 **	-0.230 **	0.139 **	0.282 ***	0.154 ***	0.270 **
	Asia	(-1)	0.030	0.017	-0.397 ***	0.046 --	-0.014 --	0.001 --	-0.342 ***	0.020 --
		(-2)	-0.065	-0.039	-0.186 ***	-0.045 --	0.037 --	0.012 --	-0.246 ***	0.019 --
		(-3)	0.000	0.010	-0.130 ***	0.019 --	-0.155 ***	-0.232 ***	-0.142 ***	-0.285 ***
		(-4)	-0.013	-0.009	-0.120 ***	-0.005 --	-0.041 --	-0.026 --	-0.045 *	-0.018 --
		(-5)	0.067 *	0.049	-0.046 *	0.049 --	-0.052 *	-0.151 ***	-0.089 ***	-0.201 ***
	EMU	(-1)	0.170 *	-0.154 *	0.138 **	-0.214 **	-0.028 --	-0.008 --	0.074 **	0.134 *
		(-2)	0.051	-0.175 **	0.102 *	-0.155 *	-0.076 --	-0.057 --	-0.063 *	-0.114 --
		(-3)	0.153 *	0.054	0.110 *	0.013 --	0.072 --	-0.032 --	-0.004 --	-0.074 --
		(-4)	0.194 **	0.030	-0.003	0.032 --	-0.032 --	0.051 --	0.021 --	0.026 --
		(-5)	-0.088	0.019	0.288 ***	0.073 --	-0.014 --	-0.118 **	-0.025 --	-0.053 --

Panel C: Granger Causality tests (crisis period) b

		America	Europe	Asia	EMU	America	Europe	Asia	EMU
	America	27.514 ***	17.246 ***		12.116 **	10.059 *	24.282 ***	21.471 ***	4.836	19.551 ***
	Europe	4.872	6.405		7.898	7.345	4.601	9.962 *	8.856	7.808
	Asia	12.518 **	8.269		13.325 **	8.660	8.156	15.033 **	6.788	18.573 ***
	EMU	7.262	3.971		10.336 *	5.509	4.627	2.448	8.318	3.711
		Panel D: Equality of transmission	c
		America	Europe	Asia	EMU	America	Europe	Asia	EMU
	America	18.762 ***	19.487 ***		55.659 ***	19.958 ***	27.329 ***	79.690 *** 46.512 *** 54.460 ***
	Europe	2.094	8.308		4.406	12.441 **	7.108	16.199 *** 11.258 **	8.443
	Asia	18.609 ***	6.036		9.900 *	5.459	16.690 ***	22.200 *** 22.514 *** 17.411 ***
	EMU	3.972	2.480		3.935	4.526	7.315	2.537	15.977 ***	3.982
	(..) Delay. ***, ** and * Significant coefficients at 1%, 5% and 10%. a Long-run mean spillover is measured as Φ = (I m -Φ 1 ⋯ -Φ 5 ) -1 .

b Under null hypothesis of non Granger Causality test statistic follow 𝜒 (5) 2 . C Under null hypothesis of Equal transmission test statistic follow 𝜒 (5) 2 .

Table V :

 V Conditional variances coefficients estimates, Variance spillover and Variance Causality estimates for BEKK model

	Panel A: Coefficients estimates
		America	Europe	Asia	EMU
	𝜔 0 2	0.140 ***	0.112 ***	0.055 **	0.123 ***
	𝛼 2	0.000 ***	0.009 ***	0.090 ***	0.014 ***
	𝛾 2	0.111 ***	0.067 ***	0.000 ***	0.073 ***
	𝛽 2	0.903 ***	0.906 ***	0.726 ***	0.914 ***
	Panel B : Persistence and Half-life
	Persistence	0.959	0.948	0.816	0.964
	Half-life	16.388	13.080	3.414	19.010
	Panel C: Variance spillover
		America	Europe	Asia	EMU
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐴𝑀𝐸𝑅𝐼𝐶𝐴,𝑖	--	3,72% ***	8,51% ***	2,70% ***
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐸𝑈𝑅𝑂𝑃𝐸,𝑖	3,82% **	--	1,78%	0,05%
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐴𝑆𝐼𝐴,𝑖	0,20%	0,14%	--	0,09%
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐸𝑀𝑈,𝑖	0,63%	0,04%	0,16%	--
	Log-likelihood	-5943.683
	AIC	5.604

a (..) Variance Causality Wald test significance. ***, ** and * Significant coefficients at 1%, 5% and 10%. 20

Table VI :

 VI Conditional variances coefficients estimates, Variance spillover and Variance Causality estimates.

	Panel A: Coefficients estimates
		America	Regime period		Europe				Asia	EMU
	2 𝜔 0 Regime 2	0.3963 *** 1.2886 ***	03/06/2005 20/07/2007 23/07/2007 11/09/2008	A	0.5849 *** 03/06/2005 20/07/2007 0.8339 *** 23/07/2007 10/09/2008	0.5017 *** 0.8337	03/06/2005 11/05/2006	10/05/2006 07/08/2006	0.7788 *** 03/06/2005 20/07/2007 0.7800 23/07/2007 02/09/2008
	Regime 3	7.9544 ***	12/09/2008 28/11/2008	B	1.5817	*** 11/09/2008 05/12/2008	0.2535 ***	08/08/2006	25/07/2007	0.7237	03/09/2008 05/12/2008
	Regime 4	1.0863 ***	01/12/2008 15/05/2009	C	0.3690	*** 08/12/2008 14/07/2009	0.3994	26/07/2007	04/09/2008	1.0130 *** 08/12/2008 14/07/2009
	Regime 5	0.7107 **	18/05/2009 06/11/2009		0.8073 *** 15/07/2009 16/12/2009	0.6497	05/09/2008	18/05/2009	0.9507 *** 15/07/2009 06/07/2011
	Regime 6	0.5022		09/11/2009 23/04/2010		0.5390	17/12/2009 23/04/2010	0.1579 ***	19/05/2009	28/07/2011	1.3647 *** 07/07/2011 29/11/2011
	Regime 7	0.9492 **	26/04/2010 31/08/2010		0.7268	26/04/2010 10/08/2010	0.1501	29/07/2011	20/12/2011	1.2551 *** 30/11/2011 02/08/2012
	Regime 8	0.4556		01/09/2010 29/07/2011		0.5585	11/08/2010 01/08/2011	0.2662 ***	21/12/2011	01/10/2013	0.9015 *** 03/08/2012 01/10/2013
	Regime 9	1.3017		01/08/2011 19/12/2011	D	0.5273	02/08/2011 14/12/2011	--	--
	Regime 10 0.3991		20/12/2011 01/10/2013		0.5696	15/12/2011 02/08/2012	--	--
	Regime 11	--					0.4227 *** 03/08/2012 01/10/2013	--	--
	𝛼 2	0.004 ***				0.015 ***			0.005 ***	0.027 ***
	𝛾 2	0.106 ***				0.077 ***			0.115 ***	0.072 ***
	𝛽 2	0.865 ***				0.849 ***			0.421 ***	0.864 ***
	Persistence	0.922					0.903				0.483	0.926
	Half-life	8.482					6.765				0.953	9.052
	Panel B: Variance spillover
						Entire Period					Crisis Period
		America	Europe			Asia		EMU		America	Europe	Asia	EMU
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐴𝑀𝐸𝑅𝐼𝐶𝐴,𝑖 --		5,08%	***	13,87% **	4,47%	***	--	15.48% ***	29.34% ***	13.92% ***
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐸𝑈𝑅𝑂𝑃𝐸,𝑖	0,27%	--			7,82%		0,02%		1.26% *	--	26.26% ***	0.20% *
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐴𝑆𝐼𝐴,𝑖	1,99%	0,77%		--			1,21%		4.20%	0.96%	--	1.66%
	(𝑡𝑐𝑣² + 𝑡𝑣²) 𝐸𝑀𝑈,𝑖	0,93%	0,17%			0,14%	--			1.74%	0.54%	0.42%	--

Panel C : Wald test of non-contagion (NSC) e America Europe Asia EMU

  Under null hypothesis of no break in variance, LR statistic follow 𝜒 (33)

	America	--	3.983 **	7.862 ***	6.145 ***
	Europe	3.446 **	--	4.256 **	4.226 **
	Asia	1.236	2.555	--	0.026
	EMU	1.786	2.265	1.993	--
	Log-likelihood	-5843.18
	AIC	5.540
	LR test	201.006
					2	(∑ 𝑁𝑅 𝑖 4 𝑖=1

c *** b (..) Variance Causality Wald test significance. c

Edwards (1998) used a GARCH model to investigate the national and regional impact of interest rate spread during the Tequila crisis.Park and Song (2000) applied a similar GARCH model to test volatility spillover across Asian markets during the Asian flu epidemic.

Φ(L) = I m -Φ 1 L ⋯ -Φ n1 L n1, where n1 is VAR process order defined by sequential LR test (see details in appendix A.).

The DCC models of[START_REF] Engle | T heoretical and Empirical Properties of Dynamic Conditional Correlation Multivariate GARCH[END_REF] and[START_REF] Tse | A Multivariate GARCH Model with Time-varying Correlations[END_REF] have the advantage of a two-step estimation. However, these templates provide a linear structure to the correlation dynamics and impose a similar dynamic conditional correlation. In addition, these models do not permit variance spillover.

Appendix A

Figure 1 : Variance regimes (May, 2005-October, 2013).

Variance regimes are obtained by ICSS-H algorithm applied to MSCI yields (mid-cap and large-cap) for America, Europe, Asia and EMU. Shading areas indicate crisis periods. We use an ICSS algorithm based on the CUSUM test to detect the structural change in variance. Following Inclan and Tiao (1994), the variance of a given series shows a structural change due to an exogenous shock. These changes mean a permanent decline in the tendency which continues until the appearance of a new significant shock. This analysis supposes a stationary variance between two points of structural change. Let M series of independent and normally distributed observations: 𝑥 𝑖,𝑡 (𝑖 = 1, 𝑀 ������ ). The non-conditional variance of each one of them is 𝜎 𝑖,𝑡 2 and 𝑁𝑅 𝑖 the number of break point in the variance. On the whole sample of N observations, we have:

Where the 𝑐 𝑖,𝑗 (j=1…𝑁𝑅 𝑖 ) are the dates of break in variance. To estimate the number of changes of variance tendency, the cumulative sum of the square residuals

is calculated. Inclan and Tiao (1994) define the statistics For each series 𝑥 𝑖,𝑡 , there are 𝑁𝑅 𝑖 break points in variance and 𝑁𝑅 𝑖 + 1 régimes of variances.

12 𝛿 ̂4 is obtained from the non parametric estimator :

, where w(l,m) is a Bartlett window. The 𝛿 ̂4 estimates depend on the choice of m parameter with the Newey-West method. It's usually equal to 2.

No Transmission tests and

𝑁𝑉𝑇

Under 𝐻 0,𝑖→𝑗 𝑁𝑉𝑇 , Wald statistic follow a 𝜒 2 (𝑝 + 𝑞).

No Strict Contagion test

There is strict contagion when transmission during crisis period is significantly different then


Nina Miolane
email: nina.miolane@inria.fr

Xavier Pennec

Biased estimators on Quotient spaces

Usual statistics are defined, studied and implemented on Euclidean spaces. But what about statistics on other mathematical spaces, like manifolds with additional properties: Lie groups, Quotient spaces, Stratified spaces etc? How can we describe the interaction between statistics and geometry? The structure of Quotient space in particular is widely used to model data, for example every time one deals with shape data. These can be shapes of constellations in Astronomy, shapes of human organs in Computational Anatomy, shapes of skulls in Palaeontology, etc. Given this broad field of applications, statistics on shapes -and more generally on observations belonging to quotient spaces-have been studied since the 1980's. However, most theories model the variability in the shapes but do not take into account the noise on the observations themselves. In this paper, we show that statistics on quotient spaces are biased and even inconsistent when one takes into account the noise. In particular, some algorithms of template estimation in Computational Anatomy are biased and inconsistent. Our development thus gives a first theoretical geometric explanation of an experimentally observed phenomenon. A biased estimator is not necessarily a problem. In statistics, it is a general rule of thumb that a bias can be neglected for example when it represents less than 0.25 of the variance of the estimator. We can also think about neglecting the bias when it is low compared to the signal we estimate. In view of the applications, we thus characterize geometrically the situations when the bias can be neglected with respect to the situations when it must be corrected.

Introduction

In Quantum Field theory, one can roughly consider a relativistic quantum field of quarks as a function from space-time to the space of colors: φ : R 4 → C 3 . The function φ associates to each point of space-time an elementary particle: the quark living in the space of colors. The space of quarks fields carries a right action of the group of symmetry of space-time and a left action of the group of quark symmetry. These group actions represent what does not change the laws of Physics.

In Computational Anatomy when we study anatomical shapes, we encounter the same mathematical structures [START_REF] Grenander | Computational anatomy: An emerging discipline[END_REF]. In others words, we have often two groups acting on the same anatomical object φ respectively on the right and on the left. These group actions represent what does not change the shape of the object.

A first example is a set of K anatomical landmarks [START_REF] Kendall | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF][START_REF] Le | The riemannian structure of euclidean shape spaces: A novel environment for statistics[END_REF], which is a function from the space of labels to the real 3D space: φ : [1, ..., K] → R 3 . In this case, the manifold M of landmarks sets carries a right action of the group of permutations of the labels and a left action of the group of rotations and translations of the 3D space SE [START_REF] Bauer | Overview of the geometries of shape spaces and diffeomorphism groups[END_REF]. Relabeling the landmarks, or rotating and translating the whole set of landmarks, does not change the anatomical shape described. A second example is a 1D-signal [START_REF] Kurtek | Signal estimation under random timewarpings and nonlinear signal alignment[END_REF] which is a function: φ : [0, 1] → R. The manifold M of 1D-signals carries a right action of the group of diffeomorphisms of [0, 1]: this action represents the reparameterizations of the signal, which do not change its shape. A third example is a parameterized surface, or more generally immersed submanifold, in R d [START_REF] Bauer | Overview of the geometries of shape spaces and diffeomorphism groups[END_REF] which is a function: φ : N → R d . The manifold M of parameterized surfaces (resp. of immersed submanifolds) carries a right action of the group of diffeomorphisms of N and a left action of the group of rotations and translations of the R d space SE(d) (poses). Again, reparameterizing, rotating or translating the surface or the submanifold does not change its shape.

The anatomical data that we observe belong to the manifold M. We are interested in the statistical analysis of their shapes. The corresponding shape data now belong to a quotient space Q, which is the manifold M quotiented by the group actions that leave the shape invariant. Thus we perform statistics on Q rather than on M, even if the data originally belong to M. But Q not a differentiable manifold in the general case. To be more precise, Q is a manifold with singularities or a stratified space.

Statistics on quotient spaces or shape spaces Q have been studied since the 1980's. Regarding shapes of landmarks, the theory has been first introduced by Kendall [START_REF] Kendall | A survey of the statistical theory of shape[END_REF]. Regarding shapes of curves one can refer to [START_REF] Kurtek | Signal estimation under random timewarpings and nonlinear signal alignment[END_REF][START_REF] Joshi | Riemannian structures on shape spaces: A framework for statistical inferences[END_REF] and, in the more general setting of manifold quotiented by isometric group actions, to [START_REF] Huckemann | Intrinsic shape analysis: Geodesic principal component analysis for riemannian manifolds modulo lie group actions[END_REF]. In all cases when one performs statistics on shape spaces, the estimation of the mean shape is central. A way of doing it is via parametric statistical estimation. By doing so, studies model the variability in the shapes of the objects (meaning that they assume some variability in Q) and the variability in reparameterization or in objects pose (meaning that they assume some variability in the "orbit" of the group actions).

However, the effect of noise on the objects themselves (which is a noise in the ambient manifold M) has not been thoroughly investigated yet. The noise on the objects in M exists as it comes from imperfect measure instruments or human lack of precision, for example while positioning the landmarks on a medical image. Usual estimators were proven to be consistent in the theory without noise. But experiments have shown that they have bias and even inconsistency in the presence of noise [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]. For example, the standard estimator of the mean shape of the data, the Fréchet mean computed in the quotient space with the max-max algorithm, is biased and even inconsistent. This bias is experimentally observed and shows to be dependent on the noise level. But it is not theoretically understood so far. There is a need of a extended statistical theory on quotient spaces that takes the noise into account and quantifies the induced bias. This paper is a first step in this direction. We suggest a geometric interpretation of the usual estimator of the mean shape. It enables to show that noise on the observations in M induces bias on the estimator in Q. We work in the case where the observations belong to a finite dimensional flat manifold M, which is quotiented by a proper isometric action of a finite dimensional Lie group G. We describe how the bias on the statistical estimator is controlled by the geometry of the quotient space Q, in particular by its singularities at the scale of the noise. Even if we work in finite dimension, we provide the intuition of the behavior for infinite dimension.

In the first section, we present notions of differential geometry on quotient spaces. In the second section, we present the statistical framework and the geometric interpretation of the estimator usually computed. In the third section, we demonstrate the geometric origin of the bias on this estimator. In the fourth section, we show explicit computations of the bias for the case of R m quotiented by SO(m), which is the common example of finite dimensional flat manifold quotiented by the isometric action of a finite dimensional Lie group.

1 Differential geometry of Quotient spaces Q

We consider M a finite dimensional flat Riemannian manifold and G a finite dimensional Lie group. For more details on this section, we refer to [1, ?,?].

Basis on Quotient Spaces

A Lie group action of G on M is a differentiable map ρ : G×M (g, X) → g•X ∈ M, such that e•X = X and g •(g•X) = (g g)•X. If G × M (g, X) → (g • X, X) ∈ M × M,
is a proper mapping, i.e. if the inverse image of every compact is compact, the action is proper. If dρ g : T X M → T g•X M with ρ g = ρ(g, .) leaves the metric of M invariant, the action is isometric. We only consider proper isometric actions.

The orbit of

X ∈ M is defined as O X = {g • X | g ∈ G}.
The orbit of X represents all the points that can be reached by moving X through the action of G. The isotropy group of X, also called the stabilizer of X, is defined as

G X = {g ∈ G | g • X = X}.
The stabilizer of X represents the elements of G that fail at moving X. The orbit O X is a submanifold of M and the stabilizer G X is a Lie subgroup of G. They are related by the orbit-stabilizer theorem as follows: O X ∼ G/G X . The orbits form a partition of M and the set of orbits is called the quotient space Q = M/G. Because the action is proper, Q is Hausdorff and inherits a Riemannian structure from M. This is precisely why we consider a proper action.

The isotropy groups G X and G X of elements X and X in the same orbit are conjugate groups in G. This property enables to define the orbit type of an orbit as the conjugacy class (H) of the isotropy groups of its elements. From the orbit-stabilizer theorem, we observe that a "large" orbit leads to a "small" isotropy group, so a "small" orbit type. The smallest orbit type is called the principal orbit type. The corresponding largest orbits are called principal orbits. In contrast, non principal orbits are called singular orbits. The corresponding points in Q are called singularities. The action of G on M is said to be free if there is no singular orbits. In this case, Q has no singularities and is a differentiable manifold. Otherwise Q is a manifold with singularities.

We investigate how the singular orbits are distributed in M, or equivalently how the singularities are distributed in Q. The partition of M into orbits is a singular Riemannian foliation called the canonical foliation of M. We can gather the orbits of same orbit type (H) and define M (H) = {X ∈ M| G X conjugate to H}. The decomposition of M into the connected components of the M (H) is a stratification called the orbit type stratification of M. In this stratification, the principal type component is open and dense.

The decomposition of Q into the corresponding components Q (H) = M (H) /G also forms a stratification, whose principal type component is also dense. Thus, singularities are "sparsely distributed" in Q in the sense that they are of null Lebesgue measure. These foliation and the stratifications are illustrated in Figure 1 in the case of M = R 2 quotiented by the action of SO(2) so that Q = R 2 /SO(2) = R + . We use this 2D example throughout the paper as it is convenient for the illustrations. But of course, the theory applies to the general case.

Principal strata Singular Strata

Manifold M = R 2 Quotient space Q = R+
Orbits as immersed submanifolds An orbit O is an immersed submanifold in M, i.e. a differentiable manifold together with an immersion I : O → M. We identify O with I(O) and we denote T X O, resp. N X O, the tangent space, resp. the normal space with respect to the metric of M, of O at X. The first fundamental form of O is defined to be the induced metric on O, i.e. the pullback of the ambient (flat) metric of M by I. The second fundamental form of O is the quadratic form h defined by T

X O × T X O (u, v) → h(v, w) = (∇ u w) ⊥ ∈ N X O,
where ∇ is the Levi-Civita connection on M and (∇ v w) ⊥ the orthogonal projection on N X O.

We can write O locally around X as the graph G of a smooth function from T X O to N X O. Let Z be a point on O near X. On T X O we parameterize Z, and thus O, by its Riemannian Logarithm z = Log X (Z). The coordinates in T X O are labeled with i, j, the coordinates in N X O are labeled with a, b. Locally around X, and equivalently around z = 0, we have the graph equation of O:

G a (z) = - 1 2 h a ij (0)z i z j + O(|z| 3 ) (1)
where we use Einstein summation convention: there is an implicit sum on indices that are repeated up and down. The second fundamental form h at 0 is thus the best approximation of O around X as the graph of a quadratic function. It quantifies the external curvature of O in M.

Here we have considered a fixed orbit O. But later in the paper we will consider different orbits. As the space of orbits is by definition the quotient space Q, we will label the orbits by the corresponding point in Q. So we will write O y for y ∈ Q. We will thus have expressions such as h a ij (y, 0), for the second fundamental form of the orbit represented by y.

Geometric interpretation of the template estimation

Here we interpret the statistical template estimation in Computational Anatomy with the quotient space framework. Such an estimation is usually parametric. It means that we assume a parametric model that has produced the observations, in our case: the anatomies. The (parametric) statistical estimation amounts to infer the parameters of the model from the observations, in our case: infer the template, or the mean anatomical shape, from the observed anatomies.

More precisely, one usually assumes that there is an underlying probability density of anatomical shapes, i.e. a probability density on Q, whose mean is called the template. We consider here the simplest generative model with a dirac distribution in Q, i.e. no variability in shapes y ∼ δ(y 0 ). Here y 0 is the template. As we do not observe the anatomical shapes in Q, but rather the anatomies in M, one has to model the poses of the shape. One usually assumes that there is a probability density in the orbits. We assume here a uniform distribution on the orbits, i.e. maximal variability in poses: z ∼ U(O). Finally, we model the noise on the observations by a Riemannian Gaussian in M with isotropic noise σ: x ∼ N M ((y, z), σ) [START_REF] Pennec | Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements[END_REF]. As we have an isometric action and a isotropic noise, this whole generative model is equivalent to the even simpler one: y ∼ δ(y 0 ) and z ∼ δ(z 0 ) and x ∼ N M ((y 0 , z 0 ), σ). Figure 2 illustrates the geometric picture for the generative model. The green dot in Q = R + is the template, or the mean shape. The green circle in M = R 2 is the orbit of y 0 . We choose the pose z 0 on this orbit to get the anatomical object represented by the black dot x 0 = (y 0 , z 0 ) on the green circle. The observed anatomical objects are noisy observations of (y 0 , z 0 ). They are represented by the black triangles generated by the bivariate Gaussian of isotropic noise σ.

y ∼ δ(y0) z ∼ δ(z0) x ∼ N ((y0, z0), σ) 2σ y0 On Q On M Noisy observations on M (y0, z0)
The goal of the statistical estimation is to compute an estimator ŷ0 of the template y 0 . A standard method is the max-max algorithm with unimodal approximation, whose two steps can be geometrically interpreted as follows. First, the registration step amounts to the projection of the observations in Q. Then, one computes the Fréchet mean ŷ0 (a generalization of the notion of mean in Riemannian manifolds) of the observations in Q in order to estimate y 0 .

Figure 3 illustrates the geometric picture for the estimation. First, the observations are registered: the triangles follows the black curve on their orbits (blue circles) in order to be "projected" on Q = R + . Ultimately, the Fréchet mean is computed from the projected triangles on Q = R + .

Registration of the observations

Fréchet mean of the registered observations The algorithm produces an estimator ŷ0 of y 0 . How good is this estimator? To answer this question, we can study its unbiasedness and its consistency. We generalize the usual definitions on linear spaces to Riemannian manifolds. Thus, the bias of an estimator ŷ0 , relative to the parameter y it is designed to estimate, is defined as: Bias y0 [ ŷ0 ] = Log y0 (E y0 [ ŷ0 ]), where Log is the Riemannian logarithm. An unbiased estimator has null bias. A consistent estimator converges in probability to the estimator is designed to estimate, as the number of observations increases. We show that ŷ0 is biased and inconsistent as an estimator of y 0 .

From M to Q Biased statistics on Q y0 ŷ0

Geometric foundations of bias and inconsistency

We consider observations in M generated with the probability density F of our model. We show that they are equivalent to observations in Q generated with a probability density f which we compute.

The probability density on M writes as follows:

F (x) = 1 C exp(- d 2 M (x,x0) 2σ 2
) where C is the normalizing constant. We write O y for the orbit corresponding to y ∈ Q. For an isometric action we have: dM(x) = dO y (x)dQ(y). Thus the probability of having an observation x projecting within the interval of quotient coordinates [y, y + dy] is: P(y ∈ [y, y + dy]) = y+dy y Oy F (x)dO y (x) dQ(y). By taking dy → 0, we have the induced probability density of the quotient space:

f (y) = Oy F (x)dO y (x) = 1 c Oy exp(- d M (x, x 0 ) 2 2σ 2 )dO y (x) (2)
To compute an approximation of the density f , we make the assumption that the noise is low, i.e. σ << 1. The low noise assumption enables to use Taylor expansions in the orbit coordinates and then integrate in O y locally on the ball B y = B(x 0 , σ) ∩ O y ⊂ O y . To this aim, we introduce momenta defined for the uniform measure µ y (z) of the ball B y as:

M i1...ip Oy,p = Oy z i1 ...z ip µ(dz) = By z i1 ...z ip dz (3)
We get in coordinates centered at y 0 ∈ Q:

f (y) = 1 2σ 2 exp(- y a y a 2σ 2 ) M Oy,0 + S(y) ij M ij Oy,2 + O(|M Oy,3 |) (4)
where: S ij (y) = -δ ij + y a h a ij (y, 0) + 1 2 h a ki (y, 0)h a kj (y, 0). The probability density f on Q differs from a normal distribution because of the y-dependent term: M Oy,0 + S(y) ij M ij Oy,2 + O(|M Oy,3 |). This term can be interpreted as a Taylor expansion of the bias with respect to the local geometry of the orbits. The first order of the bias is M Oy,0 and corresponds to the area of B y , i.e. the area of O y seen at the scale of the noise. The second order of the bias is S(y) ij M ij Oy,2 and we recognize a contraction of the second momentum with the matrix S(y) that depends on the external curvature of O y . We expect the higher order terms to be also such contractions between higher order momenta and higher order derivatives of the external curvature of the orbits.

The expectation of f on the ball B(y 0 , σ) ⊂ Q computed at the tangent space of the mean shape T y0 Q gives: Log y0 ( ŷ0 ) a = B(y0,σ) y a f (y)dQ(y). We recognize the bias of the estimator ŷ0 in the case of an infinite number of observations. It differs from 0 because the f distribution on Q is not symmetric. Because we are in the case of an infinite number of observations, this also shows the inconsistency of the estimator ŷ0 . Given the expression of f , we see that bias and inconsistency depend on the external curvature of the orbits and its first derivatives at the scale of the noise. As the external curvature of orbits generally increases when we approach singularities [START_REF] Lytchak | Curvature explosion in quotients and applications[END_REF], the nearer we are from a singularity of Q, the larger is the bias. All in all, when one performs usual statistics on Q from observations on M, the singularities in Q induce bias and inconsistency.

4 An illustration on the quotient Q = R m /SO(m)

We consider the case R m quotiented by SO(m), which is a common example of a finite dimensional flat manifold quotiented by an isometric Lie group action. We perform the computations globally without the low noise assumption.

Generating observations in M = R m with a multivariate normal law and then projecting to Q = R + is equivalent to generating observations directly in Q with the following probability density:

f m (y) = 2 1-m/2
σ m y m-1 exp(-

y 2 + y 0 2 2σ 2 ) 0 F1 ( m 2 , y 2 y 0 2 4σ 2 ) ( 5
)
where 0 F1 is a regularized hypergeometric function. Figure 4 shows that an increase in the noise σ induces that the expectation of f is shifted away from y 0 . The expectation is precisely the estimator ŷ0 for the case of an infinite number of observations. Thus we see the bias increasing with the noise level. We also note that the probability density and therefore the bias depend on the dimension m. The expectation of the probability density f m on the quotient R + is E m (α) = y 0 .e m (α) where

α 2 = y0 2 2σ 2 and e m (α) = 1 α Γ ( m+1 2 ) 1 F1 (-1 2 , m 2
, -α 2 ). Because R + is linear, the bias writes:

Bias y0 ( ŷ0 ) = E[ ŷ0 ] -y 0 = y 0 .(e m (α) -1) = y 0 .bias m (α) (6)
The function bias m (α) is more precisely the bias of any estimator ŷ0 in units of the parameter y 0 it is designed to estimate. It entirely depends in the variable α. First, α can be seen as a signal over noise ratio (SNR): we interpret y 0 as the signal we seek to recover. But we can also interpret α geometrically: it is the ratio of the distance y 0 to the singularity 0 of Q = R + over the noise σ.

The most favorable conditions are when there is no noise with respect to the signal, or equivalently when we are far away from the singularity at the scale of the noise: α → ∞. In this case, we have: bias m (α) → α→∞ 0. In contrast, the less favorable conditions are when the noise is preponderant with respect to the signal, or equivalently when we are close to the singularity at the scale of the noise: α → 0. In this case, we have: bias m (α) →

Conclusion

There is bias in statistics on quotient spaces Q = M/G for noisy observations in M. For a finite dimensional flat Riemannian manifold M, local computations at the scale of noise show that the bias depends on geometric properties of Q, more precisely on its singularities. Global computations on R m quotiented by SO(m) further emphasize that bias cannot be neglected as soon as signal and noise are of the same order, or equivalently as soon as we are close to a singularity at the scale of the noise. Additionally, the increase of the bias with the dimension leads to think that the same phenomenon exists in infinite dimension. Further developments will involve computations for non flat manifolds together with an algorithm to correct the bias. Ultimately, one should generalize the study to the infinite dimensional case.

Fig. 1 .

 1 Fig. 1. Stratification and foliation of M and stratification of Q for M = R 2 quotiented by the action of SO(2) so that Q = R 2 /SO(2) = R+.

Fig. 2 .

 2 Fig. 2. Geometric interpretation of the three steps of the generative model.

Fig. 3 .

 3 Fig. 3. Geometric interpretation of the two steps of the statistical estimation.

Fig. 4 .

 4 Fig.[START_REF] Grenander | Computational anatomy: An emerging discipline[END_REF]. Induced distributions on Q = R+ for m = 3 and noise level σ = 0.3 (left) and σ = 3 (right). In green, the original mean. In red, the estimator of the mean.

1 α 1 Fig. 5 .

 115 Fig. 5. Bias function of the SNR α (left). Influence of the dimension on the bias (right).


Xavier Pennec
email: xavier.pennec@inria.fr

Barycentric Subspaces Analysis on Spheres

This paper addresses the generalization of Principal Component Analysis (PCA) to Riemannian manifolds. Current methods like Principal Geodesic Analysis (PGA) and Geodesic PCA (GPCA) minimize the distance to a "Geodesic subspace". This allows to build sequences of nested subspaces which are consistent with a forward component analysis approach. However, these methods cannot easily be adapted to a backward analysis and alck symmetric in the parametrization of the subspaces. We propose in [10] a new and more general type of family of subspaces in manifolds, barycentric subspaces, which are implicitly defined as the locus of points which are weighted means of k + 1 reference points. Depending on the generalization of the mean that we use, we obtain the Fréchet / Karcher / Exponential barycentric subspaces (FBS / KBS / EBS). The completion of the last one is called the affine span. These definitions were shown to define locally submanifolds of dimension k.

In this paper, we investigate barycentric subspaces in one of the simplest manifold: the sphere. We show that the affine span is a great subsphere in generic conditions, i.e. also a geodesic subspace. This coincidence of spaces is due to the very high symmetry of the sphere. For second order jets, we show that we obtain subspheres of different radii as in the case of principal nested spheres (PNS) analysis. Among the points of the affine span, determining which ones belong to the Karcher barycentric subspaces (KBS) turns out to be a surprisingly difficult problem. Practical experiments show that the KBS covers in general only a small portion of the subsphere containing the reference points. This suggests that the affine span might be a much more interesting definition to work with for subspace definition purposes. We finally discuss the use to these barycentric subspaces to generalize PCA on manifolds, a procedure that we name Barycentric Subpsace Analysis (BSA). Like PGA, barycentric subspaces allow the construction of a forward nested sequence of subspaces which contains the Fréchet mean. However, the definition also allows the construction of backward nested sequence which may not contain the mean.

Introduction

In a Euclidean space, the principal k-dimensional affine subspace of the Principal Component Analysis (PCA) procedure is equivalently defined by minimizing the variance of the residuals (the projection of the data point to the subspace) or by maximizing the explained variance within that affine subspace. This is due to Pythagoras' theorem, which does not hold in more general manifolds. A second important observation is that principal components or different orders are nested, which allows to build forward and backward estimation methods by iteratively adding or removing principal components.

Generalizing PCA to manifolds first requires to define the equivalent of affine subspaces in manifolds. For the zero-dimensional subspace, intrinsic generalization of the mean on manifolds naturally comes into mind: the Fréchet mean is the set of global minima of the variance, as defined by Fréchet in general metric spaces [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF]. The set of local minima of the variance was named Karcher mean by W.S Kendall [START_REF] Kendall | Probability, convexity, and harmonic maps with small image I: uniqueness and fine existence[END_REF] after the work of Karcher et al. [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] on Riemannian centers of mass (see [START_REF] Karcher | Riemannian Center of Mass and so called karcher mean[END_REF] for a discussion of the naming and earlier works). The one-dimensional component is then quite naturally a geodesic which should pass through the mean point. Higher-order components are more difficult to define. The simplest generalization is tangent PCA (tPCA), which amounts to unfold the whole distribution in the tangent space at the mean, and to compute the principal components of the covariance matrix in the tangent space. The method is thus based on the maximization of the explained variance. tPCA was used implicitly or explicitly in a lot of statistical works on shape spaces and Riemannian manifolds because it is simple and efficient. However, if tPCA is good for analyzing data which are sufficiently centered around a central value (unimodal or Gaussian-like data), it is often not sufficient for multimodal or large support distributions (e.g. uniform on close compact subspaces).

Instead of an analysis of the covariance matrix, Fletcher et al. [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF] proposed to rely on the least square distance to subspaces which are totally geodesic at a point. These Geodesic Subspaces (GS) are spanned by the geodesics going through a point with tangent vector restricted to belong to a linear subspace of the tangent space. The procedure was coined Principal Geodesic Analysis (PGA). However, the least-square procedure was computationally expensive, so that it was approximated in practice with tPCA. A complete implementation of the original PGA procedure was only provided recently by Sommer et al. [START_REF] Sommer | Optimization over geodesics for exact principal geodesic analysis[END_REF]. PGA is intrinsic and allows to build a flag (sequences of embedded subspaces) of principal geodesic subspaces which is consistent with a forward component analysis approach: we build iteratively the components from dimension 0 (the mean point), dimension 1 (a geodesic) and higher dimensions by iteratively selecting the direction in the tangent space at the mean that optimally reduce the square distance of data point to the geodesic subspace. In this procedure, the mean always belong to geodesic subspaces even when they are not part of the support of the distribution.

To alleviate this problem, Huckemann et al. [START_REF] Huckemann | Intrinsic shape analysis: Geodesic principal component analysis for Riemannian manifolds modulo Lie group actions[END_REF] proposed to relax the fact that the base-point of the geodesic subspace has to be the Fréchet mean: they start at the first order component directly with the geodesic that best fits the data, which is not necessarily going through the mean. The second principal geodesic is chosen orthogonally to the first one, and higher order components are added orthogonally at the crossing point of the first two components. The method was named Geodesic PCA (GPCA). Further relaxing the assumption that second and higher order components should cross at a single point, Sommers [START_REF] Sommer | Horizontal Dimensionality Reduction and Iterated Frame Bundle Development[END_REF] proposed to parallel transport the second direction along the first principal geodesic to define the second coordinates, and iteratively define higher order coordinates through horizontal development along the previous modes.

All the cited methods are intrinsically forward methods that build successively larger and larger approximation spaces. A notable exception is the concept of Principal Nested Spheres (PNS), proposed by Jung, et al. [START_REF] Jung | Analysis of principal nested spheres[END_REF] as a general framework for non-geodesic decomposition of high-dimensional spheres or highdimensional planar landmarks shape spaces. Here, subsphere can be viewed as a slicing of a higher dimensional sphere by an affine hyperplane. In this process, the nested subsphere is not of radius one, unless the hyperplane is passing through the origin. The backward analysis approach, determining a decreasing family of subspace, has been recently generalized to more general manifold with the help of a nested sequence of relations [START_REF] Damon | Backwards Principal Component Analysis and Principal Nested Relations[END_REF]. However, up to know, such a sequence of relationships was only known for spheres or Euclidean spaces.

In [START_REF] Pennec | Barycentric Subspaces and Affine Spans in Manifolds[END_REF], we proposed to replace geodesic subspaces with a new and more general type of family of subspaces in manifolds: barycentric subspaces (BS). They are implicitly defined as the locus of points which are weighted means of k + 1 reference points. Depending on the generalization of the mean that we use on manifolds, Fréchet mean, Karcher mean or exponential barycenter, we obtain the Fréchet / Karcher / Exponential barycentric subspaces (FBS / KBS / EBS). These definition were shown to be included into each other. Here we call affine span the metric completion of the largest barycentric subspace. In generic conditions, barycentric subspaces are stratified spaces that are locally submanifolds of dimension k. Their singular set of dimension k -l corresponds to the case where l of the reference point belongs to the barycentric subspace defined by the k -l other reference points. In non-generic conditions, points may coalesce along certain directions, defining non local jets instead of a regular k + 1-tuple. Geodesic subspaces (in a restricted sense), which are defined by k tangent vectors at a point, correspond to the limit of the affine span when the k-tuple converges towards that jet.

In this paper, we derive the equations of barycentric subspaces in one of the simplest manifold: the sphere. We show that the affine span of k + 1 different reference points on the n-dimensional sphere is the k-dimensional great subsphere that contains the reference points. In fact, any k + 1-tuple of points of a great k-dimensional subsphere generates the same affine span, which is also a geodesic subspace. This coincidence of spaces is due to the very high symmetry of the sphere. For second order jets, we show that we obtain subspheres of different radii as used in PNS. Among the points of the affine span, determining which ones belong to the Karcher barycentric subspaces (KBS) turns out to be a surprisingly difficult problem. Practical experiments show that the index of the variance at critical points can be arbitrary, thus subdividing the affine span into many regions. As a result, the KBS covers in general only a small portion of the subsphere containing the reference points. This suggests that the affine span might be a much more interesting definition for subspace definition pur-poses. Finally, we discuss the use to these barycentric subspaces to generalize PCA on manifolds, a procedure that we name Barycentric Subspace Analysis (BSA). Barycentric subspaces can be naturally nested, by defining an ordering of the reference points. Like for PGA, this allows the construction of a forward nested sequence of subspaces which contains the Fréchet mean. However, BSA also allows the construction of backward nested sequence which may not contain the mean.

Introduction

In this section, we recall the main notations and results from of [START_REF] Pennec | Barycentric Subspaces and Affine Spans in Manifolds[END_REF]. We consider an embedding Riemannian manifold M of dimension n. The Riemannian metric is denoted . | . x on each tangent space T x M of the manifold. The expression of the the underlying norm in a chart is

v 2 x = v T G(x) v = v i v j g ij (x
) using Einstein notations for tensor contractions. We assume the manifold to be geodesically complete (no boundary nor any singular point that we can reach in a finite time). As an important consequence, the Hopf-Rinow-De Rham theorem states that there always exists at least one minimizing geodesic between any two points of the manifold.

We denote by exp x (v) the exponential map at point x which associate to each tangent vector v ∈ T x M the point of M reached by the geodesic starting at x with this tangent vector after a unit time. This map is a local diffeomorphism from 0 ∈ T x M to M, and we denote -→ xy = log x (y) its inverse: it may be defined as the smallest vector of T x M that allows to shoot a geodesic from x to y. When the tangent space is provided with an orthonormal basis, the Riemannian exponential and logarithmic maps provide a normal coordinate systems at x. A set of normal coordinate systems at each point of the manifold realize an atlas which allows to work very easily on the manifold. The implementation of exp and log maps is the basis of programming on Riemannian manifolds, and most the geometric operations needed for statistics or image processing can be rephrased based on them [START_REF] Pennec | Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements[END_REF][START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF].

(k + 1)-pointed Riemannian manifold

Let {x 0 , . . . x k } ∈ M k+1 be a set of k + 1 distinct points in the Riemannian manifold M and C(x 0 , . . . x k ) = ∪ k i=0 C(x i ) be the union of the cut loci of these points. We call (k + 1)-pointed manifold M * (x 0 , . . . x k ) = M/C(x 0 , . . . x k ) the submanifold of the non-cut points of the points.

Since the cut locus of each point is closed and has null measure, M * (x 0 , . . . x k ) is open and dense in M. Thus, it is a submanifold of M (not necessarily connected). On this submanifold M * (x 0 , . . . x k ), the distance to the points x i and the Riemannian log function -→ xx i = log x (x i ) are smooth.

Weighted moments of a (k + 1)-pointed manifold

Let (λ 0 , . . . λ k ) ∈ R k+1 be weights such that i λ i = 0. We call such weights barycentric coordinates. They are elements of projective space P k minus the orthogonal of the line element 1 = (1 : 1 : . . . 1): P * k = {(λ 0 : . . . : λ k ) ∈ R k+1 such that i λ i = 0}. Standard charts of this space are given either by the intersection of the line elements with the "upper" unit sphere S k of R k+1 with north pole 1/ √ k (unit weights) or by the k-plane of R k+1 passing through the point 1/k and orthogonal to this vector. We call normalized weights λ i = λ i /( k j=0 λ j ) this last projection. Given barycentric coordinates λ ∈ P * k , we can consider the distribution (or 0-current) µ(x) = i λ i δ xi (x) on M. As it is not normalized and weights can be negative, it is generally not a probability distribution. To define the the weighted n-order moment of that distribution, we have to restrict to the (k + 1)-pointed Riemannian manifold M * (x 0 , . . . x k ) because the Riemannian log and distance functions are not defined (resp. smooth) at the cut-locus of the points {x i }:

M n (x, λ) = i λ i -→ xx i ⊗ -→ xx i . . . ⊗ -→ xx i n times (1)
The 0-th order moment (the mass) M 0 (λ) = i λ i = 1 T λ is constant. All other moment are homogeneous of degree 1 in λ and can be normalized by dividing by the mass M 0 (λ). The first order moment M 1 (x, λ) = i λ i -→ xx i is a smooth vector field on the manifold M * (x 0 , . . . x k ). The second and higher order moments are smooth (n, 0) tensor fields that will be used through their contraction with the Riemannian curvature tensor.

Barycentric subspaces of k + 1 points

Let (M, dist) be a metric space and (x 0 , . . . x k ) ∈ M k be k + 1 distinct reference points. The (normalized) weighted variance at point x with weight λ ∈ P * k is:

σ 2 (x, λ) = 1 2 k i=0 λ i dist 2 (x, x i ) = 1 2 k i=0 λ i dist 2 (x, x i )/( k j=0 λ j ).
The Fréchet barycentric subspace is the locus of weighted Fréchet means of these points, i.e. the set of absolute minima of the weighted variance:

F BS(x 0 , . . . x k ) = {arg min x∈M σ 2 (x, λ), λ ∈ P * k }.
The Karcher barycentric subspace KBS(x 0 , . . . x k ) is defined similarly with local minima instead of global ones. A point x ∈ M * (x 0 , . . . x k ) is a weighted exponential barycenters of the reference points if we can find barycentric coordinates λ ∈ P * k such that

M 1 (x, λ) = k i=0 λ i -→ xx i = 0. ( 2
)
The Exponential barycentric subspace EBS(x 0 , . . . x k ) is the set of weighted exponential barycenters of the reference points in M * (x 0 , . . . x k ):

EBS(x 0 , . . . x k ) = {x ∈ M * (x 0 , . . . x k )|∃λ ∈ P * k : M 1 (x, λ) = 0}.
This definition is only valid on M * (x 0 , . . . x k ) and may hide some discontinuities or continuity on the union of the cut locus of the reference points. In order to ensure the completeness of the subspace and potentially reconnect components, we define consider the closure of this set: we call affine span of the points (x 0 , . . . x k ) ∈ M k the closure of the EBS in M:

Aff(x 0 , . . . x k ) = EBS(x 0 , . . . x k ).
Because we assumed that M is geodesically complete, this is equivalent to the metric completion of the EBS.

Outside the cut locus of the reference points, which is of null measure, the gradient of the squared distance d 2 xi (x) = dist 2 (x, x i ) is well defined and is equal to ∇d 2 xi (x) = -2 log x (x i ). Thus, one recognizes that Eq.( 2) defines nothing else than the critical points of the variance σ 2 (x, λ) = 1 2 i λ i dist 2 (x, x i ). The local minima of the variance which are potnetially located on the cut-locus of the reference points are not part of the EBS but they are recovered in the affine span thanks to the metric completion. FBS and KBS are thus included in the affine span, and the affine span is the largest of the barycentric subspaces.

SVD Characterization of the EBS

Let us consider field of n×(k+1

) matrices Z(x) = [ --→ xx 0 , . . . --→ xx k ] on M * (x 0 , . . . x k ).
We can rewrite Eq.( 2) in matrix form: M 1 (x, λ) = Z(x)λ = 0. Thus, we see that the EBS is controled by the kernel of the matrix field Z(x). Let now Z(x) = U (x).S(x).V (x) T be a singular decomposition with singular values sorted in decreasing order. The barycentric subspace EBS(x 0 , . . . x k ) is the zero levelset of the k + 1 singular value s k+1 (x) and the subspace of valid barycentric weights is spanned by the right singular vectors corresponding to the l vanishing singular values: Span(v k-l , . . . v k ) (it is empty if l = 0).

Karcher barycentric subspace and positive span

A critical point of the variance x ∈ EBS(x 0 , . . . x k ) is said non-degenerated (resp. positive) if the Hessian matrix H(x, λ) = -k i=0 λ i D x log x (x i ) is invertible (resp. positive definite) for all λ in the right singular space of the zero singular value of Z(x). The set of degenerate (resp. non-degenerate or positive) exponential barycenter is called the degenerate EBS and denoted EBS 0 (x 0 , . . . , x k ) (resp. non-degenerate EBS * (x 0 , . . . , x k ) or positive EBS + (x 0 , . . . x k )). In Euclidean spaces, all the points are positive and non-degenerated. However, in Riemannian manifolds, we generally have degenerated points and non-degenerated points which are non-positive, as we will see with with the example of spheres. Thus, we can conclude that the KBS is the positive EBS plus potentially some degenerate points of the affine span and some points of the cut locus of the reference points.

Geodesic subspaces as limit case of the affine span

The usual definition of the geodesic subspaces GS(W x ) = {exp x (w) , w ∈ W x } is too large in certain cases to be useful (e.g. in torus). We call restricted geodesic submanifold GS * (W x ) = {exp x (w) , w ∈ W x ∩ D(x)} its restriction to the points that are reached without going through the cut locus of x. This is a well defined submanifold of M whose points are described by homogeneous coordinates at infinity (or on the equator of 1/ √ k depending of the chart we chose for P * k ) of the affine span Aff(x, x 1 , . . . x k ) when the points x i = exp x (ηw i ) are converging to x at first order along the tangent vectors w i defining the k-dimensional subspace W x ⊂ T x M.

Example on spheres

Let us consider the unit sphere as our base manifold. We represent points of M = S n as unit vectors in R n+1 . The tangent space at x is naturally represented by the linear space of vectors orthogonal to x:

T x S n = {v ∈ R n+1 , v T x = 0}.
The natural Riemannian metric on the unit sphere is inherited from the Euclidean metric of the embedding space R n+1 . With these conventions, the Riemannian distance is the arc-length d(x, y) = arccos(x T y). Let us denote f (θ) = 1/sinc(θ) = θ/sin(θ). The spherical exp and log maps are

exp x (v) = cos( v )x + sin( v ) v v = cos( v )x + sinc( v )v log x (y) = θ sin(θ) (y -cos(θ)x) = f (θ) (y -cos(θ)x) with θ = arccos(x T y).
Notice that f (θ) is a smooth function from ] -π; π[ to R that is always greater than one and is locally quadratic at zero:

f (θ) = 1 + θ 2 /6 + O(θ 4 ).
Let us pick k+1 points on the sphere that we put in a matrix X = [x 0 ; . . . x k ]. In the sequel, we use the same notation for the matrix and the set of points. The cut locus of x i is its antipodal point -x i so that M * (X) = S n / -X. Let us denote θ i = arccos(x T i x). The log at a point x is -→

xx i = ( Id -xx T )f (θ i )x i , so that the first weighted moment is M 1 (x, λ) = ( Id -xx T ) i λ i f (θ i )x i = ( Id -xx T )XF (X, x)λ
where F (X, x) = Diag(f (θ i )) is a diagonal matrix with entries that are always greater than one for x ∈ M * (X). Thus the matrix Z(x) = ( Id -xx T )XF (X, x) has the same kernel as the matrix Z(x) = ( Id -xx T )X. This corresponds to a renormalized λ = F (X, x)λ of the weights which is linear in λ but non-linear in x and X through the function F (X, x). The solutions of the equation Z(x) λ = 0 under the constraint x = 1 are given by (x T X λ)x = X λ or more explicitly x = ±X λ/ X λ . This means that the point x ∈ M * (X) has to belong to the Euclidean span of the reference vectors. Conversely, any unit vector x = X λ of the Euclidean span of X verifies Z(x) λ = 0. Thus, the unit vectors x = ±X λ/ λ have barycentric coordinates λ = F (X, x) (-1) λ if x is not at the cut locus of the reference points. This shows that Aff(X) = Span{x 0 , . . . x k } ∩ S n .

(

Following the same principle, we can orthogonalize the reference points x i : let us denote by X = U SV T a singular value decomposition of the matrix of reference vectors. All the singular values s i are positive since the vectors x i are assumed to be linearly independent. Thus, Z(x) has the same kernel as Y (x) = ( Id -xx T )U . This shows that the exponential barycentric subspace generated by the original and orthogonalized points is the same, except at the cut locus of all these points which however belongs to the closure: the affine span. Thus, for spherical data as for Euclidean data, the affine span only depend on the reference points through the point of the Grassmanian that they define.

When the reference points x i are not linearly independent, the matrix X has l ≥ 1 vanishing singular values. A singular value decomposition X = U SV T shows that the value of λ (and thus of λ) is in that case unconstrained in the vector space generated by the right singular vectors v k-l , . . . v k associated to the l vanishing singular values s k-l , . . . s k . Thus, the space of admissible weights at each point of the EBS is of dimension l, and the affine span itself is still the subsphere generated the Euclidean span of the reference points (minus their cut-locus), which is of dimension k -l.

In conclusion, the affine span Aff(X) of k + 1 different reference unit points X = [x 0 ; . . . x k ] on the n-dimensional sphere S n provided with the canonical Euclidean metric of the embedding space R n+1 is the great subsphere of dimension rank(X) -1 that contains the reference points.

Reference points coalescing at order 1

Assume now all the reference points coalesce to a single point x i = exp x0 ( w i ) along the tangent vectors w i which are satisfying x T 0 w i = 0 (to belong to the tangent space at x 0 ) and i w i = 0. This amounts to say that we are following the curve X = X 0 + W in the space of barycentric subspaces, with X T 0 W = 0 and W X T 0 = 0, where here X 0 = x 0 1 T . As previously, the points of the EBS of X are solution of the equation Z(x)λ = 0 and since Z(x) has the same kernel as the matrix Z (x) = ( Idxx T )X , x is solution if we can find some α and λ such that α x = X λ. The additional constraint x 2 = 1 tells us that α 2 = X λ 2 , which gives x = X λ/α when we take the positive root and reintegrate the sign into λ.

In our case, thanks to the orthogonality of X 0 and W , we have

α = X 0 λ + W λ = X 0 λ 2 + 2 W λ 2 ,
Assuming that X 0 λ = 0 (which is in particular the case when X 0 = x 0 1 T since 1 T λ is assumed not to vanish), we have α

-1 = 1 X0λ - 2 2 W λ 2 X0λ 3 + O( 4 ),

Barycentric subspace analysis

We turn in this section to the generalization of principal component analysis itself. PCA can be viewed as the search for a sequence of nested linear spaces that best approximate the data at each level, for instance by minimizing the variance of the residues. In a Euclidean space, this process boils down to an independent optimization of orthogonal subspaces at each level of approximation, thanks to the Pythagorean theorem. This allows to build each subspace of the sequence by adding (resp. subtracting) the optimal one-dimensional subspace iteratively in a forward (resp. backward) analysis. Of course, this property does not scale up to manifolds, for which subspaces have no reason to be orthogonal (even this notion is not well defined).

Flags of barycentric subspaces in manifolds

Nestedness of approximation spaces has been argued to be one of the most important characteristics for generalizing PCA to more general spaces [START_REF] Damon | Backwards Principal Component Analysis and Principal Nested Relations[END_REF]. Barycentric submanifolds can easily be nested, for instance by adding or removing one or several points at a time, which corresponds to put the barycentric weight of this (or these) point(s) to zero. One obtains in that case a family of embedded submanifolds which we call a flag because this generalizes flags of vector spaces. Indeed a flag of a vector space V is a filtration of subspaces (an increasing sequence of subspaces, where "increasing" means each subspace is a proper subspace of the next):

{0} = V 0 ⊂ V 1 ⊂ V 2 ⊂ • • • ⊂ V k = V . Denoting by d i = dim(V i ) the dimension of the subspaces, we have 0 = d 0 < d 1 < d 2 < • • • < d k = n,
where n is the dimension of V. Hence, we must have k ≤ n. A flag is called a complete flag if d i = i, otherwise it is called a partial flag.

With barycentric subspaces of an n-dimensional manifold M, an ordering of n + 1 distinct points x 0 , . . . x n defines a complete flag of barycentric subspaces in the sense that: BS(

x 0 ) = {x 0 } ⊂ • • • BS(x 0 , x 1 , x k ) • • • ⊂ BS(x 0 , . . . x n ).
Grouping points together in the addition/removal process generates a partial flag of barycentric subspaces. Among the barycentric subspaces, the affine span seems to be the most interesting definition to use because the affine span of n + 1 distinct points covers the full manifold: Aff(x 0 , . . . x n ) = M while we only generate a submanifold with the Fréchet or Karcher barycentric subspaces, as we have seen with the example of spheres.

Forward and backward barycentric subspaces analysis

In the classical PCA, the flag of linear subspaces can be built in a forward way, by computing the best 0-th order approximation (the mean), then the best first order approximation (the first mode), etc. It can also be built backward, by removing the direction with the minimal residue from the current affine subspace. In a manifold, we can use similar forward and backward analysis, but they have no reason to give the same result.

With a forward analysis, we compute iteratively the flag of affine spans by adding one point at a time and keeping the previous ones fixed. Thus, we begin by computing the optimal barycentric subspace of dimension 0: Aff(x 0 ) = {x 0 }. Since there is only one weight and it should be unit, the optimal point x 0 found by minimizing the unexplained variance is a Karcher mean. Adding one more point amounts to compute the geodesic passing through the mean that best approximate the data. Adding a second points now differ from PGA, unless the three points coalesce to a single one. The procedure is continued point by point, which mean that the Fréchet mean always belong to the barycentric subspace. In practice, the forward analysis should be stopped when the variance of the residues reaches the noise level of the data, hopefully with k much lower than the dimension n of the embedding manifold, which allows to have an efficient dimension reduction.

The backward analysis consists in iteratively removing one dimension, thus one point in our case. One theoretically should start with a full set of points x 0 , . . . x n which generates the full manifold and chose which one to remove. However, as all the sets of n + 1 distinct points generate the full manifold, the optimization really begin with the set of n points x 0 , . . . x n-1 . Actually this should normally be the only time when we perform an optimization for the point positions, since one should afterward only test for which of the n points we should remove, and this optimization is particularly ill-posed and inefficient in very large dimensional spaces! In order to get around this problem, we propose to run a nonnested forward analysis until we reach the noise level of the data for a dimension k n. Since the goal is only to characterize the optimal k-dimensional subspace, we may optimize the point positions at each step to better fit the data. Then, a backward sweep at the end only reorders the points if necessary by iteratively selecting the one that least increase the unexplained variance. With this process, there is no reason why the Fréchet mean should belong to the reference points (and even to any of the barycentric subspaces). For instance, if we have clusters of points, one expects the reference points to localize within these clusters rather than at the Fréchet mean.

Discussion

We have first investigated in the paper barycentric subspaces in spheres and shown that they encompass both principal geodesic subspaces and nested subspheres as limit cases. It would be interesting to see if we can obtain other types of subspaces with higher order and non-local jets.

The second study point of this paper concerns the generalization of PCA to manifolds using Barycentric Subspace Analysis (BSA). We showed that an ordering of the reference points naturally defines a flag of nested barycentric subspaces. We proposed a first optimization procedure, but the lack of symmetry between the forward and the backward estimations calls for a proper global criterion to be optimized by all k-tuple for k = 0 . . . n together and not just a greedy approach as done by the classical forward and backward approaches.

Other potential practical issues include the fact that the optimization on ktuple might have multiple solutions, as in the case of spheres. Here, we need to find a suitable quotient space similar to the quotient definition of Grassmanians. The optimization might also converge towards a non-local jet instead on a ktuple, and good renormalization techniques need to be designed to guaranty the numerical stability.

so that finally

Thus, we see that the space Af f (x) is the intersection of the sphere with the Euclidean hyperplane generated by X 0 and W , which is once again the geodesic subspace GS(X 0 , W ).

Coalescence at order 2 and link with principal nested spheres

Principal nested spheres were proposed by Jung, Dryden and Marron as a general framework for non-geodesic decomposition of high-dimensional spheres or highdimensional planar landmarks shape spaces [START_REF] Jung | Generalized PCA via the Backward Stepwise Approach in Image Analysis[END_REF][START_REF] Jung | Analysis of principal nested spheres[END_REF]. A subsphere A n-1 of S n is defined as the set of points which are at a fixed distance θ ∈ (0, π/2] of a point x ∈ S n : A n-1 (x, θ) = {y ∈ S n / d(x, y) = θ}. The subsphere A n-1 (x, θ) can be viewed as a slicing of S n by the n-dimensional affine hyperplane P (x, θ) = {y ∈ R n+1 / y T x = cos(θ)}. Notice that the coordinates (x, cos(θ)) of the affine hyperplane parametrize all the possible subspheres of dimension n -1. In this process, the nested subsphere is not of radius one, unless one takes θ = π/2, in which case the hyperplane is passing through the origin.

In order to figure out how nested subspheres and barycentric subspaces are related, we consider the top circle of radius r ∈ [0; 1] around the axis e 3 on the 3-sphere implicitly described by the equation x T e 3 = √ 1 -r 2 . The computations would be exactly the same for the bottom circle x T e 3 = -√ 1 -r 2 . The explicit equation of our circle is: x(ψ) = r cos(ψ)e 1 + r sin(ψ)e 2 + √ 1 -r 2 e 3 . We consider the three points on this cicle at angle ψ = 0, ψ = and ψ = -: X = [x(0), x( ), x(-)]. The spherical affine span consists of the points x = Xλ/ Xλ for λ ∈ S 3 , with:

Thus, the equation x = Xλ can only describe the hyperplane x T e 3 = √ 1 -r 2 when goes to zero (up to a scaling factor s that we can freely choose to be 1 thanks to the homogeneous coordinates), and its intersection with the sphere can only describe a circle of radius r.

Iterating the process, one can generalize the above construction to subspheres of arbitrary dimensions. Thus, we see that Nested Spheres as a limit case of the affine span when the k reference points tend to a 2-jet. It would be interesting to determine which types of subspaces could be obtained by such limits for more general non-local and higher order jets.