Thomas Despois
email: thomas.despois@psemail.eu

Catherine Doz
email: catherine.doz@univ-paris1.fr.

Marie Bessec

Domenico Giannone

Siem Sylvia Kaufmann

Jan Koop- Man

Eric Renault

Mirco Rubin

Takashi Yamagata

Identifying and interpreting the factors in factor models via sparsity: Different approaches

Keywords: factor model, identification, factor interpretation, sparsity, sparse PCA, factor rotation factor model, identification, factor interpretation, sparsity, sparse PCA, factor rotation

With the usual estimation methods of factor models, the estimated factors are notoriously difficult to interpret, unless their interpretation is imposed via restrictions. This paper considers different methods to identify the factor structure and interpret the factors without imposing their interpretation: sparse PCA and factor rotations. We establish a new consistency result for the factors estimated by sparse PCA. Monte Carlo simulations show that our exploratory methods accurately estimate the factor structure, even in small samples. We also apply them on two standard large datasets about international business cycles and the US economy: for each empirical application, they identify the same factor structure, offering a clear economic interpretation of the estimated factors. These exploratory methods can be useful to justify or complement approaches in which the factor structure is imposed a priori.

Introduction

Factor models have been successfully used in various applications, and have become increasingly popular in economics and finance. Despite all their qualities, they suffer from a weakness: with the usual estimation methods, like for instance principal component analysis (PCA), the estimated factors are difficult to interpret, if not impossible. However, identifying the factor structure and interpreting the factors is essential in many economic and financial applications, in particular to allow a structural analysis (e.g. to study international business cycles).

An individual interpretation of the factors can be attempted through the factor loadings. They represent the link between each factor and each variable. A variable with a high loading (in absolute value) on a factor is highly correlated (in absolute value) with this factor. If a large number of loadings of a given factor are small (or even zero), it means that this factor is mainly related to a small number of variables, which generally allows to give an economic interpretation to it. Nevertheless, with the usual estimation methods, the estimated loading matrix tends to have a complex structure with many large loadings for each factor, making the interpretation of the estimated factors difficult, if not impossible.

Most of the literature dealing with the issue of having to identify the factor structure and interpret the factors place restrictions in the loadings in order to impose a specific interpretation to each factor, whether via hierarchical factor models, or via factor rotations with identifying restrictions. As its name suggests, a hierarchical factor model [START_REF] Kose | International business cycles: World, region, and country-specific factors[END_REF] enforces a hierarchy in the factors, e.g. a world factor, several regional factors, and many country-specific factors. The factor rotations with identifying restrictions were introduced in economics by [START_REF] Gürkaynak | Do actions speak louder than words? The response of asset prices to monetary policy actions and statements[END_REF]. The goal was to impose the following interpretation for the two factors of their monetary policy application: surprise changes in the current federal funds rate target, and moves in interest rate expectations over the coming year that are not driven by changes in the current funds rate. Their approach has then been used by several other papers in monetary policy (e.g. [START_REF] Miranda-Agrippino | US monetary policy and the Global Financial Cycle[END_REF]. This methodology was then extended by [START_REF] Bai | Principal components estimation and identification of static factors[END_REF], who proposed two rotations with identifying restrictions such that the first variable is affected by the first factor only, the second variable is affected by the first two factors only (resp. by the second factor only), the third variable is affected by the first three factors only (resp. by the third factor only), etc.

An alternative and more exploratory approach to identify the factor structure and interpret the factors consists in using sparse factor models, introducing sparsity in the factor loadings via purely statistical methods, with very little or no a priori about the factor structure. Sparse factor models have been used for some time in statistics, especially biostatistics (see for instance [START_REF] West | Bayesian factor regression models in the "large p, small n" paradigm[END_REF], among many others), but they have been introduced in economics and finance only very recently. They turned out to help the factor interpretation in these fields too, whether implemented via penalization methods [START_REF] Croux | Sparse and robust factor modelling[END_REF][START_REF] Andreou | Inference in group factor models with an application to mixed frequency data[END_REF][START_REF] Pelger | Interpretable sparse proximate factors for large dimensions[END_REF]Uematsu and Yamagata, 2020a,b), via Bayesian endogenous clustering [START_REF] Francis | An endogenously clustered factor approach to international business cycles[END_REF], or via Bayesian methods using a mixture prior inducing sparsity in the estimated loadings [START_REF] Kaufmann | Bayesian estimation of sparse dynamic factor models with order-independent and ex-post mode identification[END_REF]Beyeler and Kaufmann, forthcoming). On a related issue, sparse factor models can also be used to identify the irrelevant variables for the estimation of the factor model [START_REF] Kaufmann | Identifying relevant and irrelevant variables in sparse factor models[END_REF]. In addition, they have been proved to be able to slightly improve the macroeconomic nowcasting and forecasting performance, for both developed and emerging economies [START_REF] Croux | Sparse and robust factor modelling[END_REF][START_REF] Kristensen | Diffusion indexes with sparse loadings[END_REF][START_REF] Kim | Mining big data using parsimonious factor, machine learning, variable selection and shrinkage methods[END_REF][START_REF] Cepni | Nowcasting and forecasting GDP in emerging markets using global financial and macroeconomic diffusion indexes[END_REF]Uematsu and Yamagata, 2020a).

The main goal of this paper is to compare the performance of conceptually very different methods to implement sparse factor models, allowing to get interpretable factors while remaining deliberately agnostic with respect to the factor structure. We consider two approaches: sparse PCA and factor rotations. Sparse PCA is a machine learning technique which consists in penalizing or constraining the PCA problem so as to induce sparsity in the loadings. We use the SPCA algorithm [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF], which adds an elastic net penalty to the PCA problem. We also consider factor rotations that optimize a simple structure criterion concerning the loading matrix (instead of placing identifying restrictions). These long-established methods, originating from statistics and psychometrics, essentially tend to eliminate the medium loadings (in absolute value), so that each estimated factor tends to have only very large and very small loadings (in absolute value). We use two standard factor rotations: varimax [START_REF] Kaiser | The varimax criterion for analytic rotation in factor analysis[END_REF], which maximizes the variance of squared loadings within each factor, and quartimin [START_REF] Carroll | An analytical solution for approximating simple structure in factor analysis[END_REF][START_REF] Jennrich | Rotation for simple loadings[END_REF], which minimizes the covariance of squared loadings between factors. In addition, we compare the results of these approaches to the findings in [START_REF] Francis | An endogenously clustered factor approach to international business cycles[END_REF] and Beyeler and Kaufmann (forthcoming), who use Bayesian methods. Our approaches are even more agnostic (they do not require a prior or covariate series), are easier to implement, and have a lower computational cost (even sparse PCA). However, the posterior of the Bayesian methods allows to precisely measure the uncertainty in the estimates, e.g. the probability of inclusion of a given variable in each cluster.

The contributions of the paper are both theoretical and empirical. We establish a new consistency result for sparse PCA. Under the standard set of assumptions of approximate dynamic factor models, we establish that sparse PCA consistently estimates the factor space if the 1 penalty parameter is

O(1/ √ N ),
where N is the number of series. Monte Carlo simulations show that sparse PCA and factor rotations accurately estimate the factor structure, even in small samples. The improvement over standard PCA is very substantial when the true structure is sparse. We also study whether these techniques, which are designed to allow to recover forms of sparse structures in the loading matrix, would not fallaciously recover a sparse structure when the true structure is not sparse. It turns out they do not, and can even marginally improve upon PCA, especially when the factors are correlated. The quartimin rotation and sparse PCA (with the SPCA algorithm) seem to be the preferable methods to recover the correct factor structure (at least among the methods considered here). We apply our methods on two standard large datasets about international business cycles and the US economy. For each empirical application, they identify the same factor structure, offering a clear economic interpretation of the estimated factors. The four factors of the application on international business cycles can be interpreted as measures of regional business cycles related to Europe, Latin America, Northern America, and Developed Asia. [START_REF] Francis | An endogenously clustered factor approach to international business cycles[END_REF] get very comparable results with their Bayesian sparse factor model with endogenous clustering. Our eight US macroeconomic factors are related to output, prices, spreads, interest rates, housing, labor, stock market, money and credit. Similarly, Beyeler and Kaufmann (forthcoming) get very comparable results

with their Bayesian sparse FAVAR. This robustness supports the validity of our respective approaches.

It is quite remarkable that such long-established statistical methods as factor rotations, and stateof-the-art machine learning and Bayesian methods, lead to extremely similar results, on two different applications. This suggests that standard estimators like PCA, while computationally convenient, usually recover a quite unnatural factor structure. As soon as an exploratory method allowing for a sparse structure in the loadings is used (no matter the method), it yields the same factor structure. This robustness could also constitute an interesting addition to the results of Giannone, Lenza and Primiceri (forthcoming): there is an illusion of sparsity in many predictive models, but not necessarily in factor models. These exploratory methods provide new lens to study different issues with factor models. In particular, they can be useful to justify or complement approaches in which the factor structure is imposed a priori. For example, most of the literature studying international business cycles with large factor models impose the interpretation of each factor, and that each country is influenced by a single regional business cycle, corresponding to its geographic region. Our methods allow each country to be influenced by several regional business cycles, which proves to be relevant. Several economies appear to be mostly driven by the regional business cycle of another geographic region than their own (e.g. Japan), or to be largely influenced by several regional business cycles (e.g. the United Kingdom).

The rest of the paper is structured as follows. Section 2 introduces the econometric framework.

Section 3 presents the results of Monte Carlo simulations. Section 4 provides two empirical applications.

Section 5 concludes. The online Appendix collects the proof of consistency of the factors estimated by sparse PCA, the detailed presentation of our two datasets, and additional empirical results.

Throughout the paper, the 1 norm of a vector v ∈ R n is denoted by v 1 = n i=1 |v i |, and its

2 norm by v 2 = n i=1 v 2 i .
The 0 norm of v refers to its cardinality, i.e. its number of nonzero elements. For a matrix A = (a ij ) ∈ R m×n , we denote the maximum and minimum eigenvalues of A A by σ max (A A) and σ min (A A), the Frobenius norm of A by

A F = tr(A A) = m i=1 n j=1 a 2
ij , and the induced Euclidean norm of A by A = σ max (A A). The all-ones vector of length n is denoted

1 n .

Econometric framework

Dynamic factor model

We observe N stationary time series over T periods. It is assumed that the observed variables are driven by a fixed number r of unobserved common factors, with r << N . Let X be the T × N panel data matrix (with all the variables standardized beforehand), F the T × r matrix including the r factors, Λ the N × r factor loading matrix, and e the T × N idiosyncratic errors matrix. Denoting i the cross-section unit (i = 1, . . . , N ), t the time unit (t = 1, . . . , T ), and x it the observation of x i at time t, F t = (f 1t , . . . , f rt ) , λ i = (λ i1 , . . . , λ ir ) , the model is written as:

x it = λ i F t + e it , for i = 1, . . . , N, t = 1, . . . , T.

(1)

In matrix notations, with X t = (x 1t , . . . , x N t ) and e t = (e 1t , . . . , e N t ) , it is also written as

X t = ΛF t + e t , for t = 1, . . . , T, (2)
or

X = F Λ + e. (3)
Factor models are characterized by the fact that (F t ) and (e t ) are unobserved stochastic processes with zero mean, which are assumed to be uncorrelated at all leads and lags. Note that although the relation between the observed variables X t and the common factors F t is static (no lag of F t enters equation ( 2)), the model is a dynamic factor model, since (F t ) and (e t ) can be autocorrelated processes.

Further, we assume that it is an approximate factor model, which means that the (e it )'s may be weakly cross-correlated and that the bulk of the comovements of the (x it )'s is captured by the factors. This is obtained through a standard set of assumptions, including in particular the fact that the r eigenvalues of Λ Λ diverge at rate N whereas the covariance matrix of the idiosyncratic errors e t stays bounded when N goes to infinity (see e.g. Stock and Watson, 2002a[START_REF] Bai | Determining the number of factors in approximate factor models[END_REF][START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF], for some standard sets of assumptions).

It is well known that even if the factor space and the common components ΛF t are uniquely identified, F and Λ are identified only up to an invertible matrix. Indeed, for any r × r invertible matrix R, equation (2) can also be written as

X t = Λ F t + e t (4)
with Λ = ΛR and F t = R -1 F t . Similarly, equation (3) can also be written as

X = F Λ + e with F = F (R ) -1 . ( 5
)
Due to this indeterminacy, only the subspace of R r spanned by the factors, and the subspace of R N spanned by the columns of Λ are identifiable and can be estimated. This indeterminacy issue will be the core of our paper. We aim to estimate the loading matrix Λ, which possibly has sparsity properties, in which case the factors should be interpretable. Indeed, the factor loading λ ij measures the influence of the jth factor on the ith variable, so the jth column of Λ can be used to get an interpretation of the jth estimated factor. With the usual estimation methods (e.g. PCA), the estimated loading matrix Λ tends to have a complex structure with many large loadings for each factor, making the interpretation of the estimated factors very difficult, if not impossible. In this paper, we build on the PCA estimator of the model, and propose several approaches allowing to recover a sparse structure in Λ.

Principal component analysis

Stock and Watson (2002a) and [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] proved that the factors in model ( 1) are consistently estimated by principal components under general sets of assumptions associated with approximate dynamic factor models. PCA can be used to estimate model (1) by solving the following nonlinear least squares problem:

min F,Λ N i=1 T t=1 (x it -λ i F t ) 2 = min F,Λ X -F Λ 2 F . (6)
It is well known that the factors obtained as solutions of ( 6) are rescaled versions of the r first principal components of the data covariance matrix, and that the factor loadings obtained as solutions of ( 6) are rescaled versions of the corresponding eigenvectors. It is also well known that this minimization problem is associated to the same indeterminacy issue as the theoretical model: for any r × r

invertible matrix R, if ( F , Λ) is a solution, then ( F , Λ) is also a solution if Λ = ΛR and F = F (R ) -1 .
Normalization constraints have thus to be added in order to get a unique solution for the estimates of and

F P CA = X P D -1/2 (8)
where D = diag( d 1 , . . . , d r ) is the r × r diagonal matrix whose entries correspond to the r largest eigenvalues of the data covariance matrix in decreasing order of magnitude (they are also such that

d k = X p k 2 2 /T ∀k ∈ {1, .
. . , r}), and P = ( p 1 , . . . , p r ) is the N × r matrix containing the corresponding eigenvectors.

Factor rotations

The original method to get interpretable factors is to use factor rotations. 1 The idea is to take advantage of the indeterminacy of the estimated loadings and factors. The initially estimated loading matrix Λ ( Λ P CA here), leading to uninterpretable factors, is multiplied by an r × r invertible matrix R, chosen 1 Factor rotations were originally designed to be applied on the maximum likelihood estimates of exact factor models. When N is large, under the assumptions of approximate factor models, it can be shown that this estimator is asymptotically equivalent to the PCA estimator (for more on this topic, see [START_REF] Chamberlain | Arbitrage, factor structure, and mean-variance analysis on large asset markets[END_REF][START_REF] Chamberlain | Arbitrage, factor structure, and mean-variance analysis on large asset markets[END_REF][START_REF] Schneeweiss | Factor analysis and principal components[END_REF][START_REF] Schneeweiss | Factor analysis and principal components[END_REF]. Applying a factor rotation after a PCA is also valid, and commonplace in many disciplines (see Jolliffe, 2002, chapter 11).

such that the transformed loading matrix optimizes a simple structure criterion:

Λ rotated = ΛR. (9)
The initially estimated factors are then transformed accordingly:

F rotated = F (R ) -1 . ( 10
)
The literature offers numerous rotation methods, which can be either orthogonal or oblique. Orthogonal rotations constrain the factors to stay orthogonal, whereas oblique rotations do not. 2

Even if orthogonal rotations are routinely used, it is highly recommended to apply oblique rotations, at the very least as a starting point for the analysis [START_REF] Fabrigar | Evaluating the use of exploratory factor analysis in psychological research[END_REF][START_REF] Costello | Best practices in exploratory factor analysis: Four recommendations for getting the most from your analysis[END_REF]. Indeed, perfectly uncorrelated factors are quite unrealistic in most fields, and oblique rotations are almost always necessary to recover the correct structure in the loadings (Cattell, 1978, p. 128). Moreover, oblique rotations do not force the correlation of the factors, but rather simply allow for it. If the factors are truly uncorrelated, orthogonal and oblique rotations produce nearly identical results, with the correlations among obliquely rotated factors being close to zero. We chose to consider both orthogonal and oblique rotations to ensure the robustness of our results.

Among orthogonal rotations, varimax [START_REF] Kaiser | The varimax criterion for analytic rotation in factor analysis[END_REF] is generally regarded as the best one. It is overwhelmingly the most widely used rotation. This method consists in finding the rotation matrix R maximizing the variance of squared loadings within each factor: 3

R varimax = arg max R 1 N r k=1 N i=1 ( ΛR) 4 ik - r k=1 1 N N i=1 ( ΛR) 2 ik 2 s.t. R R = I r . (11)
We compute R varimax with the gradient projection algorithm for orthogonal rotations proposed by [START_REF] Jennrich | A simple general procedure for orthogonal rotation[END_REF].

Concerning oblique rotations, the oblimin family is the most popular due to its great flexibility.

Within this family, the quartimin rotation [START_REF] Carroll | An analytical solution for approximating simple structure in factor analysis[END_REF][START_REF] Jennrich | Rotation for simple loadings[END_REF] tends to be recommended because of its simplicity and its results (see Costello and Osborne, 2005, among others).

It searches for the matrix R minimizing the covariance of squared loadings between factors: 4

R quartimin = arg min R j =k N i=1 ( ΛR) 2 ij ( ΛR) 2 ik s.t. diag((R R) -1 ) = 1 r . ( 12
)
2 Though standard in the literature, the term "oblique rotation" is a misnomer. A rotation should preserve angles and distances. Since the factors are orthogonal prior to a rotation, they should remain so afterwards. "Oblique rotations" are invertible transformations which allow the angles between the rotated factors to be different from 90 degrees, hence their name.

3 Hence the name: VARIance MAXimization of squared loadings within each factor. 4 Hence the names: the oblimin rotations all consist in OBLIque rotations via the MINimization of a simplicity criterion, and the quartimin rotation more specifically involves the minimization of a criterion including fourth degree terms.

We compute R quartimin with the gradient projection algorithm for oblique rotations proposed by [START_REF] Jennrich | A simple general method for oblique rotation[END_REF].

The covariance matrix of the rotated estimated factors is

1 T F rotated F rotated = 1 T R -1 F F (R ) -1 = R -1 (R ) -1 for 1 T F F = I r = (R R) -1 . ( 13
)
For quartimin, the constraint diag((R R) -1 ) = 1 r maintains the scale of the estimated factors prior to the rotation. For varimax, the stronger constraint R R = I r also maintains the orthogonality of the estimated factors.

Some points should be clarified about the effects of a rotation. First, a factor rotation does not deteriorate nor improve the fit of the factor model: the estimated factor space, the total variance explained by the estimated factors, and all the commonalities5 remain the same. In particular, since PCA consistently estimates the factor space, this result is not lost after a factor rotation. Secondly, factor rotations like varimax or quartimin achieve near-sparsity rather than sparsity: they do not yield any zero loading, but instead ensure that each estimated factor tends to have only very large and very small loadings in absolute value. Thirdly, when it comes to include the estimated factors in a diffusion index or another model, using the unrotated or rotated estimated factors yields exactly the same predicted values y: a rotation does not change the estimated factor space, so it does not change either the orthogonal projection of the predicted variable y onto this space.

We have tried several other orthogonal and oblique rotations, yielding very similar results. For the sake of brevity, we only consider the standard varimax and quartimin rotations in the following.

Sparse principal component analysis

With sparse PCA, the machine learning literature offers an alternative to factor rotations in order to get interpretable factors, starting with [START_REF] Jolliffe | A modified principal component technique based on the LASSO[END_REF]. The general idea is to penalize or constrain the PCA problem with respect to the norms 0 or 1 so as to induce sparsity in the estimated loadings. The resulting estimated factors are easier to interpret, at the cost of a slightly lower explained variance. This constitutes a trade-off between interpretability and statistical fidelity (explaining most of the variance in the data). However, it is considered that in many applications, the decrease in statistical fidelity implied by the sparse factors is "small and relatively benign" (d'Aspremont, Bach and El Ghaoui, 2008, among others).

Sparse PCA has been successfully used in various fields, especially biostatistics and finance. Indeed, this method offers two key advantages over factor rotations: it yields exact zero loadings (instead of very small loadings in absolute value), 6 and the output can here be governed by one or several hyperparameters, allowing for more flexibility and a finer control. In addition, contrary to the PCA estimates, the factors estimated by sparse PCA are allowed to be correlated, like with oblique rotations.

There has been extensive research in machine learning about sparse PCA, leading to several dozens of competing algorithms (see the survey by [START_REF] Trendafilov | From simple structure to sparse components: a review[END_REF]. We chose to work with the SPCA algorithm [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF], one of the two algorithms recommended by [START_REF] Trendafilov | From simple structure to sparse components: a review[END_REF]. SPCA consists in formulating PCA as a regression-type problem, and adding an elastic net penalty [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], which is a convex combination of the 1 (lasso) and the 2 (ridge) penalties associating their comparative advantages.

The lasso penalty enforces sparsity in the estimated coefficients by performing automatic variable selection. Nevertheless it comes up with limitations, especially when the variables are highly correlated (which is obviously the case with macroeconomic and financial variables). When T > N , if there are high correlations between the variables, the prediction performance of the lasso tends to be dominated by ridge, and thus also by the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Moreover, if there is a group of variables with very high pairwise correlations, then the lasso tends to select only one variable from this group and does not care which one is selected. Using a strictly convex penalty instead of the lasso (e.g. by adding a ridge penalty) overcomes this problem by encouraging what is called the grouping effect in machine learning: the regression coefficients of a group of highly correlated variables tend to be equal (up to a change of sign if negatively correlated). We consider the grouping effect as crucial for a macroeconomic factor model, since it precisely considers groups of highly correlated variables (economic activity variables, prices, interest rates, etc.). SPCA is one of the very few sparse PCA algorithms using a strictly convex penalty.

It contains an elastic net penalty in order to produce modified principal components, possibly sparse. This leads to the following regression-type optimization problem:

7 min A,B 1 T X -XBA 2 F + r k=1 κ 1k b k 1 + κ 2 r k=1 b k 2 2 s.t. A A = I r , B = (b 1 , . . . , b r ). (14)
It is worth emphasizing the links between this minimization problem and PCA. First, elementary calculations show that the solution of the following problem:

min P 1 T X -XP P 2 F s.t. P P = I r ( 15
)
is P , the matrix of the r normalized eigenvectors associated to the r largest eigenvalues of X X/T .

Thus equation ( 14) is a generalization of (15) where P P has been replaced by BA , with the addition of an elastic net penalty. Secondly, if the penalty terms are removed in (14), i.e. if κ 1k = 0, ∀k = 1, . . . , r and κ 2 = 0, then the solution of ( 14) is A = B = P . 8 Thus, even if it is formulated in a different way (due to the presence of the two matrices A and B, which facilitates the construction of the estimation algorithm), ( 14) is indeed a penalized version of the PCA minimization problem (15). The kth column of B, denoted b k , is used to build the kth modified principal component, κ 1k is the 1 penalty parameter controlling the level of sparsity in b k , and κ 2 is the 2 penalty parameter. If desired, the κ 1k hyperparameters can be distinct so as to impose a substantially different level of sparsity for each modified principal component.

SPCA algorithm and estimators. In order to be self-contained, we present the main steps of the SPCA algorithm, for given values of the hyperparameters κ (whose tuning is explained below). We refer the reader to [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF] for a complete presentation.

The algorithm is initialized by setting A (1) = P , where P is the N × r matrix whose columns are the normalized eigenvectors associated to the r largest eigenvalues of the empirical covariance (or correlation) matrix of the data, in decreasing order. It satisfies A (1) A (1) = P P = I r .

Then, if A (k) is the current value of the A matrix, the kth iteration of the algorithm consists in the following two steps.

•

Step 1 (elastic net regression)

If A (k) = a (k) 1 , . . . , a (k) r
, compute

B (k) = b (k) 1 , . . . , b (k) r
such that:

B (k) = argmin B    1 T X -XBA (k) 2 F + r j=1 κ 1j b j 1 + κ 2 r j=1 b j 2 2    .
It can be shown (see [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF] for details) that this amounts to solve separately the r minimization problems:

b (k) j = argmin b 1 T Xa (k) j -Xb 2 2 + κ 1j b 1 + κ 2 b 2 2 .
• Step 2 (reduced rank Procustes rotation) k+1) is defined by:

Given B (k) , A (
A (k+1) = argmin A/A A=Ir 1 T X -XB (k) A 2 F .
It can be shown (see Theorem 4 in [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF] that A (k+1) can be obtained in the following way: if the SVD of

X X T B (k) is given by X X T B (k) = U k ∆ k V k , then A (k+1) = U k V k .
These two steps are iterated until convergence.9

In [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF], denoting B the matrix obtained from the algorithm, the loadings are obtained by the normalization of each column of B i.e. b j is replaced by b j / b j 2 . We choose to adopt a different normalization, in order to be consistent with our PCA normalization, scaling

A = P .
the estimated factors such that diag( F SP CA F SP CA /T ) = 1 r . The estimated loadings and factors are then obtained with

Λ SP CA = B D 1/2 and F SP CA = X B D -1/2 , where d k = X b k 2 2 /T ∀k ∈ {1, . . . , r}, with B and Λ sparse.
Tuning the SPCA hyperparameters κ. First of all, we set all the 1 penalization parameters equal: κ 1k = κ 1 ∀k ∈ {1, . . . , r}. This leaves us with the tuning of only two hyperparameters (κ 1 and κ 2 ) instead of r + 1, drastically reducing the computational cost of the tuning (which can still be high), and also the risk of overfitting.10 Tuning hyperparameters is always a challenge, especially when the variables are time series. This can be done by cross-validation (CV), or by minimizing an information criterion. The standard CV techniques like k-fold CV or leave-one-out CV, assuming that the subsamples are independent and identically distributed, are invalidated when using time series due to the inherent serial correlation.11 A time series CV (simulated out-of-sample exercise) would be relevant, but it seems preferable to use an information criterion instead. First, SPCA is not only used for forecasting applications. Secondly, the computational time can be much smaller using an information criterion. Thirdly, [START_REF] Smeekes | Macroeconomic forecasting using penalized regression methods[END_REF] showed that tuning the penalization parameters with the BIC criterion rather than with time series CV leads to a better forecasting performance in the context of macroeconomic forecasting using penalized regression methods. Thus, we chose to tune the hyperparameters by minimizing a BIC-type criterion, as Kristensen ( 2017) did with the unique tuning parameter of another sparse PCA algorithm (sPCA-rSVD, [START_REF] Shen | Sparse principal component analysis via regularized low rank matrix approximation[END_REF]:

( κ 1 , κ 2 ) = arg min κ 1 ,κ 2 log 1 N T X -F SP CA (κ 1 , κ 2 ) Λ SP CA (κ 1 , κ 2 ) 2 F + m(κ 1 , κ 2 ) log(N T ) N T ( 16
)
where m is the number of nonzero loadings in Λ SP CA . Obviously, Λ SP CA , F SP CA , and m depend on (κ 1 , κ 2 ). We solve ( 16) via a grid search, by minimizing the BIC-type criterion over the following grid:

κ 1 ∈ {0, 0.1, . . . , 1}, κ 2 ∈ {0, 0.1, . . . , 1}.
Consistency of the estimated factors. To the best of our knowledge, the issue of the consistency of the factors estimated by sparse PCA has only been addressed by [START_REF] Kristensen | Diffusion indexes with sparse loadings[END_REF], who proved that sparse PCA consistently estimates the factor space in the case of an approximate dynamic factor model. 12 However, it is worth noticing that in Kristensen (2017) the 1 penalty parameter is supposed to tend to zero quickly, since it vanishes asymptotically at rate 1/T. Consequently, the resulting estimates may not be really sparse. Here, we prove that, in the approximate dynamic factor model framework, under a set of assumptions ensuring the consistency of the PCA estimator, sparse PCA consistently estimates the factor space if

κ 1 = O(1/ √ N ).
This condition is far less restrictive than the one in [START_REF] Kristensen | Diffusion indexes with sparse loadings[END_REF] since in most economic applications N is of the same order of magnitude than T , or even smaller. Different sets of assumptions have been considered in the literature to define approximate dynamic factor models (see for instance Stock and Watson, 2002a[START_REF] Bai | Determining the number of factors in approximate factor models[END_REF][START_REF] Forni | The generalized dynamic factor model consistency and rates[END_REF][START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF]. Here we use the Doz, [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF] set of assumptions, but we could have obtained similar results with, for instance, the Stock and Watson (2002a) or [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] sets of assumptions. In particular, we use the following assumptions of [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF].

(CR1) lim inf n→∞ 1 n σ min (Λ Λ) > 0. (CR2) lim sup n→∞ 1 n σ max (Λ Λ) is finite. (CR3) lim sup n→∞ h∈Z E(e t e t-h
) is finite, which implies in particular that lim sup n→∞ Σ e is finite.

(CR4) inf n σ min (Σ e ) > 0.

We then obtain the following theorem.

Theorem 1. Under assumptions CR1-CR4, if κ 1 = O(1/ √ N ), if F t is the PCA estimator of F t
and if F t is the SPCA estimator of F t , there exists an invertible matrix M , with M = O P (1) and

M -1 = O P (1) such that: 1 T T t=1 F t -M F t 2 2 = O P 1 N .
Proof. See online Appendix A.

Thus we see that the factors estimated by SPCA asymptotically span the same subspace as the factors estimated by PCA. Since the factors estimated by PCA asymptotically span the true factor space, the same applies for the factors estimated by SPCA. Our current result is a mean result: consistency for any t is left for future research.

3 Monte Carlo simulations

Design

In this section, we run Monte Carlo simulations to assess the performance of the different estimators in recovering the factor structure. The data generating process (DGP) is detailed below.

x it = λ i F t + √ θe it = r j=1 λ ij f jt + √ θe it ( 17
)
F t = Φ F F t-1 + u t , Φ F = φ f I r , u t iid ∼ N (0 r , Σ u ) (18) e t = Φ e e t-1 + v t , Φ e = φ e I r , v t iid ∼ N (0 r , Σ v ) (19)
This DGP can generate serially and cross-sectionally correlated factors, as well as serially and crosssectionally correlated errors. All of the experiments are repeated 1,000 times. We set N = 60 and T = 100, r = 3, φ f = 0.5, and φ e ∈ {0, 0.5}.

The covariance (and correlation) matrix of the factors is defined as

Σ F =     1 ρ 1 ρ 2 ρ 1 1 ρ 3 ρ 2 ρ 3 1     . (20)
We set ρ = (ρ 1 , ρ 2 , ρ 3 ) = {(0, 0, 0), (-0.4, 0.2, 0)}, and deduce

Σ u with Σ u = Σ F -Φ F Σ F Φ F = (1 -φ 2 f )Σ F .
The covariance (and correlation) matrix of the errors is a Toeplitz matrix defined as Σ e = (τ |i-j| ) ij .

We set τ ∈ {0, 0.5}, and deduce

Σ v with Σ v = Σ e -Φ e Σ e Φ e = (1 -φ 2
e )Σ e . The case τ = 0 corresponds to an exact factor model.

Concerning the loadings, we consider two structures: dense and by block.

1. Dense structure in the loadings. The loadings are simply generated according to

λ ij iid ∼ N (0, 1). (21)
2. Block structure in the loadings. Each of the three columns of Λ is partitioned into subvectors of whether large or small (possibly zero) loadings:

Λ = Λ 1 Λ 2 Λ 3 =     Λ 11 Λ 12 Λ 13 Λ 21 Λ 22 Λ 23 Λ 31 Λ 32 Λ 33     , (22)
with Λ ij being vectors of length 20 and

Λ ij iid ∼ N (µ l .1 20 , σ 2 l .I 20 ) if i = j, Λ ij iid ∼ N (µ s .1 20 , σ 2 s .I 20 ) if i = j. (23)
We set µ s = 0, and (σ 2 l , σ 2 s ) = {(0.2, 0), (0.2, 0.2), (0.5, 0.2)}. When σ 2 s = 0, the small loadings are zero, and the sparsity level in Λ is 2/3. In order to be able to compare the results between the dense structure and the block structure, the expected norm of the loadings should be the same in both cases.

To do that, µ l is deduced from µ s , σ 2 l , and σ 2 s such that E(Λ k Λ k /N ) = 1 like with the dense structure. Finally, for both the dense structure and the block structure cases, θ is set such that the common component and the idiosyncratic component both explain half of the total variation in X when neither the factors nor the errors are cross-sectionally correlated.

Before computing the metrics assessing the estimation precision, the estimated loadings and factors need to undergo a three-step post-processing.

1. The columns of Λ are rescaled so that their norm is equal to the expected norm of the simulated loadings:

˙ Λ k = √ N Λ k Λ k 2 ∀k ∈ {1, . . . , r}. (24)
This yields

˙ Λ k ˙ Λ k /N = 1, which is equal to E(Λ k Λ k /N ).
The metrics we use in the following to measure the estimation precision are invariant to the scale of the estimated factors, so we do not rescale them.

2. Since F and Λ are identified up to column permutations, when necessary, the columns of F and ˙ Λ are reordered to correspond to the ordering in F and Λ. The columns of Λ and ˙ Λ (taken in absolute value) with the smallest distance according to the 2 norm are iteratively associated.

The order of the estimated factors is then switched accordingly.

3. Since F and Λ are also identified up to column sign changes, when necessary, column sign changes are operated in the matrix containing the reordered estimated factors so that each simulated factor is positively correlated with its estimation. Column sign changes in the matrix containing the rescaled and reordered loadings are then operated accordingly. The resulting loading matrix is denoted Λ.

For the loadings, the estimation precision is measured with the root mean squared error (RMSE)

and the mean absolute error (MAE).

RM SE = 1 N r N i=1 r j=1 (λ ij -λ ij ) 2 (25) M AE = 1 N r N i=1 r j=1 |λ ij -λ ij | (26)
For the factors, the estimation precision is measured with the correlations among the post-processed estimated factors (denoted ρ), and the R 2 of the multivariate regression of the true factors onto the estimated factors:

R 2 F, F = P F F 2 F F 2 F , (27)
where

P F = F ( F F ) -1
F is the projection matrix onto the space spanned by the estimated factors.

An R 2 closer to one implies a more precise estimation of the space spanned by the true factors. This metric is invariant to the chosen rotation and the estimation post-processing.

Results

Tables 1 and2 summarize the results of the simulations with a dense structure in the loadings. The Tables 3 and4 summarize the results of the simulations with a block structure in the loadings.

PCA is now vastly outperformed by the three other estimation methods concerning the estimation of the structure in the loadings. As with the dense structure case, PCA varimax and PCA quartimin lead to extremely similar results when the factors are uncorrelated, but PCA quartimin now dominates PCA varimax more clearly when the factors are correlated. Rotated PCA outperforms SPCA when the structure is near-sparse, but it is the opposite when the structure is sparse. However, the magnitude of the performance improvement with SPCA when the structure is sparse should be read in the light of the high level of sparsity in the model in this case. Concerning the factors, the four estimation methods still lead to similar R 2 s, close to one. The correlations among factors are accurately recovered with SPCA, though slightly underestimated in absolute value. The underestimation is a bit more substantial with PCA quartimin (especially when the structure is near-sparse instead of sparse), but the results are still satisfactory.

For each DGP and estimation method, the introduction of serial correlation or cross-correlation in the errors only has a marginal influence on the estimation performance (except when adding crosscorrelation in the errors to a DGP with a block structure in the loadings: the R 2 s are then slightly more affected, but not the other metrics). None of the competing estimation methods seem to be more affected than the others. For a given DGP, the introduction of correlation among the factors improves the general performance of PCA quartimin and SPCA relative to standard PCA and PCA varimax, which force the estimated factors to be uncorrelated.

In brief, rotated PCA and SPCA slightly improve upon PCA in recovering the correct factor structure if that structure is dense, and vastly if that structure is sparse or near-sparse. We consider the quartimin rotation and SPCA as the preferable approaches to recover the correct factor structure (at least among the methods considered here).

In our simulations, PCA quartimin always performs as well or better than PCA varimax (and there is usually no reason to impose the estimated factors to be uncorrelated as it is done with orthogonal rotations), and outperforms SPCA when the factor structure is near-sparse. SPCA outperforms PCA quartimin when the factor structure is sparse, but this comes at several expenses. First, the computational cost of SPCA can be considerable, which constitutes a serious inconvenience in some applications, especially forecasting exercises. In our Monte Carlo study, the 1,000 repetitions of the experiment for a given DGP require less than one minute for PCA quartimin, versus several days for SPCA. Another problem for forecasting applications is that SPCA cannot handle missing observations per se, and its computational cost tends to discourage the use of an EM algorithm to treat them like with PCA (Stock and Watson, 2002b). A possible solution could be to first treat the missing observations, e.g. with PCA and an EM algorithm, and then apply SPCA on the resulting balanced panel.

Finally, SPCA implies a small loss of explained variance with respect to PCA (rotated or not). In view of this, the quartimin rotation can still be a good compromise over SPCA when the factor structure is sparse. 4 Empirical applications

International business cycles

For comparability purposes, we consider a dataset very similar to other studies of international business cycles with large factor models, especially [START_REF] Kose | International business cycles: World, region, and country-specific factors[END_REF], [START_REF] Francis | An endogenously clustered factor approach to international business cycles[END_REF], and [START_REF] Kaufmann | Identifying relevant and irrelevant variables in sparse factor models[END_REF]. Our measure of business cycle conditions is the annual growth rate of the real GDP at constant national prices. We include 60 countries from 7 regions: North America, Latin America, Europe, Africa, Emerging Asia, Developed Asia, and Oceania (see online Appendix B.1 for further details). The series are retrieved from the version 9.1 of the Penn World Tables [START_REF] Feenstra | The next generation of the Penn World Table[END_REF]. We use the maximum time span possible with this sample , and standardize the series.

The most standard methods to select the number of factors are the information criteria proposed by [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF]. Their six main information criteria select between 1 and 4 factors for this dataset. We choose to include 4 factors in the model, like Francis, Owyang and Savascin (2017) and [START_REF] Kaufmann | Identifying relevant and irrelevant variables in sparse factor models[END_REF]. For r = 4, the grid search of the SPCA hyperparameters yields

( κ 1 , κ 2 ) = (0.6, 0.8).
As routinely reported with PCA, the estimated factors are very difficult to interpret (the estimated factor loadings are presented in online Appendix B.2). Concerning the first factor, the countries with the largest loading are North American, European, and Developed Asian countries. Thus, it could be interpreted as a factor specific to developed countries. Nevertheless, such an interpretation seems imperfect. Several emerging or developing countries also have a large loading (Brazil, Colombia, Costa Rica, El Salvador, Guatemala, Mexico, and South Africa), while several developed countries only have a medium loading (Australia, Iceland, Korea, Malaysia, Singapore, and Thailand), or even a small loading (Ireland and New Zealand). Similarly, the second factor is rather associated with Latin America, but many countries outside this region also have a quite strong loading in absolute value.

The remaining two estimated factors are impossible to interpret.

On the opposite, the loadings estimated by rotated PCA or SPCA systematically offer a straightforward economic interpretation of the factors (Figure 1). The interpretation tends to be slightly more straightforward with SPCA due to the sparsity in the estimated loadings, but the three methods lead to a very similar pattern of estimated loadings. The largest loadings on the first factor correspond to European countries so the first factor can be interpreted as a measure of the European business cycle.

Japan appears to be strongly associated with it. The second factor is the Latin American business cycle: the largest loadings are associated with many Latin American countries and Mexico. South

Africa also has a quite large loading. A possible explanation for the presence of South Africa in this group is that all these countries are highly dependent on oil prices (whether as export or import). The third factor can be interpreted as the Northern American business cycle: Canada and the USA have a large loading. Several Latin American countries that are closely linked to the economy of the USA, in particular because their economy is or was dollarized, also have a large loading (Costa Rica, El Salvador, Honduras) or a medium loading (Chile, Guatemala). It is worth noting that several other former British Commonwealth countries have a medium loading in absolute value (Australia, Kenya, the United Kingdom, and Zimbabwe). This is especially true with PCA varimax and PCA quartimin.

The fourth factor is clearly the Developed Asian business cycle. Indonesia, one of the richest Emerging Asian countries, also has a large loading. Japan has a small loading with PCA varimax and PCA quartimin, and even a zero loading with SPCA. It is not really surprising to find that Japan is more synchronized with the European economies than with the other Asian economies. In particular, Japan started its development earlier than the Asian Tigers and Tiger Cub Economies.

Despite the high level of sparsity in the loadings estimated by SPCA (74.58%), the loss in explained variance is quite limited: the factor space estimated by PCA (rotated or not) accounts for 45.62% of the data variance, versus 44.25% for the factor space estimated by SPCA. Note. The factors are labeled with their economic interpretation. Since F and Λ are identified up to a column sign change, for each factor, we impose the largest loading in absolute value to be positive. The proportions of nonzero loadings for the factors estimated with SPCA are, respectively: 0.42, 0.27, 0.22, 0.12.

Since PCA varimax, PCA quartimin, and SPCA lead to very similar estimated loadings, they also lead to very similar estimated factors (see online Appendix B.2). The dynamics of the estimated factors are consistent with the economic history (see Figure 2 for an illustration with the factors estimated with SPCA). For example, the impact of the Great Recession appears to have been more pronounced on the Northern American and European business cycles. In addition, the Northern American business cycle accurately captures the early 1980s recession, while the Developed Asian business cycle accurately captures the 1997 Asian financial crisis. The estimated factors are strongly correlated with the real GDP growth rates of some representative countries. For example, the correlation between the SPCA estimated Northern American business cycle and the real GDP growth rate of the USA is 85%. These correlations tend to be higher with the SPCA estimates than with the PCA quartimin estimates, and higher with the PCA quartimin estimates than with the PCA varimax estimates. With PCA quartimin and SPCA, which allow the estimated factors to be correlated contrary to PCA and PCA varimax, the pattern of the correlation matrices of the estimated factors is similar (see online Appendix B.2).

For example, the four estimated factors are all positively correlated. However, the correlations tend to be quite substantially higher in absolute value with SPCA. This is essentially in line with the results of our Monte Carlo simulations with a block structure in the loadings and correlated factors. The two most correlated factors are the Northern American and European business cycles (59% with the SPCA estimates, versus 34% with the PCA quartimin estimates). Figure 2: Factors estimated with SPCA, compared with the real GDP growth rates of some representative countries (international business cycles).

Note. The correlations between the factors estimated with SPCA and the real GDP growth rates of the representative countries are, respectively: 0.94, 0.66, 0.85, 0.82.

These interpretable factors provide new lens to study different issues. For instance, we can use the factors in a variance decomposition like many papers of the literature studying international business cycles with large factor models, in order to measure the relative importance of each regional business cycle in driving national business cycles. developing and emerging economies tend to be less driven by international business cycles. Among the European economies, the countries which are part of the European Union tend to be more synchronized with the European business cycle. Most of the literature studying international business cycles with large factor models impose the interpretation of each factor, and that each country is influenced by a single regional business cycle, corresponding to its geographic region. Our methods allow each country to be influenced by several regional business cycles, which proves to be relevant. Several economies appear to be mostly driven by the regional business cycle of another geographic region than their own, or to be largely influenced by several regional business cycles. As seen before with the loadings, Japan is mostly associated with the European business cycle, while several Latin American economies are largely, or even predominantly, driven by the Northern American business cycle (especially Costa Rica, El Salvador, and Honduras). Mexico is more influenced by the Latin American business cycle than by the Northern American business cycle, and the United Kingdom is more influenced by the Northern American business cycle than by the European business cycle. It is also worth noting that the North American economies seem to be more affected by the European business cycle than the European economies are affected by the Northern American business cycle.

Interestingly, Francis, Owyang and Savascin (2017) get very similar results concerning the interpretation of their four estimated factors and the variance decomposition, even though they use a different model (a Bayesian sparse factor model with endogenous clustering) and a slightly different dataset (their prior requires some covariate series in addition to the annual real GDP growth rates of the 60 countries considered). More precisely, they build their model with one global factor and three cluster factors, without imposing the composition of the clusters. The first factor, though supposed to be global, is mainly related to the Developed Asian countries (see Table 4 in [START_REF] Francis | An endogenously clustered factor approach to international business cycles[END_REF]. It explains more than half of the variance of the real GDP growth rate for Indonesia, Korea, Malaysia, Singapore, and Thailand, and far less for the other countries (including Japan). The first cluster includes most of the countries in Europe and also Japan. The second cluster consists of the USA and former British Commonwealth countries (especially Canada and the United Kingdom).

The third cluster comprises most of the countries in Latin America, and Mexico. This similarity with our results supports the validity of our respective approaches. Furthermore, we also extend and complement their results. For example, their model constrain each country to be influenced by only a single cluster factor (in addition to the global factor and the idiosyncratic factor), while our methods imply that each country is allowed to be influenced by several regional business cycles. Moreover, our techniques are very easy to implement, and have a very low computational cost, even for SPCA here (about 0.01 seconds for the rotations, and about 1 minute for SPCA, including the tuning of the hyperparameters). Note. This table shows for each country the percentage of the variance of its real GDP growth rate explained by each estimated factor, and the commonality. Since the factors estimated by SPCA are correlated, we also report an adjusted measure of the percentage of the variance explained by each estimated factor, controlling for the influence of the other estimated factors. The variance decomposition method is detailed in Appendix A.

US economy

The reference large dataset of US macroeconomic series, in particular to implement factor models, is the so-called Stock-Watson dataset (Stock and Watson, 2002a,b) 2016), then standardized.

The six main [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] information criteria select between 8 and 10 factors. We choose to include 8 factors in the model, like [START_REF] Mccracken | FRED-MD: A monthly database for macroeconomic research[END_REF] and [START_REF] Kristensen | Diffusion indexes with sparse loadings[END_REF], among others.

For r = 8, the grid search of the SPCA hyperparameters yields ( κ 1 , κ 2 ) = (0.1, 0.6).

With PCA, the estimated factors are difficult to interpret (the estimated factor loadings are presented in online Appendix C.2). Only four of them can be given an economic interpretation. Factor 1 can be interpreted as a real economic activity factor (even if the interest rates also have a quite strong loading in absolute value). Factor 4 can be interpreted as an interest rates factor (even if the housing, forex, and stock market variables also seem quite important). Factor 7 is obviously a stock market factor. Factor 8 is a money and credit factor. Some variables outside the money and credit group have a quite strong loading in absolute value, but they actually correspond to determinants of the demand for money and credit: average hourly earnings in the goods-producing, construction, and manufacturing sectors (the three largest loadings in absolute value in the labor market group), real personal consumption expenditures, retail and food services sales (the two largest loadings in absolute value in the consumption, orders and inventories group). The remaining estimated factors are impossible to interpret.

On the opposite, the loadings estimated by rotated PCA or SPCA systematically offer a straightforward economic interpretation of the factors: output, prices, spreads, interest rates, housing, labor, stock market, money and credit. More precisely, the three methods lead to the same pattern of estimated loadings (Figure 3). We can also notice that the loadings of the interest rates, stock market, and money and credit factors estimated with these methods are very close to the loadings of their PCA estimated counterparts (especially for the stock market, and money and credit factors). It should be 13 https://files.stlouisfed.org/files/htdocs/fred-databases/fredmdchanges.pdf 14 https://research.stlouisfed.org/econ/mccracken/fred-databases/ stressed that the 8 interpretable factors do not simply consist in recovering the 8 groups of variables of the dataset, or some subgroups of them. For example, there is neither a consumption, orders and inventories factor, nor a forex factor. Besides, the prices, labor, and money and credit factors appear to have very heterogeneous influences on the variables from the prices, labor, and money and credit groups. Moreover, the money and credit factor involves several variables outside the money and credit group.

Despite the high level of sparsity in the loadings estimated by SPCA (70.33%), the loss in explained variance is absolutely marginal: the factor space estimated by PCA (rotated or not) accounts for 47.95%

of the data variance, versus 47.75% for the factor space estimated by SPCA.

Since PCA varimax, PCA quartimin, and SPCA lead to almost identical estimated loadings, they also lead to almost identical estimated factors (see online Appendix C.2). The dynamics of the estimated factors are consistent with the economic history. For example, the housing factor recovers the slow-building of the last housing bubble, its progressive burst starting in early 2006, and the slow recovery of the housing sector starting in early 2009. The estimated factors are strongly correlated with some representative observed variables. For example, the correlation between the SPCA estimated output factor and the growth rate of the industrial production is 96%. These correlations tend to be higher with the SPCA estimates than with the PCA quartimin estimates, and higher with the PCA quartimin estimates than with the PCA varimax estimates. The correlation matrices of the estimated factors are extremely similar between PCA quartimin and SPCA, more than in the previous empirical application (see online Appendix C.2). The correlations tend to be only slightly higher in absolute value with SPCA. The most correlated factors are the three real economic activity factors: output, labor, and housing. For example the output factor and the labor factor have a correlation of 61% with the SPCA estimates, versus 41% with the PCA quartimin estimates.

Beyeler and Kaufmann (forthcoming) get very similar interpretations for the eight US macroeconomic factors of their baseline model, even though they use a different model (a Bayesian sparse FAVAR), and another database (FRED-QD, a quarterly frequency companion to the monthly frequency FRED-MD, with a different composition). This supports the validity of our respective approaches. Note. The factors are labeled with their economic interpretation. Since F and Λ are identified up to a column sign change, for each factor, we impose the largest loading in absolute value to be positive. The proportions of nonzero loadings for the factors estimated with SPCA are, respectively: 0.36, 0.23, 0.25, 0.27, 0.26, 0.35, 0.29, 0.37.

Conclusions

With the usual estimation methods of factor models, the estimated factors are notoriously difficult to interpret, unless their interpretation is imposed via restrictions. This paper considers different methods to identify the factor structure and interpret the factors without imposing their interpretation: sparse PCA (a machine learning technique) and factor rotations (originating from statistics and psychometrics).

We prove that sparse PCA consistently estimates the factor space if the 1 penalty parameter is O(1/ √ N ), where N is the number of series. Monte Carlo simulations show that our exploratory methods accurately estimate the factor structure, even in small samples. They slightly improve upon PCA in recovering the correct factor structure if that structure is dense, and vastly if that structure is sparse or near-sparse. The quartimin rotation and sparse PCA (with the SPCA algorithm) seem to be the preferable methods to recover the correct factor structure (at least among the methods considered here). In particular, PCA quartimin outperforms SPCA when the structure is near-sparse, while SPCA outperforms PCA quartimin when the structure is sparse (but PCA quartimin can then still be a good compromise). We apply our methods on two standard large datasets about international business cycles and the US economy: for each empirical application, they identify the same factor structure, offering a clear economic interpretation of the estimated factors. These interpretations are very close to the results of Francis, Owyang and Savascin (2017) and Beyeler and Kaufmann (forthcoming).

It is quite remarkable that such long-established statistical methods as factor rotations, and stateof-the-art machine learning and Bayesian methods, lead to extremely similar results, on two different applications. This suggests that standard estimators like PCA, while computationally convenient, usually recover a quite unnatural factor structure. As soon as an exploratory method allowing for a sparse structure in the loadings is used (no matter the method), it yields the same factor structure.

These exploratory methods provide new lens to study different issues with factor models. In particular, they can be useful to justify or complement approaches in which the factor structure is imposed a priori.

• The percentage of data variance explained by the r estimated factors is R 2 X, F

. For PCA, several other computations can be used, for example r k=1 d k / N k=1 d k , where d k is the kth largest eigenvalue of the data covariance matrix X X/T .

• The percentage of the variance of the ith variable explained by the r estimated factors (the

commonality of the ith variable) is R 2 X i , F
, where X i denotes the ith column of X.

• The percentage of the variance of the ith variable explained by the kth estimated factor is R 2

X i , F k ,
where F k denotes the kth column of F . It also corresponds to the squared correlation between the ith variable and the kth estimated factor. Of course, when all the estimated factors are uncorrelated (as with PCA and PCA varimax), R 2

X i , F = r k=1 R 2 X i , F k .
• When the estimated factors are correlated (as allowed by PCA quartimin and SPCA), it can be useful to adjust the above measure. The percentage of the variance of the ith variable explained by the kth estimated factor, controlling for the influence of the other estimated factors on the kth estimated factor, is R 2

X i , F -R 2 X i , F -k
, where F -k denotes F except its kth column. This can be represented by a Venn diagram. It also corresponds to the squared semi-partial correlation between the ith variable and the kth estimated factor, controlling for the influence of the other estimated factors on the kth estimated factor. Of course, when all the estimated factors are uncorrelated (as with PCA and PCA varimax), this is simply the percentage of the variance of the ith variable explained by the kth estimated factor.

It is worth noting that R 2 X, F and R 2 X i , F
are invariant to a factor rotation (orthogonal or oblique).

Indeed, given that F rotated = F (R ) -1 , it is trivial to show that P F = P F rotated .

Identifying and interpreting the factors in factor models via sparsity:

Different approaches

Online Appendix

Thomas Despois and Catherine Doz

A Proof of Theorem 1

Preliminary results and notations.

We use the same set of assumptions as in [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF]. We suppose that the normalization condition is the same as in Doz, [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF], which is innocuous, since F t is defined up to a rotation matrix. In particular, the true value of the loadings matrix Λ, satisfies Λ Λ = D, where D is a diagonal r × r matrix whose terms go to infinity linearly with N .1

PCA notations:

We denote:

• d 1 ≥ • • • ≥ d n the ordered eigenvalues of X X T
• p 1 , . . . , p n an associated orthonormal set of eigenvectors

• P = ( p 1 , . . . , p r ) • D = diag( d 1 , . . . , d r )
Other notations:

• We use the same notations as in the description of the SPCA algorithm.

• For any matrix A, we take A = σ max (A A), where σ max (A A) is the maximum eigenvalue of A A.

For any symmetric positive matrix A, we then have A = σ max (A).

• For any matrix A, we say that

A = O P (f (N )) if A = O P (f (N )).
• For any symmetric matrix A, we denote by σ min (A) the smallest eigenvalue of A.

To prove our theorem, we use properties which have been obtained by [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF] under their set of assumptions or which can be easily derived from properties which they have stated. In particular, we use the following properties:

Preliminary properties (PCA properties): [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF] have shown that (see Lemma 2 in their Appendix):

i) D = O P (N ) ii) D -1 = O P 1 N iii) X X T -1 = O P (1) Proof: i) Doz,
1 N D -D = O 1 N + O P 1 √ T
We can thus write:

D = D + O(1) + O P N √
T . As all the terms of D tend to infinity linearly with N , we thus obtain D = O P (N ).

ii) Doz, [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF] have also shown in the same Lemma that:

N D -1 -D -1 = O 1 N + O P 1 √ T
We can thus write:

D -1 = D -1 + O 1 N 2 + O P 1 N √
T . As all the terms of D tend to infinity linearly with N , we thus obtain [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF] have shown that:

D -1 = O P 1 N . iii) If we denote Σ e = X X T -Λ Λ ,
2 σ ij,e -σ ij,e = O P 1 √ N + O P 1 √ T (A1)
and that the result is uniform w.r.t. (i, j).

By assumption (CR4) we know that σ min (Σ e ) = c > 0.

By Weyl theorem, we know that:

σ min Σ e ≥ σ min Σ e -Σ e + σ min (Σ e ) (A2)
Further:

σ min Σ e -Σ e = min x x=1
x Σ e -Σ e x

= min

x x=1 i,j

x i x j ( σ ij,e -σ ij,e )

≤ min

x x=1 i,j |x i ||x j | | σ ij,e -σ ij,e | ≤ max i,j | σ ij,e -σ ij,e | min x x=1 i,j |x i ||x j | ≤ max i,j | σ ij,e -σ ij,e | by Schwarz inequality = O P 1 √ N + O P 1 √ T
Using equation (A2), we thus have:

σ min Σ e ≥ c + O P 1 √ N + O P 1 √ T .
As σ min X X T ≥ σ min Σ e , we also get:

σ min X X T ≥ c + O P 1 √ N + O P 1 √ T so that: σ max X X T -1 ≤ 1 c + O P 1 √ T and X X T -1 = O P (1).
Before proving our theorem, we now provide a few intermediary results concerning the algorithm, which are stated in the following properties.

Proposition 1:

At iteration k, B (k) can be written as:

B (k) = X X T + κ 2 I N -1 X X T A (k) - κ 1 2 Z (k) (A3)
where Z (k) is a N × r matrix whose general term satisfies |z

(k) ij | ≤ 1.
Proof : We know from [START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF] (see the result which we have recalled in our presentation of the algorithm, Step 1) that B (k) is obtained as

B (k) = b (k) 1 , . . . , b (k) r
where:

∀j = 1, . . . , r, b (k) j = argmin b (a (k) j -b) X X T (a (k) j -b) + κ 1 b 1 + κ 2 b 2 2 (A4) For a given k, let us denote f j (b) = (a (k) j -b) X X T (a (k) j -b) + κ 1 b 1 + κ 2 b 2 2 and b (k) j = argmin b f j (b).
Standard convex analysis results (see [START_REF] Rockafellar | Convex Analysis[END_REF], part VI, section 27) imply that 0 belongs to the sub-

differential of f j in b (k) j , denoted as ∂f j (b (k) j ).
Let us then calculate ∂f j (b) the subdifferential of f j in b, for any b.

If we denote g j (b) = (a

(k) j -b) X X T (a (k) j -b) + κ 2 b 2 2 , g j is differentiable, with gradient ∇g j (b) = -2 X X T a (k) j + 2 X X T b + 2κ 2 b.
Lemma:

If we suppose that A (k) = P M k +O P 1 N where M k is a r ×r matrix which satisfies:

M k M k = I r +O P 1 N , then: i) B (k) = P M k + O P 1 N -κ1 2 X X T + κ 2 I N -1 Z (k) ii) X X T B (k) = P DM k + O P (1)
Proof of the Lemma:

i) Using proposition 1, we know from equation (A3) that:

B (k) = X X T + κ 2 I N -1 X X T A (k) - κ 1 2 Z (k)
with Z (k) is a N × r matrix whose general term satisfies |z

(k) ij | ≤ 1. If A (k) = P M k + O P 1
N , this can be written as:

B (k) = X X T + κ 2 I N -1 X X T P M k + O P 1 N - κ 1 2 X X T + κ 2 I N -1 Z (k) (A10)
or equivalently:

B (k) = X X T + κ 2 I N -1 X X T P M k + X X T + κ 2 I N -1 X X T ×O P 1 N - κ 1 2 X X T + κ 2 I N -1 Z (k) (A11)
Let us now study the two first terms of (A11).

We know that the columns of P are eigenvectors of X X T , associated to the eigenvalues contained in the diagonal matrix D. It follows that the columns of P are also eigenvectors of X X T + κ 2 I N associated to the eigenvalues contained in the matrix D + κ 2 I r , and eigenvectors of ( X X T + κ 2 I N ) -1 associated to the eigenvalues contained in the matrix ( D + κ 2 I r ) -1 . We thus get:

X X T + κ 2 I N -1 X X T P M k = P ( D + κ 2 I r ) -1 DM k (A12)
Further, the diagonal terms of ( D + κ 2 I r ) -1 D are di di+κ2 for i = 1, . . . , r, with

d i d i + κ 2 = 1 1 + κ2 di ∼ 1 - κ 2 d i = 1 + O P 1 N (A13) so that ( D + κ 2 I r ) -1 D = I r + O P 1 N .
As P P = I r , we have P = O P (1). As M k M k = I r + O P 1 N , we also have M k = O P (1). We can then write:

P ( D + κ 2 I r ) -1 DM k = P (I r + O P 1 N )M k = P M k + O P 1 N (A14)
Using (A12), we thus get:

X X T + κ 2 I N -1 X X T P M k = P M k + O P 1 N . (A15)
Turning to the second term of (A11), we can notice that the eigenvalues of

X X T + κ 2 I N -1 X X T are di di+κ2
, i = 1, . . . , n, which are all smaller than 1.

We thus get:

X X T + κ 2 I N -1 X X T ≤ 1, so that this second term is O P 1 N .
Part (i) of the Lemma then follows.

ii) It follows from (i) that:

X X T B (k) = X X T P M k + O P 1 N - κ 1 2 X X T + κ 2 I N -1 Z (k) (A16) = P DM k + X X T × O P 1 N - κ 1 2 X X T X X T + κ 2 I N -1 Z (k) (A17)
As X X T = O P (N ), we thus get:

X X T B (k) = P DM k + O P (1) - κ 1 2 X X T X X T + κ 2 I N -1 Z (k) (A18) Further: • X X T X X T + κ 2 I N -1
has the same eigenvalues than X X T + κ

2 I N -1 X X T so that X X T X X T + κ 2 I N -1 = X X T + κ 2 I N -1 X X T ≤ 1 • Z (k) 2 ≤ tr Z (k) Z (k) = r j=1 N i=1 z (k) ij 2 ≤ rN so that Z (k) = O P ( √ N ).
It then follows that:

X X T X X T + κ 2 I N -1 Z (k) ≤ X X T X X T + κ 2 I N -1 Z (k) = O p ( √ N ) As κ 1 = O 1 √ N , part (ii) of the Lemma follows.

Proof of Proposition 2:

We prove Proposition 2 by mathematical induction.

Using results (i) and (ii) of the lemma, and using the fact that D = O P (N ), P = O P (1) and M k = O P (1), we get that:

B X X T B = P M k + O P 1 N - κ 1 2 X X T + κ 2 I N -1 Z (k) P DM k + O P (1) = M k DM k + O P (1) - κ 1 2 Z (k) X X T + κ 2 I N -1 P DM k - κ 1 2 Z (k) X X T + κ 2 I N -1 × O P (1) = M k DM k + O P (1) - κ 1 2 Z (k) P D + κ 2 I N -1 DM k - κ 1 2 Z (k) X X T + κ 2 I N -1 × O P (1)
We have seen that

Z (k) = O P ( √ N ). As we assume that κ 1 = O 1 √ N , we get κ 1 Z (k) = O P (1).
We also have seen that

D + κ 2 I N -1 D ≤ I r so that D + κ 2 I N -1 D = O P (1)
.

Further X X T + κ 2 I N -1 ≤ X X T -1
= O P (1). We thus get:

B X X T B = M k DM k + O P (1) If we denote M k = m (k) 1 . . . m (k) r , the diagonal terms of M k DM k are m (k) i Dm (k) i = r j=1 d j m (k) ij 2 As M k M k = I r + O P 1 N , we have r j=1 d j m (k) ij 2 ∈ [ d r + O P (1), d 1 + O P (1)
] for i = 1, . . . , r so that all the terms of ∆ diverge linearly with N like the terms of D. We thus have:

∆ -1/2 = O P 1 √ N
As F t = ∆ -1/2 B x t , we can write, using result (i) of the lemma:

F t = ∆ -1/2 M k P + O P 1 N - κ 1 2 Z (k) X X T + κ 2 I N -1
x t and we know that the O P 1 N term comes from equation (A13) and is the same one for any t.

Using the fact that F t = D -1/2 P x t , we then get:

F t = ∆ -1/2 M k D 1/2 F t + O P 1 N -∆ -1/2 κ 1 2 Z (k) X X T + κ 2 I N -1 x t (A25) If we denote M = ∆ -1/2 M k D 1/2
, it is clear that M is invertible and that M = O P (1).

Let us now show that 1

T T t=1 F t -M F t 2 2 = O P 1 N .
We can write:

1 T T t=1 F t -M F t 2 2 = 1 T T t=1 O P 1 N -∆ -1/2 κ1 2 Z (k) X X T + κ 2 I N -1 x t O P 1 N -∆ -1/2 κ1 2 Z (k) X X T + κ 2 I N -1
x t

We have seen that:

. the O P 1 N term does not depend on t

.

Z (k) = O P ( √ N ) . ∆ -1/2 = O P 1 √ N . X X T + κ 2 I N -1 ≤ X X T -1 = O P (1)
We can decompose 1

T T t=1 F t -M F t 2 2
into three terms that we study separately.

• The first term has the form 1

T T t=1 O P 1 N × O P 1 N and the O P 1 N term does not depend on t, so this term is O P 1 N 2 .
• The second term is:

2 × O P 1 N × 1 T T t=1 ∆ -1/2 κ 1 2 Z (k) X X T + κ 2 I N x t (A26)
which can be also written as

κ 1 × O P 1 N × ∆ -1/2 Z (k) X X T + κ 2 I N 1 T T t=1 x t (A27) Since Z (k) = O P ( √ N ), ∆ -1/2 = O P 1 √ N , X X T + κ 2 I N -1
= O P (1), we know that:

κ 1 × O P 1 N × ∆ -1/2 Z (k) X X T + κ 2 I N = O P 1 N √ N . Now 1 T T t=1 x t = Λ 1 T T t=1 F t + 1
T T t=1 e t , with:

-Λ = O( √ N ) -1 T T t=1 F t = O P 1 √ T since (F t ) is a stationary process -1 T T t=1 e t = O P √ N √ T
by assumption (CR3) since:

E   1 √ T T t=1 e t 2   = E 1 T T t=1 h∈Z e t e t-h = 1 T T t=1 h∈Z
tr E e t e t-h ≤ N h∈Z E(e t e t-h ) = O(N ).

We thus get that 1

T T t=1 x t = O P √ N √ T
and, using (A27), we obtain that (A26) is

O P 1 N √ T .
• The third term is:

κ 2 1 4 1 T T t=1 x t X X T + κ 2 I N -1 Z (k) ∆ -1 Z (k) X X T + κ 2 I N -1 x t (A28)
which can also be written as:

κ 2 1 4 1 T T t=1 tr ∆ -1/2 Z (k) X X T + κ 2 I N -1 x t x t X X T + κ 2 I N -1 Z (k) ∆ -1/2 = κ 2 1 4 tr ∆ -1/2 Z (k) X X T + κ 2 I N -1 1 T T t=1 x t x t X X T + κ 2 I N -1 Z (k) ∆ -1/2 = κ 2 1 4 tr ∆ -1/2 Z (k) X X T + κ 2 I N -1 X X T X X T + κ 2 I N -1 Z (k) ∆ -1/2 (A29)
Now, we can write:

X X T + κ 2 I N -1/2 X X T X X T + κ 2 I N -1/2 ≤ I N so that: X X T + κ 2 I N -1 X X T X X T + κ 2 I N -1 ≤ X X T + κ 2 I N -1 ≤ X X T -1 = O P (1) Since κ 1 = O 1 √ N , Z (k) Z (k) = O P (N ) and ∆ -1/2 = O P 1 √
N , we then get:

κ 2 1 4 ∆ -1/2 Z (k) X X T + κ 2 I N -1 X X T X X T + κ 2 I N -1 Z (k) ∆ -1/2 = O P 1 N
As this matrix is r × r, its trace is also

O P 1 N , so that (A29) is O P 1 N .
It then follows that the summation of the three terms of our decomposition is also

O P 1 N i.e. that 1 T T t=1 F t -M F t 2 2 = O P 1 N .
B Supplementary material for the application on international business cycles

B.1 Data

The regional definitions follow [START_REF] Kose | International business cycles: World, region, and country-specific factors[END_REF]. The Developing Asia region is renamed Emerging Asia.

• North America: Canada, Mexico, USA.

• Latin America: Argentina, Bolivia, Brazil, Chile, Colombia, Costa Rica, Dominican Republic, Ecuador, El Salvador, Guatemala, Honduras, Jamaica, Panama, Paraguay, Peru, Trinidad and Tobago, Uruguay, Venezuela.

• Europe: Austria, Belgium, Denmark, Finland, France, Germany, Greece, Iceland, Ireland, Italy, Luxembourg, Netherlands, Norway, Portugal, Spain, Sweden, Switzerland, United Kingdom.

• Africa: Cameroon, Côte d'Ivoire, Kenya, Morocco, Senegal, South Africa, Zimbabwe.

• Emerging Asia: Bangladesh, India, Indonesia, Pakistan, Philippines, Sri Lanka.

• Developed Asia: Hong Kong SAR, Japan, Korea, Malaysia, Singapore, Thailand.

• Oceania: Australia, New Zealand.

B.2 Additional results

Table B1: Selected number of factors (international business cycles). Note. Number of factors selected by each of the six main [START_REF] Bai | Determining the number of factors in approximate factor models[END_REF] information criteria, allowing for at most 4 factors. Note. The interpretable factors are labeled with their economic interpretation, and the uninterpretable factors with a question mark. Since F and Λ are identified up to a column sign change, for each factor, we impose the largest loading in absolute value to be positive. Note. The correlation matrices of the factors estimated with PCA and PCA varimax are not displayed since they are simply the identity matrix, by construction. C Supplementary material for the application on the US economy

C.1 Data

We use the 2018:9 vintage of FRED-MD. The composition of FRED-MD is subject to minor changes over time.

The 128 variables of this vintage are listed below. Most of them are retrieved from the Federal Reserve Economic Database (FRED), so their mnemonic is the same as in FRED. Some series require adjustments to the raw data available in FRED (see [START_REF] Mccracken | FRED-MD: A monthly database for macroeconomic research[END_REF]. The transformation code (TC) indicates how the series was transformed to ensure stationarity: (1) no transformation, (2) ∆x t , (3) ∆ 2 x t , (4) log(x t ), (5) ∆ log(x t ), (6)

∆ 2 log(x t ), (7) ∆(x t /x t-1 -1). In order to work with a balanced panel, we drop 5 of 128 variables. The binary entry BP indicates whether that variable is included in the balanced panel 1960:1-2018:4.

C.2 Additional results

Table C14: Selected number of factors (FRED-MD).

Criterion ICp1 ICp2 ICp3 P Cp1 P Cp2 P Cp3 Note. The correlation matrices of the factors estimated with PCA and PCA varimax are not displayed since they are simply the identity matrix, by construction.

F

  and Λ (up to column permutations and column sign changes). The two standard normalizations are F F /T = I r with Λ Λ diagonal, and Λ Λ/N = I r with F F diagonal. The first normalization sets the scale of the estimated factors, while the second sets the scale of the estimated loadings. Since most factor rotations require the unrotated estimated factors to be such that F F /T = I r , we adopt the first normalization. Under this normalization, the estimates of the loadings and factors are obtained with Λ P CA = P D 1/2(7)

  four estimation methods lead to comparable results. PCA varimax, PCA quartimin and SPCA do not fallaciously recover a sparse or near-sparse structure in the loadings. They even tend to slightly improve upon standard PCA, especially when the factors are correlated. The factor space is precisely estimated: the R 2 s are close to one. With PCA quartimin and SPCA, the correlations among factors are very accurately estimated when they are uncorrelated, but are underestimated in absolute value when they are correlated.

F3Figure 1 :

 1 Figure1: Factor loadings estimated with PCA varimax, PCA quartimin, and SPCA (international business cycles).

F3

Figure 3 :

 3 Figure 3: Factor loadings estimated with PCA varimax, PCA quartimin, and SPCA (FRED-MD).

Figure B1 :

 B1 Figure B1: Grid search of the SPCA hyperparameters κ 1 and κ 2 (international business cycles).

F3Figure B2 :

 B2 FigureB2: Factor loadings estimated with PCA (international business cycles).

F3Figure B3 :

 B3 Figure B3: Factors estimated with PCA varimax, PCA quartimin, and SPCA (international business cycles).

Figure B4 :

 B4 Figure B4: Correlation matrices of the estimated factors (international business cycles).

Figure C5 :

 C5 Figure C5: Grid search of the SPCA hyperparameters κ 1 and κ 2 (FRED-MD).

Figure C7 :Figure C8 :

 C7C8 Figure C7: Factors estimated with SPCA, compared with some representative observed variables (FRED-MD).Note. Shaded areas indicate periods of recession as defined by the National Bureau of EconomicResearch. The correlations between the factors estimated with SPCA and the representative observed variables are, respectively: 0.96, 0.89, 0.95, 0.96, 0.99, 0.91, 0.94, 0.61.

Figure C9 :

 C9 Figure C9: Correlation matrices of the estimated factors (FRED-MD).

Table 1 :

 1 Dense structure in the loadings with uncorrelated factors.

	Estimation method φe	τ	RMSE MAE	ρ	R 2 F, F
	PCA	0	0	0.65	0.52	(0,0,0)	0.90
	PCA varimax	0	0	0.63	0.51	(0,0,0)	0.90
	PCA quartimin	0	0	0.63	0.51	(0,0,0)	0.90
	SPCA	0	0	0.66	0.53	(0,-0.01,0.01)	0.89
	PCA	0.5	0	0.66	0.53	(0,0,0)	0.91
	PCA varimax	0.5	0	0.65	0.52	(0,0,0)	0.91
	PCA quartimin	0.5	0	0.65	0.52	(0,0,0)	0.91
	SPCA	0.5	0	0.66	0.53	(0,0,0)	0.91
	PCA	0	0.5	0.66	0.53	(0,0,0)	0.88
	PCA varimax	0	0.5	0.65	0.52	(0,0,0)	0.88
	PCA quartimin	0	0.5	0.65	0.52	(0,0,0)	0.88
	SPCA	0	0.5	0.67	0.54	(0,0,0)	0.88
	PCA	0.5 0.5	0.67	0.54	(0,0,0)	0.90
	PCA varimax	0.5 0.5	0.65	0.53	(0,0,0)	0.90
	PCA quartimin	0.5 0.5	0.65	0.53	(0,0,0)	0.90
	SPCA	0.5 0.5	0.67	0.54	(0,0,0)	0.90

Table 2 :

 2 Dense structure in the loadings with correlated factors.

	Estimation method φe	τ	RMSE MAE	ρ	R 2 F, F
	PCA	0	0	0.72	0.58	(0,0,0)	0.89
	PCA varimax	0	0	0.67	0.54	(0,0,0)	0.89
	PCA quartimin	0	0	0.65	0.52	(-0.09,0.04,0.01)	0.89
	SPCA	0	0	0.70	0.56	(-0.13,0.07,0.02)	0.89
	PCA	0.5	0	0.73	0.59	(0,0,0)	0.91
	PCA varimax	0.5	0	0.68	0.55	(0,0,0)	0.91
	PCA quartimin	0.5	0	0.67	0.54	(-0.09,0.04,0.01)	0.91
	SPCA	0.5	0	0.71	0.57	(-0.09,0.05,0.01)	0.91
	PCA	0	0.5	0.73	0.59	(0,0,0)	0.88
	PCA varimax	0	0.5	0.68	0.55	(0,0,0)	0.88
	PCA quartimin	0	0.5	0.67	0.54	(-0.09,0.04,0.01)	0.88
	SPCA	0	0.5	0.71	0.57	(-0.12,0.05,0.01)	0.88
	PCA	0.5 0.5	0.74	0.59	(0,0,0)	0.90
	PCA varimax	0.5 0.5	0.69	0.56	(0,0,0)	0.90
	PCA quartimin	0.5 0.5	0.68	0.55	(-0.08,0.04,0.01)	0.90
	SPCA	0.5 0.5	0.72	0.58	(-0.08,0.04,0.01)	0.90

Table 3 :

 3 Block structure in the loadings with uncorrelated factors.

	Estimation method σ 2 l	σ 2 s	φe	τ	RMSE MAE	ρ	R 2 F, F
	PCA	0.2	0	0	0	0.67	0.57	(0,0,0)	0.92
	PCA varimax	0.2	0	0	0	0.23	0.18	(0,0,0)	0.92
	PCA quartimin	0.2	0	0	0	0.22	0.17	(0,0,0)	0.92
	SPCA	0.2	0	0	0	0.15	0.07	(0,0,0)	0.93
	PCA	0.2 0.2	0	0	0.65	0.55	(0,0,0)	0.92
	PCA varimax	0.2 0.2	0	0	0.24	0.20	(0,0,0)	0.92
	PCA quartimin	0.2 0.2	0	0	0.24	0.19	(0,0,0)	0.92
	SPCA	0.2 0.2	0	0	0.30	0.24	(0,0,0)	0.92
	PCA	0.5 0.2	0	0	0.66	0.54	(0,0,0)	0.91
	PCA varimax	0.5 0.2	0	0	0.29	0.23	(0,0,0)	0.91
	PCA quartimin	0.5 0.2	0	0	0.28	0.22	(0,0,0)	0.91
	SPCA	0.5 0.2	0	0	0.32	0.26	(0,0,0)	0.91
	PCA	0.2	0	0.5	0	0.69	0.58	(0,0,0)	0.94
	PCA varimax	0.2	0	0.5	0	0.25	0.20	(0,0,0)	0.94
	PCA quartimin	0.2	0	0.5	0	0.24	0.18	(0,0,0)	0.94
	SPCA	0.2	0	0.5	0	0.16	0.08	(0,0,0)	0.94
	PCA	0.2 0.2 0.5	0	0.67	0.56	(0,0,0)	0.93
	PCA varimax	0.2 0.2 0.5	0	0.26	0.21	(0,0,0)	0.93
	PCA quartimin	0.2 0.2 0.5	0	0.25	0.20	(0,0,0)	0.93
	SPCA	0.2 0.2 0.5	0	0.32	0.26	(0,0,0)	0.93
	PCA	0.5 0.2 0.5	0	0.67	0.56	(0,0,0)	0.93
	PCA varimax	0.5 0.2 0.5	0	0.31	0.24	(0,0,0)	0.93
	PCA quartimin	0.5 0.2 0.5	0	0.31	0.24	(0,0,0)	0.93
	SPCA	0.5 0.2 0.5	0	0.35	0.28	(0,0,0)	0.93
	PCA	0.2	0	0	0.5	0.67	0.58	(0,0,0)	0.85
	PCA varimax	0.2	0	0	0.5	0.24	0.19	(0,0,0)	0.85
	PCA quartimin	0.2	0	0	0.5	0.22	0.17	(0.01,0,0)	0.85
	SPCA	0.2	0	0	0.5	0.15	0.07	(0.01,0,0)	0.85
	PCA	0.2 0.2	0	0.5	0.66	0.55	(0,0,0)	0.85
	PCA varimax	0.2 0.2	0	0.5	0.25	0.20	(0,0,0)	0.85
	PCA quartimin	0.2 0.2	0	0.5	0.24	0.19	(0,0,0)	0.85
	SPCA	0.2 0.2	0	0.5	0.30	0.24	(0,0,0)	0.85
	PCA	0.5 0.2	0	0.5	0.67	0.55	(0,0,0)	0.84
	PCA varimax	0.5 0.2	0	0.5	0.30	0.23	(0,0,0)	0.84
	PCA quartimin	0.5 0.2	0	0.5	0.29	0.23	(0,0,0)	0.84
	SPCA	0.5 0.2	0	0.5	0.33	0.26	(0,0,0.01)	0.84
	PCA	0.2	0	0.5 0.5	0.69	0.58	(0,0,0)	0.88
	PCA varimax	0.2	0	0.5 0.5	0.25	0.20	(0,0,0)	0.88
	PCA quartimin	0.2	0	0.5 0.5	0.24	0.18	(0.01,0,0)	0.88
	SPCA	0.2	0	0.5 0.5	0.17	0.08	(0.01,0,0)	0.88
	PCA	0.2 0.2 0.5 0.5	0.67	0.56	(0,0,0)	0.88
	PCA varimax	0.2 0.2 0.5 0.5	0.26	0.21	(0,0,0)	0.88
	PCA quartimin	0.2 0.2 0.5 0.5	0.25	0.20	(0,0,0)	0.88
	SPCA	0.2 0.2 0.5 0.5	0.33	0.26	(0,0,0)	0.88
	PCA	0.5 0.2 0.5 0.5	0.68	0.56	(0,0,0)	0.87
	PCA varimax	0.5 0.2 0.5 0.5	0.32	0.25	(0,0,0)	0.87
	PCA quartimin	0.5 0.2 0.5 0.5	0.31	0.24	(0,0,0)	0.87
	SPCA	0.5 0.2 0.5 0.5	0.35	0.28	(0,0,0)	0.87

Table 4 :

 4 Block structure in the loadings with correlated factors.

	Estimation method σ 2 l	σ 2 s	φe	τ	RMSE MAE	ρ	R 2 F, F
	PCA	0.2	0	0	0	0.72	0.61	(0,0,0)	0.92
	PCA varimax	0.2	0	0	0	0.29	0.23	(0,0,0)	0.92
	PCA quartimin	0.2	0	0	0	0.23	0.18	(-0.32,0.16,0.01)	0.92
	SPCA	0.2	0	0	0	0.16	0.07	(-0.38,0.18,0)	0.93
	PCA	0.2 0.2	0	0	0.71	0.60	(0,0,0)	0.92
	PCA varimax	0.2 0.2	0	0	0.30	0.24	(0,0,0)	0.92
	PCA quartimin	0.2 0.2	0	0	0.26	0.21	(-0.23,0.11,0.01)	0.92
	SPCA	0.2 0.2	0	0	0.32	0.26	(-0.35,0.18,0.01)	0.92
	PCA	0.5 0.2	0	0	0.72	0.60	(0,0,0)	0.91
	PCA varimax	0.5 0.2	0	0	0.34	0.27	(0,0,0)	0.91
	PCA quartimin	0.5 0.2	0	0	0.30	0.24	(-0.23,0.11,0.01)	0.91
	SPCA	0.5 0.2	0	0	0.34	0.27	(-0.36,0.18,0.01)	0.91
	PCA	0.2	0	0.5	0	0.73	0.62	(0,0,0)	0.94
	PCA varimax	0.2	0	0.5	0	0.31	0.25	(0,0,0)	0.94
	PCA quartimin	0.2	0	0.5	0	0.25	0.20	(-0.32,0.16,0.01)	0.94
	SPCA	0.2	0	0.5	0	0.18	0.09	(-0.38,0.19,0)	0.94
	PCA	0.2 0.2 0.5	0	0.72	0.60	(0,0,0)	0.93
	PCA varimax	0.2 0.2 0.5	0	0.31	0.25	(0,0,0)	0.93
	PCA quartimin	0.2 0.2 0.5	0	0.28	0.22	(-0.22,0.11,0.01)	0.93
	SPCA	0.2 0.2 0.5	0	0.36	0.29	(-0.31,0.16,0.02)	0.93
	PCA	0.5 0.2 0.5	0	0.73	0.61	(0,0,0)	0.93
	PCA varimax	0.5 0.2 0.5	0	0.36	0.28	(0,0,0)	0.93
	PCA quartimin	0.5 0.2 0.5	0	0.33	0.26	(-0.22,0.11,0.01)	0.93
	SPCA	0.5 0.2 0.5	0	0.38	0.30	(-0.32,0.16,0.02)	0.93
	PCA	0.2	0	0	0.5	0.72	0.61	(0,0,0)	0.85
	PCA varimax	0.2	0	0	0.5	0.28	0.23	(0,0,0)	0.85
	PCA quartimin	0.2	0	0	0.5	0.24	0.18	(-0.30,0.15,0.01)	0.85
	SPCA	0.2	0	0	0.5	0.16	0.08	(-0.34,0.17,0.01)	0.86
	PCA	0.2 0.2	0	0.5	0.72	0.60	(0,0,0)	0.85
	PCA varimax	0.2 0.2	0	0.5	0.29	0.23	(0,0,0)	0.85
	PCA quartimin	0.2 0.2	0	0.5	0.26	0.20	(-0.21,0.11,0.01)	0.85
	SPCA	0.2 0.2	0	0.5	0.32	0.26	(-0.33,0.16,0.02)	0.85
	PCA	0.5 0.2	0	0.5	0.73	0.60	(0,0,0)	0.85
	PCA varimax	0.5 0.2	0	0.5	0.34	0.27	(0,0,0)	0.85
	PCA quartimin	0.5 0.2	0	0.5	0.31	0.24	(-0.22,0.11,0.01)	0.85
	SPCA	0.5 0.2	0	0.5	0.35	0.28	(-0.33,0.17,0.01)	0.84
	PCA	0.2	0	0.5 0.5	0.73	0.62	(0,0,0)	0.88
	PCA varimax	0.2	0	0.5 0.5	0.30	0.24	(0,0,0)	0.88
	PCA quartimin	0.2	0	0.5 0.5	0.25	0.19	(-0.30,0.15,0.01)	0.88
	SPCA	0.2	0	0.5 0.5	0.18	0.09	(-0.35,0.18,0)	0.88
	PCA	0.2 0.2 0.5 0.5	0.72	0.61	(0,0,0)	0.88
	PCA varimax	0.2 0.2 0.5 0.5	0.30	0.24	(0,0,0)	0.88
	PCA quartimin	0.2 0.2 0.5 0.5	0.27	0.22	(-0.21,0.11,0.01)	0.88
	SPCA	0.2 0.2 0.5 0.5	0.36	0.29	(-0.29,0.15,0.02)	0.88
	PCA	0.5 0.2 0.5 0.5	0.74	0.61	(0,0,0)	0.88
	PCA varimax	0.5 0.2 0.5 0.5	0.36	0.28	(0,0,0)	0.88
	PCA quartimin	0.5 0.2 0.5 0.5	0.33	0.26	(-0.21,0.11,0.01)	0.88
	SPCA	0.5 0.2 0.5 0.5	0.39	0.31	(-0.29,0.16,0.02)	0.87

Table 5 :

 5 Variance decomposition with the factors estimated by SPCA (international business cycles).

	Country	F1 (Europe)	F2 (Latin Am.) F3 (Northern Am.) F4 (Dev. Asia) Common.
		Unadj. Adj. Unadj.	Adj.	Unadj.	Adj.	Unadj.	Adj.
	Canada	49.9	9.3	12.4	0.0	61.1	20.5	4.6	0.6	71.1
	Mexico	26.8	5.5	32.4	13.9	14.9	0.2	4.9	0.1	41.8
	USA	35.1	3.7	0.7	7.8	72.1	42.8	4.0	0.0	81.6
	Argentina	4.5	0.7	24.4	16.8	8.9	2.3	3.1	0.1	26.8
	Bolivia	0.4	6.0	7.2	9.4	4.0	6.2	1.5	2.7	19.7
	Brazil	32.1	7.2	44.2	20.4	10.5	0.7	15.4	1.7	56.2
	Chile	2.3	4.6	6.9	1.7	20.9	17.9	5.5	2.2	27.6
	Colombia	15.5	0.1	49.2	27.7	14.9	0.9	23.0	7.5	58.9
	Costa Rica	19.0	1.2	15.4	1.4	71.6	49.5	5.2	0.0	73.7
	Dominican Republic	4.5	0.0	8.3	3.7	4.9	0.7	3.1	0.5	10.5
	Ecuador	3.9	0.1	39.5	35.9	1.3	0.8	3.4	0.2	41.2
	El Salvador	14.5	0.1	8.0	0.3	49.5	33.8	1.4	0.4	50.3
	Guatemala	27.6	0.6	35.3	12.2	38.8	11.0	8.0	0.1	55.3
	Honduras	6.5	6.1	10.1	1.3	49.8	43.1	5.9	1.1	56.6
	Jamaica	14.3	7.0	0.7	1.2	4.6	0.0	9.7	3.8	18.9
	Panama	3.6	0.3	39.4	40.8	0.0	4.4	0.2	1.1	45.9
	Paraguay	1.3	1.4	18.6	15.1	1.3	0.0	8.1	4.2	23.4
	Peru	0.3	2.7	43.5	51.4	0.8	0.1	0.1	2.2	52.5
	Trinidad and Tobago	8.9	3.1	24.4	19.1	2.0	0.9	0.0	3.7	30.0
	Uruguay	0.8	12.0	28.3	36.9	1.5	1.5	0.0	0.1	42.5
	Venezuela	12.9	0.9	33.0	20.2	6.8	0.0	6.2	0.3	34.9
	Austria	68.0	41.8	12.3	0.0	16.9	0.9	15.0	0.7	69.6
	Belgium	83.6	47.3	23.5	1.1	22.2	1.0	16.1	0.2	85.9
	Denmark	54.8	34.0	5.5	0.8	22.5	0.4	3.1	1.1	57.2
	Finland	55.6	43.7	8.7	0.0	10.5	1.9	3.7	0.8	58.4
	France	87.8	52.8	14.2	0.1	28.1	0.0	11.8	0.0	87.9
	Germany	65.4	30.2	12.4	0.0	27.5	0.3	16.6	1.2	67.0
	Greece	50.9	32.7	6.9	0.2	16.8	0.0	4.2	0.4	51.6
	Iceland	17.9	2.6	15.0	4.4	20.4	5.1	0.4	2.2	30.1
	Ireland	6.1	4.4	0.5	4.0	6.2	2.7	0.3	1.9	14.4
	Italy	75.8	45.3	16.0	0.1	17.0	1.8	21.2	2.0	79.6
	Luxembourg	25.4	15.2	0.0	5.5	13.2	1.3	3.8	0.1	31.5
	Netherlands	75.5	39.1	10.7	0.4	37.7	1.9	6.0	0.8	78.4
	Norway	26.8	6.0	14.8	2.3	19.3	2.0	6.4	0.1	32.1
	Portugal	61.8	38.3	6.4	1.0	17.8	0.2	14.1	1.0	63.9
	Spain	68.1	52.4	5.7	1.1	18.4	0.2	3.0	1.7	71.8
	Sweden	50.7	40.6	10.8	0.5	11.4	1.0	0.3	5.7	57.6
	Switzerland	50.1	34.5	13.1	0.9	11.6	1.0	3.0	1.2	52.8
	United Kingdom	26.8	7.3	0.0	9.9	36.5	17.0	3.6	0.0	50.5
	Cameroon	0.2	0.7	0.5	1.0	0.1	0.3	0.3	0.3	1.9
	Côte d'Ivoire	10.8	5.4	5.0	1.1	5.0	0.1	0.0	1.6	13.3
	Kenya	0.3	1.4	0.5	0.0	8.3	9.9	0.0	0.3	10.7
	Morocco	3.8	17.4	1.3	1.7	5.6	19.0	2.6	6.1	29.7
	Senegal	2.7	0.6	5.0	1.3	7.1	3.3	6.8	3.5	12.6
	South Africa	30.8	8.8	37.5	17.4	12.2	0.1	4.7	0.3	48.3
	Zimbabwe	3.7	0.2	0.0	1.9	22.7	21.7	0.2	0.2	26.1
	Bangladesh	6.9	6.6	1.0	0.0	0.2	1.7	1.8	0.2	8.8
	India	6.1	0.7	4.8	1.0	3.7	0.1	5.5	1.8	9.6
	Indonesia	0.9	9.6	1.5	0.6	0.6	0.6	54.5	67.4	71.9
	Pakistan	14.2	11.7	1.4	0.2	0.8	2.6	6.9	2.0	18.9
	Philippines	4.0	0.1	10.0	5.2	1.7	0.1	7.7	3.4	14.0
	Sri Lanka	1.2	0.6	14.4	15.9	1.1	7.0	0.8	0.0	21.8
	Hong Kong SAR	27.2	2.7	13.9	0.7	18.2	1.2	38.6	18.9	50.4
	Japan	66.2	41.2	11.4	0.0	13.7	1.9	20.7	2.8	70.8
	Korea	14.3	2.6	0.2	10.8	7.2	0.5	51.1	42.7	63.3
	Malaysia	8.8	0.3	4.7	0.0	8.0	0.9	67.0	56.7	68.0
	Singapore	16.0	1.0	9.3	0.4	5.4	0.2	60.1	43.7	62.0
	Thailand	15.0	2.7	4.1	0.1	3.2	0.7	48.1	35.2	50.8
	Australia	20.1	3.5	4.2	0.0	32.0	14.5	0.0	3.3	37.5
	New Zealand	4.4	2.9	0.3	0.1	1.9	0.1	0.1	0.2	4.8

  . Its composition and time coverage have been updated several times. The latest modification is the publicly available Federal Reserve Economic Data Monthly Database (FRED-MD), designed by[START_REF] Mccracken | FRED-MD: A monthly database for macroeconomic research[END_REF]. It is now the classic benchmark dataset for macroeconomic factor models. The monthly updates and revisions are taken care of by the data specialists at the St. Louis Fed 13 using mostly the time series from the Federal Reserve Economic Data (FRED), the St. Louis Fed's main publicly available economic database. Resulting vintage databases are all available on its website.14 The variables are classified in

	8 groups: 1) output and income, 2) labor market, 3) housing, 4) consumption, orders and inventories,
	5) money and credit, 6) interest rates, spreads, and exchange rates, 7) prices, and 8) stock market.
	We use the 2018:9 vintage of FRED-MD: 128 variables, covering the 1959:1-2018:7 time span. In
	order to get a balanced panel, the analysis is restricted to 1960:1-2018:4, with 123 variables (see
	online Appendix C.1 for further details). The series were stationarized using the transformations
	recommended by McCracken and Ng (

Table B2 :

 B2 Correlations between the estimated factors (international business cycles).

	Factors	PCA varimax vs. PCA quartimin PCA varimax vs. SPCA PCA quartimin vs. SPCA
	F1 (Europe)	1.00	0.97	0.99
	F2 (Latin America)	1.00	0.90	0.93
	F3 (Northern America)	0.96	0.85	0.95
	F4 (Developed Asia)	0.99	0.87	0.90

Table B3 :

 B3 Correlations between the estimated factors and the real GDP growth rates of some representative countries (international business cycles).

	Pairs					PCA varimax PCA quartimin SPCA
	F1 (Europe) vs. France			0.91	0.92	0.94
	F2 (Latin America) vs. Brazil		0.47	0.52	0.66
	F3 (Northern America) vs. USA		0.69	0.79	0.85
	F4 (Developed Asia) vs. Malaysia	0.65	0.67	0.82
		e r ic
	( E u r o p e )	( L a t in A m
	F 1	F 2			1
	1	0.13	0.34	0.16	0.8
					0.6
					0.4
	0.13	1	0.07	0.07	0.2
					0
	0.34	0.07	1	0.08	-0.2
					-0.4
					-0.6
	0.16	0.07	0.08	1	-0.8
					-1

Table B4 :

 B4 Variance decomposition with the factors estimated by PCA varimax (international business cycles).

	Country	F1 (Europe) F2 (Latin Am.) F3 (Northern Am.) F4 (Dev. Asia) Common.
	Canada	39.3	0.8	30.2	0.5	70.7
	Mexico	22.1	18.2	3.0	0.2	43.5
	USA	26.5	4.0	46.9	0.0	77.5
	Argentina	1.5	15.4	6.9	2.0	25.7
	Bolivia	2.6	15.8	12.5	3.0	33.9
	Brazil	26.5	21.9	0.0	8.3	56.7
	Chile	0.1	4.0	22.3	4.9	31.2
	Colombia	8.7	26.8	4.8	16.3	56.6
	Costa Rica	8.6	3.3	50.4	2.7	65.0
	Dominican Republic	2.8	7.0	4.8	0.0	14.5
	Ecuador	1.7	31.7	0.0	4.2	37.6
	El Salvador	8.1	0.7	37.4	0.0	46.2
	Guatemala	17.1	18.5	20.5	1.7	57.8
	Honduras	1.1	2.6	40.9	6.9	51.5
	Jamaica	16.7	1.4	0.0	5.5	23.6
	Panama	3.3	40.8	4.3	0.1	48.6
	Paraguay	0.1	18.3	0.4	8.8	27.6
	Peru	0.0	46.4	0.6	2.5	49.5
	Trinidad and Tobago	8.5	22.9	0.0	4.4	35.8
	Uruguay	3.8	32.7	3.9	0.7	41.0
	Venezuela	9.8	16.9	0.4	6.7	33.7
	Austria	63.4	0.5	0.6	3.5	67.9
	Belgium	77.4	3.5	1.0	1.4	83.3
	Denmark	53.4	0.4	4.0	0.5	58.3
	Finland	57.7	0.7	0.0	1.7	60.1
	France	83.1	0.7	2.0	0.7	86.5
	Germany	55.8	0.2	5.0	5.0	66.0
	Greece	53.5	0.1	0.4	0.0	53.9
	Iceland	13.2	9.8	15.1	3.6	41.7
	Ireland	6.7	3.2	8.2	9.9	27.9
	Italy	73.0	1.1	0.1	6.9	81.0
	Luxembourg	22.8	6.1	9.4	0.1	38.3
	Netherlands	68.1	0.0	9.3	0.0	77.4
	Norway	21.7	3.7	6.4	1.0	32.8
	Portugal	56.8	0.4	0.9	7.2	65.4
	Spain	69.2	0.0	0.9	1.1	71.2
	Sweden	55.0	1.3	0.2	9.4	66.0
	Switzerland	48.5	3.2	0.0	0.1	51.8
	United Kingdom	24.9	10.8	20.5	0.0	56.2
	Cameroon	0.1	2.9	0.0	3.7	6.7
	Côte d'Ivoire	11.4	1.4	0.7	1.2	14.7
	Kenya	0.2	1.0	19.1	1.1	21.3
	Morocco	7.6	2.5	15.4	5.4	31.0
	Senegal	0.8	3.5	5.7	6.2	16.1
	South Africa	27.2	25.9	1.1	0.6	54.9
	Zimbabwe	2.2	2.5	21.6	0.5	26.8
	Bangladesh	6.9	0.0	0.6	2.3	9.8
	India	3.8	1.3	0.2	9.9	15.2
	Indonesia	2.7	0.3	0.0	56.4	59.5
	Pakistan	17.7	0.0	4.0	3.6	25.3
	Philippines	2.0	10.5	1.3	1.1	15.0
	Sri Lanka	2.1	16.1	9.2	0.1	27.5
	Hong Kong SAR	20.0	3.3	5.3	22.0	50.6
	Japan	67.0	0.1	0.3	8.1	75.5
	Korea	12.1	11.1	2.8	35.0	60.9
	Malaysia	4.2	0.3	6.7	42.0	53.1
	Singapore	12.3	1.7	0.1	39.6	53.7
	Thailand	14.7	0.0	1.7	39.7	56.0
	Australia	17.4	0.1	16.2	2.4	36.1
	New Zealand	5.8	0.0	1.4	5.8	13.0

Table B5 :

 B5 Variance decomposition with the factors estimated by PCA quartimin (international business cycles).

	Country	F1 (Europe)	F2 (Latin Am.) F3 (Northern Am.) F4 (Dev. Asia) Common.
		Unadj. Adj. Unadj.	Adj.	Unadj.	Adj.	Unadj.	Adj.
	Canada	43.2	20.0	2.1	0.2	49.6	26.2	0.0	1.2	70.7
	Mexico	25.3	13.5	21.8	15.8	9.9	2.0	1.3	0.0	43.5
	USA	29.3	10.1	2.2	5.8	63.0	43.1	0.2	0.2	77.5
	Argentina	2.6	0.0	16.8	14.3	9.2	6.2	2.8	1.4	25.7
	Bolivia	1.7	6.8	14.6	16.0	9.1	12.9	3.1	3.4	33.9
	Brazil	30.4	18.3	26.8	19.1	3.3	0.1	12.4	6.2	56.7
	Chile	0.5	2.0	4.7	3.4	22.5	21.7	5.4	4.1	31.2
	Colombia	12.4	2.1	31.0	24.0	10.6	3.8	20.1	13.5	56.6
	Costa Rica	12.2	0.3	4.9	2.2	60.1	47.5	4.4	1.7	65.0
	Dominican Republic	3.7	0.6	7.8	6.3	7.1	4.2	0.1	0.0	14.5
	Ecuador	2.8	0.4	33.4	30.5	0.4	0.0	5.3	3.3	37.6
	El Salvador	10.4	1.0	1.4	0.3	44.7	35.2	0.1	0.1	46.2
	Guatemala	21.8	5.2	22.3	15.7	32.5	17.8	3.6	0.8	57.8
	Honduras	2.7	1.2	3.5	1.9	43.2	39.5	8.1	5.7	51.5
	Jamaica	16.9	14.1	0.6	2.0	1.2	0.1	7.6	4.8	23.6
	Panama	4.1	3.4	42.1	40.0	1.6	5.1	0.4	0.0	48.6
	Paraguay	0.4	0.3	19.4	17.6	0.9	0.3	9.5	7.8	27.6
	Peru	0.1	0.3	45.4	46.4	0.9	0.5	2.0	3.3	49.5
	Trinidad and Tobago	9.3	7.2	24.2	21.8	0.8	0.1	2.7	5.7	35.8
	Uruguay	2.2	8.2	31.5	33.1	2.4	4.0	0.5	0.4	41.0
	Venezuela	12.1	5.2	19.7	15.1	2.9	0.1	9.1	5.3	33.7
	Austria	65.6	51.0	1.9	0.0	9.6	0.1	7.6	2.1	67.9
	Belgium	80.5	61.9	6.8	1.9	12.7	0.2	4.9	0.5	83.3
	Denmark	53.7	42.7	0.0	1.1	15.4	2.5	0.0	1.1	58.3
	Finland	57.1	52.7	1.9	0.2	5.5	0.0	0.2	2.8	60.1
	France	85.5	66.6	2.5	0.1	16.0	0.8	3.6	0.2	86.5
	Germany	59.2	39.1	1.3	0.0	18.9	3.2	9.5	3.4	66.0
	Greece	53.8	45.9	0.7	0.0	7.2	0.0	0.8	0.1	53.9
	Iceland	15.3	5.6	11.5	8.4	23.0	13.3	1.9	5.1	41.7
	Ireland	6.2	4.5	2.7	3.7	10.9	7.5	8.1	10.6	27.9
	Italy	75.2	62.2	3.3	0.3	5.2	0.6	12.7	4.9	81.0
	Luxembourg	23.2	14.8	4.1	7.7	17.7	7.9	0.6	0.0	38.3
	Netherlands	70.4	49.5	0.3	0.4	27.6	6.6	0.9	0.2	77.4
	Norway	24.3	12.4	5.6	2.6	14.8	5.0	2.4	0.4	32.8
	Portugal	58.6	44.8	0.0	1.2	9.7	0.3	12.2	5.5	65.4
	Spain	68.6	60.5	0.2	0.3	10.5	0.2	0.0	2.0	71.2
	Sweden	53.7	51.4	2.6	0.6	6.3	0.0	4.9	11.9	66.0
	Switzerland	49.3	42.3	5.4	2.1	4.9	0.0	0.2	0.5	51.8
	United Kingdom	25.6	13.9	7.9	13.2	31.8	18.3	0.3	0.0	56.2
	Cameroon	0.0	0.0	2.6	3.0	0.0	0.0	3.6	3.9	6.7
	Côte d'Ivoire	11.7	9.1	2.0	1.0	3.1	0.4	0.5	1.7	14.7
	Kenya	0.0	2.5	1.0	0.9	16.5	19.2	1.0	1.3	21.3
	Morocco	6.2	14.8	2.7	2.5	9.0	16.7	4.1	5.9	31.0
	Senegal	1.5	0.0	4.2	3.0	7.2	5.3	7.0	5.4	16.1
	South Africa	30.1	19.6	29.8	23.1	7.0	0.5	0.0	1.4	54.9
	Zimbabwe	2.7	0.1	2.0	3.1	23.0	20.9	0.2	0.6	26.8
	Bangladesh	7.0	6.6	0.1	0.0	0.0	0.8	3.1	1.9	9.8
	India	4.8	1.9	2.0	0.9	1.2	0.1	11.4	9.1	15.2
	Indonesia	1.6	5.4	0.5	0.3	0.1	0.0	53.2	57.1	59.5
	Pakistan	17.4	19.2	0.1	0.1	0.4	5.1	5.3	3.0	25.3
	Philippines	2.9	0.5	11.6	9.8	2.8	1.1	1.7	0.7	15.0
	Sri Lanka	2.2	3.7	16.6	15.9	5.5	10.0	0.2	0.0	27.5
	Hong Kong SAR	24.0	9.2	5.6	2.2	13.6	4.0	27.3	19.1	50.6
	Japan	68.4	58.6	1.2	0.0	3.5	1.1	13.9	6.1	75.5
	Korea	13.4	6.2	7.8	13.2	7.1	2.1	38.5	33.9	60.9
	Malaysia	6.3	0.3	1.0	0.1	10.9	5.9	44.9	39.7	53.1
	Singapore	15.1	6.5	3.3	1.0	2.4	0.0	44.5	36.9	53.7
	Thailand	16.3	11.3	0.4	0.0	0.0	2.4	44.4	37.7	56.0
	Australia	18.8	9.1	0.3	0.0	25.0	14.4	1.0	3.3	36.1
	New Zealand	5.5	4.8	0.1	0.0	3.0	1.1	4.5	6.5	13.0

Table C6 :

 C6 Output and income.

	Mnemonic	Description	TC BP
	RPI	Real Personal Income	5	1
	W875RX1	Real Personal Income Excluding Transfer Receipts	5	1
	INDPRO	IP Index	5	1
	IPFPNSS	IP: Final Products and Nonindustrial Supplies	5	1
	IPFINAL	IP: Final Products (Market Group)	5	1
	IPCONGD	IP: Consumer Goods	5	1
	IPDCONGD IP: Durable Consumer Goods	5	1
	IPNCONGD IP: Nondurable Consumer Goods	5	1
	IPBUSEQ	IP: Business Equipment	5	1
	IPMAT	IP: Materials	5	1
	IPDMAT	IP: Durable Materials	5	1
	IPNMAT	IP: Nondurable Materials	5	1
	IPMANSICS IP: Manufacturing (Standard Industrial Classification)	5	1
	IPB51222s	IP: Residential Utilities	5	1
	IPFUELS	IP: Fuels	5	1
	CUMFNS	Capacity Utilization: Manufacturing	2	1

Table C7 :

 C7 Labor market.

Table C9 :

 C9 Consumption, orders, and inventories.

	Mnemonic	Description	TC BP
	DPCERA3M086SBEA Real Personal Consumption Expenditures	5	1
	CMRMTSPLx	Real Manufacturing and Trade Industries Sales	5	1
	RETAILx	Retail and Food Services Sales	5	1
	ACOGNO	New Orders for Consumer Goods	5	0
	AMDMNOx	New Orders for Durable Goods	5	1
	ANDENOx	New Orders for Nondefense Capital Good	5	0
	AMDMUOx	Unfilled Orders for Durable Goods	5	1
	BUSINVx	Total Business Inventories	5	1
	ISRATIOx	Total Business: Inventories to Sales Ratio	2	1
	UMCSENTx	Consumer Sentiment Index	2	0

Table C10 :

 C10 Money and credit.

	Mnemonic	Description	TC BP
	M1SL	M1 Money Stock	6	1
	M2SL	M2 Money Stock	6	1
	M2REAL	Real M2 Money Stock	5	1
	AMBSL	St. Louis Adjusted Monetary Base	6	1
	TOTRESNS	Total Reserves of Depository Institutions	6	1
	NONBORRES	Reserves Of Depository Institutions	7	1
	BUSLOANS	Commercial and Industrial Loans	6	1
	REALLN	Real Estate Loans at All Commercial Banks	6	1
	NONREVSL	Total Nonrevolving Credit	6	1
	CONSPI	Nonrevolving Consumer Credit to Personal Income	2	1
	MZMSL	MZM Money Stock	6	1
	DTCOLNVHFNM Consumer Motor Vehicle Loans Outstanding	6	1
	DTCTHFNM	Total Consumer Loans and Leases Outstanding	6	1
	INVEST	Securities in Bank Credit at All Commercial Bank	6	1

Table C11 :

 C11 Interest rates, spreads, and exchange rates.

	Mnemonic	Description	TC BP
	FEDFUNDS	Effective Federal Funds Rate	2	1
	CP3Mx	3-Month AA Financial Commercial Paper Rate	2	1
	TB3MS	3-Month Treasury Bill: Secondary Market Rate	2	1
	TB6MS	6-Month Treasury Bill: Secondary Market Rate	2	1
	GS1	1-Year Treasury Constant Maturity Rate	2	1
	GS5	5-Year Treasury Constant Maturity Rate	2	1
	GS10	10-Year Treasury Constant Maturity Rate	2	1
	AAA	Moody's Seasoned Aaa Corporate Bond Yield	2	1
	BAA	Moody's Seasoned Baa Corporate Bond Yield	2	1
	COMPAPFFx	CP3Mx -FEDFUNDS	1	1
	TB3SMFFM	TB3MS -FEDFUNDS	1	1
	TB6SMFFM	TB6MS -FEDFUNDS	1	1
	T1YFFM	GS1 -FEDFUNDS	1	1
	T5YFFM	GS5 -FEDFUNDS	1	1
	T10YFFM	GS10 -FEDFUNDS	1	1
	AAAFFM	AAA -FEDFUNDS	1	1
	BAAFFM	BAA -FEDFUNDS	1	1
	TWEXMMTH Trade Weighted U.S. Dollar Index: Major Currencies	5	0
	EXSZUSx	Switzerland / U.S. Foreign Exchange Rate	5	1
	EXJPUSx	Japan / U.S. Foreign Exchange Rate	5	1
	EXUSUKx	U.S. / U.K. Foreign Exchange Rate	5	1
	EXCAUSx	Canada / U.S. Foreign Exchange Rate	5	1

Table C12 :

 C12 Prices.

	Mnemonic	Description	TC BP
	WPSFD49207	PPI: Finished Goods	6	1
	WPSFD49502	PPI: Finished Consumer Goods	6	1
	WPSID61	PPI: Intermediate Materials	6	1
	WPSID62	PPI: Crude Materials	6	1
	OILPRICEx	Crude Oil, Spliced WTI and Cushing	6	1
	PPICMM	PPI: Metals and Metal Products	6	1
	CPIAUCSL	CPI: All Items	6	1
	CPIAPPSL	CPI: Apparel	6	1
	CPITRNSL	CPI: Transportation	6	1
	CPIMEDSL	CPI: Medical Care	6	1
	CUSR0000SAC	CPI: Commodities	6	1
	CUSR0000SAD	CPI: Durables	6	1
	CUSR0000SAS	CPI: Services	6	1
	CPIULFSL	CPI: All Items Less Food	6	1
	CUSR0000SA0L2	CPI: All Items Less Shelter	6	1
	CUSR0000SA0L5	CPI: All Items Less Medical Care	6	1
	PCEPI	Personal Consumption Expenditures: Chain Index	6	1
	DDURRG3M086SBEA Personal Consumption Expenditures: Durable Goods	6	1
	DNDGRG3M086SBEA Personal Consumption Expenditures: Nondurable Goods	6	1
	DSERRG3M086SBEA	Personal Consumption Expenditures: Services	6	1

Table C13 :

 C13 Stock market.

	Mnemonic	Description	TC BP
	S&P 500	S&P's Common Stock Price Index: Composite	5	1
	S&P: indust	S&P's Common Stock Price Index: Industrials	5	1
	S&P div yield S&P's Composite Common Stock: Dividend Yield	2	1
	S&P PE ratio S&P's Composite Common Stock: Price-Earnings Ratio	5	1
	VXOCLSx	VXO	1	0

Table C15 :

 C15 Correlations between the estimated factors (FRED-MD).

	Factors	PCA varimax vs. PCA quartimin PCA varimax vs. SPCA PCA quartimin vs. SPCA
	F1 (Output)	0.99	0.95	0.99
	F2 (Prices)	1.00	1.00	1.00
	F3 (Spreads)	0.99	0.99	1.00
	F4 (Interest rates)	0.99	0.98	1.00
	F5 (Housing)	0.99	0.97	1.00
	F6 (Labor)	0.93	0.87	0.99
	F7 (Stock market)	0.99	0.96	0.98
	F8 (Money & credit)	0.99	0.95	0.96

Table C16 :

 C16 Correlations between the estimated factors and some representative observed variables (FRED-MD).

	Pairs	PCA varimax PCA quartimin SPCA
	F1 (Output) vs. INDPRO	0.92	0.95	0.96
	F2 (Prices) vs. CPIAUCSL	0.89	0.89	0.89
	F3 (Spreads) vs. T10YFFM	0.93	0.94	0.95
	F4 (Interest rates) vs. GS1	0.93	0.95	0.96
	F5 (Housing) vs. HOUST	0.97	0.99	0.99
	F6 (Labor) vs. PAYEMS	0.72	0.87	0.91
	F7 (Stock market) vs. S&P 500	0.90	0.91	0.94
	F8 (Money & credit) vs. M2SL	0.51	0.51	0.61

The commonality of a variable refers to the part of its variance explained by the estimated factors.

This is crucial in fields like biostatistics (so as to focus on only a few genes), or finance (in order to limit the number of assets in the portfolio, reducing the transaction costs).

[START_REF] Zou | (Housing) vs. HOUST F6 (Labor) vs. PAYEMS F3 (Spreads) vs. T10YFFM F4 (Interest rates) vs. GS1 F1 (Output) vs. INDPRO F2 (Prices) vs[END_REF], the first term of the optimization problem is not premultiplied by 1/T . However, it is clear that they have normalized the data by a 1/ √ T term, since they define the data covariance matrix as X X. Here, we keep the standard notation and define the data covariance matrix as X X/T , hence the premultiplication by 1/T .

This can be proved using the following steps: first, for a given A, the derivation of the objective function w.r.t. B shows that it is minimum for B = A; then the minimization with respect to A, with the condition A A = Ir, gives

It is clear that, by construction, the objective function decreases at each iteration so that the algorithm does converge.

The introduction of this constraint has almost no impact on the results of our two empirical applications.

Except for the very specific case of purely autoregressive models with uncorrelated errors[START_REF] Bergmeir | A note on the validity of cross-validation for evaluating autoregressive time series prediction[END_REF].

The paper considers the case where there is no ridge penalty, i.e. when κ2 = 0.

Here we drop the index 0 which Doz, Giannone and Reichlin (2011) use for the true values of the parameters.

Σe is denoted Ψ0 in Doz,[START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF] and Σe is denoted Ψ.

This work was supported by the French National Research Agency with the EUR grant ANR-17-EURE-0001. Earlier versions of this paper circulated under the title "Another look into the factor model black box: factor interpretation and structural (in)stability."

A Variance decomposition method

Let R 2 A,B denote the R 2 of the multivariate regression of the N A variables in the T × N A matrix A onto the N B variables in the T × N B matrix B: R 2

A,B =

where P B = B(B B) -1 B is the projection matrix onto the space spanned by the columns of B. . . , z n ) , we know (see [START_REF] Rockafellar | Convex Analysis[END_REF], part V, section 23) that z belongs to ∂h(b), the subdifferential of h in b, if and only if:

Thus, the elements of ∂f j (b) can be written as:

where z satisfies equations (A5).

We then obtain that b

if and only if it exists a z j such that:

and :

and Z (k) = (z 1 , . . . , z r ), the result of proposition 1 follows.

Proposition 1 gives the general form of B (k) as a function of A (k) at step 1 of iteration k.

We now study the general form of A (k) , and state the following proposition:

Proposition 2:

Under standard assumptions A1-A3 and CR1-CR4 of Doz, [START_REF] Doz | A two-step estimator for large approximate dynamic factor models based on Kalman filtering[END_REF], if κ 1 = O 1 √ N , the following result holds: for any k, A (k) can be written as

where M k is a r × r matrix which satisfies:

The proof of proposition 2 will be done by mathematical induction, and relies on the following intermediary lemma:

For k = 1 : we have A (1) = P and the result is true with M 1 = I r .

If we suppose that the result is true for k, then the lemma applies.

Besides, we know that A (k+1) is obtained from step 2 of iteration k in the following way:

Further, the SVD of X X T B is obtained in the following way :

• V k is a r × r matrix whose columns are normalized eigenvectors associated to the r eigenvalues of

• ∆ 2 k is the diagonal r × r matrix whose diagonal terms are the associated eigenvalues

From the lemma (ii), we know that:

Given the fact that P P = I r and that D = O P (N ), we first get:

so that:

where V k is a r × r matrix and

Thus ∆ 2 k , the diagonal matrix which contains the eigenvalues of B k) , has terms which all diverge linearly with N 2 , like the terms of D 2 , and ∆ k has terms which all diverge linearly with N . We also get:

Now, we have seen that U k is defined by:

As we know that A (k+1) = U k V k , we have A (k+1) A (k+1) = I r by construction.

Using (A22), we can write A (k+1) as:

Using equation (A19), and the fact that ∆ -1 k = O P 1 N , we then obtain:

As A (k+1) A (k+1) = I r , we get:

As P P = I r we know that P = O P (1).

.

N by assumption, so that M k = O P (1).

We thus get: M k+1 = O P (1).

It then follows from equation (A23) that:

which completes the proof of Proposition 2.

We are now ready to prove our consistency theorem for the estimated factors.

Proof of the Theorem:

We suppose that the algorithm converges at iteration k and we build the estimated loadings and factors as mentioned in the description of the algorithm. More precisely, we use the following notations: Note. The interpretable factors are labeled with their economic interpretation, and the uninterpretable factors with a question mark. Since F and Λ are identified up to a column sign change, for each factor, we impose the largest loading in absolute value to be positive.