Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
102,608 characters · 12 sections · 71 citation commands
Cover It Up! Bipartite Graphs Uncover Identifiability in Sparse Factor Analysis
A popular technique of dimension reduction in multivariate analysis is principal component analysis (PCA) which relies on the single value decomposition (SVD) of the sample covariance matrix $S_y$ of realizations of an $m$-dimensional random variable $Y$, see e.g. and:int. More specifically, only the eigenvectors $U_1$ corresponding to the $r$ largest eigenvalues $D_1$ are kept, while the variation explained by the eigenvectors $U_2$ corresponding to the remaining eigenvalues $D_2$ is ignored, i.e.
where ${\Lambda}= U_1 D_1^{1/2}$ is a $m \times r$ matrix, typically with $r \ll m$. In its original form, PCA is purely a data reduction technique without much insight into the data generating process.
A statistical modelling framework derived from PCA is probabilistic PCA (tip-bis:pro) which adds random noise to account for the unexplained variance due to dropping the term $U_2 D_2 U_2^ \top$ in (ref). Assuming w.o.l.g. that $Y$ is centered, $Y \sim N(0,\Omega)$ is assumed to arise from a zero-mean Gaussian distribution with covariance matrix $\Omega= {\Lambda} {\Lambda} ^\top + \sigma^2 I_{m}$. In this model, the covariance between the components $Y_i$ and $Y_\ell$ of $Y$ is explained by the inner product of row $i$ and $\ell$ of ${\Lambda}$, while a single parameter, $\sigma^2$, is present to control the fraction of unexplained variance, $\sigma^2/\Omega_{ii}$, for all components of $Y$. Considering for illustration a random variable $Y$ which is not only centered, but also standardized (i.e. $\Omega_{ii}=1$), it becomes apparent that probabilistic PCA relies on the rather strict assumption that the fraction of unexplained variance is the same for all components of $Y$.
More flexibility in this regard is obtained by the multi-factor model introduced by thu:vec which found numerous applications in applied multivariate analysis and will be the focus of the present paper. The model introduces an idiosyncratic variance $\sigma^2 _i$ to account for the unexplained variance of each components $Y_i$ of $Y$ and decomposes the covariance matrix $\Omega$ as
where ${\Lambda}$ is the $m \times r$ factor loading matrix, $\Sigma_0=\text{diag}(\sigma_1^2,\ldots,\sigma_m^2)$ is diagonal and $r$ is the so-called factor dimension. As discussed by the comprehensive textbooks of gor:fac and and:int, the multi-factor model is interesting both from a mathematical and a statistical perspective.
Starting with the pioneering work of rei:ide and and-rub:sta, the mathematical analysis centers around the question of identifiability of the parameter ${\Lambda}$ and $ \Sigma_0$ for a given $\Omega$, both in situations where the assumed factor dimension is equal to or different from the true factor dimension $r$, see e.g. sha:ide and bek-ten:gen. Many mathematical conditions have been proposed to address the various types of unidentifiability inherent in any factor model, see Fruehwirth2023When for a recent review.
One such condition is variance identification which ensures that the decomposition of $\Omega$ in ((ref)) is unique in the following sense. For any pair $({\beta}, \Sigma_{r})$, where ${\beta}$ is an $m \times r$ factor loading matrix and $\Sigma_{r}$ is a diagonal matrix such that $\Omega={\beta} {\beta} ^\top + \Sigma_{r}$, it follows that $\Sigma_{r}=\Sigma_0$ and hence the cross-covariance matrices ${\beta} {\beta} ^\top = {\Lambda} {\Lambda} ^\top$ are identical. In this case, the underlying loading matrix can be identified up to rotational invariance (and-rub:sta), i.e. ${\beta} = {\Lambda} P$ where $P$ is a permutation matrix. and-rub:sta provide following sufficient condition for variance identification, also known as row deletion property: after deleting any row from ${\Lambda}$, the remaining matrix contains two disjoint submatrices of rank $r$. In the present paper, we contribute to the mathematical aspects of multi-factor analysis by proving a sufficient condition for the row-deletion property based on the zero-nonzero pattern of the factor loading matrix and show how it can be verified in practice through an efficient algorithm.
Variance identification becomes vital when the factor dimension is unknown. A typical example where variance identification is violated are so-called spurious factors, where only a single non-zero factor loading is present in the corresponding column of ${\Lambda}$. Such spurious factors emerge in particular, when a factor model of dimension $k> r$ is employed to explain the covariance matrix $\Omega$ emerging from model (ref). rei:ide shows that if there is a solution to the decomposition (ref) with factor dimension $r$, then there exist infinitely many solutions with a larger factor dimension $k >r$. tum-sat:ide give a representation of these solutions, characterized by the $m \times k$ loading matrix ${\beta}$ and the diagonal matrix $\Sigma_{k}$. They show that the true loading matrix ${\Lambda}$ is embedded within ${\beta}$, but disguised by spurious factors and a rotation $P$. For $k= r+1$, for instance,
where $ M $ is a spurious column with a single non-zero factor loading. Obviously, the pair $({\beta}, \Sigma_{k})$ implies the same covariance matrix $\Omega$ as the true model and yields the same predictive distribution for $Y$ as the pair $({\Lambda}, \Sigma_0)$. On the other hand, model $({\beta}, \Sigma_{k})$ features a factor loading matrix of dimension $k> r$ and overestimates (inflates) the true factor dimension. However, a quick check immediately reveals that the $m \times k$ loading matrix ${\beta}$ violates variance identification, even if ${\Lambda}$ satisfies the conditions for variance identification: once the row corresponding to the single non-zero element in $M$ is removed, the remaining matrix has rank $r $ and does not contain two distinct sub-matrices of rank $k > r$. Consequently, any pair $({\beta}, \Sigma_{k})$ that violates variance identification should not be considered a reliable representation for recovering the number of factors, as it very likely overestimates the factor dimension.
This example clearly indicates that variance identification is also relevant for statistical factor analysis, in particular when the factor dimension is unknown. Statistical analysis for factor models centers around fitting a suitable model to realizations $y_1, \ldots,y_T$ of $Y$, typically using ML estimation (and:bay,liu-rub:max,rub-tha:em_alg) or a Bayesian inference. The Bayesian approach combines the likelihood derived from the factor model (ref) with a prior distribution on ${\Lambda}$ and $\Sigma_0$ and offers several attractive features. The use of proper priors on the idiosyncratic variances $\sigma^2 _i$, for instance, avoids Heywood problems common in ML estimation, where some of the estimated $\sigma^2 _i$ are negative, see e.g. Fruehwirth2024Sparse.
ML and Bayesian approaches differ fundamentally when the factor dimension $r$ is unknown. To choose $r$, ML estimation employs an incremental approach where a factor model with increasing factor dimension is refitted to the data and BIC-type model selection criteria are applied to estimate $r$, see e.g. bai-ng:det2002. In recent years, sparse Bayesian factor analysis became extremely popular in dealing with uncertainty regarding the factor dimension, see among many others luc-etal:spa,wes:bay_fac,gha-etal:bay,bha-dun:spa,con-etal:bay_exp,roc-geo:fas,kau-sch:bay,zha-etal:bay_gro,Fruehwirth2024Sparse. Sparse Bayesian factor analysis recovers the number of factors $r$ from the data in a one-sweep algorithm. It combines an overfitting factor model where the factor dimension is potentially bigger than $r$ with a prior on the factor loading matrix that introduces prior column sparsity in ${\Lambda}$. This allows us to learn the number of factors on the fly and the resulting posterior distribution $p(r\mid y_1, \ldots,y_T)$ allows uncertainty quantification with respect to $r$.
To illustrate one major difference between sparse Bayesian factor analysis and PCA, Figure (ref) compares the posterior distribution $p(r\mid y_1, \ldots,y_T)$ of the number factors for the first data set considered in Section (ref), with the scree plot obtained from PCA. The data are 52 weekly returns of 17 currencies of big trading partners of the Eurozone between January and December 2005. As common for financial data, a strong market factor is present in the scree plot (shown on the right hand side) in addition to several weaker factors which explain a small fraction of variance. The posterior distribution (shown on the left hand side) translates the ambiguity in the scree plot into a posterior distribution $p(r\mid y_1, \ldots,y_T)$ that puts considerable mass on the presence of 4 or 5 factors.
The rest of the paper is organized as follows. Our mathematical results are summarized in Section (ref) and applied to sparse Bayesian factor analysis in Section (ref). While our motivation comes from sparse Bayesian factor analysis, our mathematical insights are of a purely structural nature and potentially useful beyond this specific application. After a brief review of variance identification in Section (ref), we provide a counting rule on the zero-nonzero pattern of the factor loading matrix ${\Lambda}$, summarized in the binary matrix $\delta$, to check variance identification. This counting rule was introduced by sat:stu and only shown to be a necessary condition which has to be checked for ${\Lambda}$ and all possible rotations ${\beta} ={\Lambda} P$. Recently, Fruehwirth2023When were able to prove that this counting rule is sufficient for variance identification, provided that ${\Lambda}$ exhibits a so-called generalized lower triangular (GLT) structure. However, their proof heavily relies on assuming a GLT structure and is not easily extended to alternative structures or unconstrained loading matrices. As a first major contribution, we prove in Theorem (ref) in Section (ref) that this counting rule is sufficient for the row deletion property of and-rub:sta except for a set of Lebesque measure zero. As opposed to Fruehwirth2023When, Theorem (ref) does not require any structural constraints and can be applied to constrained and unconstrained loading matrices alike. Our proof relies on matching the binary zero-nonzero pattern $\delta$ of the factor loading matrix to a bipartite graph which is a mathematical object that captures the structure of $\delta$ but is invariant to permutations of the rows and columns of $\delta$, similar to the counting rule being invariant to permutations of the rows and columns of $\delta$.
Mathematically, the counting rule is a condition on all non-empty submatrices of $\delta$ with $q \in \{ 1, \ldots, r\}$ columns requiring that this submatrix has at least $2q+1$ non-zero rows. This condition could be checked for all $2^r -1 $ submatrices which becomes infeasible for increasing $r$. As a second major contribution of this paper, we design in Section (ref) an efficient algorithm for checking the counting rule and prove in Theorem (ref) that this condition can be verified in polynomial time. This algorithm is available as open source code\footnote{The source code is available at \url{https://hdarjus.github.io/sparvaride/}.} and can be applied regardless whether the loading matrix is constrained as in Fruehwirth2024Sparse or unconstrained as in kau-sch:bay.
In Section (ref), we return to sparse Bayesian factor analysis. Posterior inference in sparse Bayesian factor analysis is typically performed using simulation techniques such as Markov chain Monte Carlo, see Fruehwirth2024Sparse among many others. During sampling from the joint posterior distribution of all unknowns, identifiability conditions are either partially or completely ignored. Estimates of the quantities of interest are obtained by post-processing the posterior draws from such an unidentified model. In Section (ref), we investigate specifically the impact of conditions ensuring variance identification on recovering both the covariance matrix $\Omega$ (which is essential for prediction) as well as the true number of factors during such a post-processing step. We illustrate for simulated data that using only those posterior draws during post-processing for which a sufficient condition for variance identification based on the counting rule of sat:stu is fulfilled is instrumental in recovering the unknown factor dimension. On the other hand, recovering $\Omega$ and, hence, prediction is fairly robust to ignoring this condition. In addition, we estimate the number of factors in a financial application dynamically using a moving window approach over 100 overlapping periods. Whereas prediction is again robust to the presence of posterior draws that are not variance identified, including this condition avoids overfitting solutions that inflate the number of factors and leads to a sparser number of factors explaining the observed variation of the data. We conclude the paper with discussions included in Section (ref).
Let ${y}=(y_1,\ldots,y_T)$ be a sequence of $m$-dimensional observations, which are {centered around zero and} assumed to arise from a latent linear factor model with $r$ factors,
where $f_t$ is the $r$-vector of latent factors, ${\beta}$ is the $m\timesr$-dimensional matrix of factor loadings $\beta_{ij}$ with full column rank, and $\epsilon_t$ is the $m$-vector of idiosyncratic errors, for $t=1,\ldots,T$. In the basic factor model, the idiosyncratic errors are assumed to be iid $m$-variate Gaussian {random variable} $\epsilon_t\sim N_m(0,\Sigma_0)$, where $\Sigma_0=\text{diag}(\sigma_1^2,\ldots,\sigma_m^2)$ is diagonal. Furthermore, the latent factors are iid $r$-variate Gaussian {random variable} $f_t\sim N_r(0,I_r)$, where $I_r$ is the $r\timesr$-dimensional identity matrix, and the factors are independent from the idiosyncratic errors $\epsilon_s$ for all $s=1,\ldots,T$. When the latent factors are integrated out, this specification gives rise to the matrix decomposition of the covariance matrix $\mathbb{V}(y_t)=\Omega={\beta}{\beta}^\top+\Sigma_0$.
It is well known that any unitary matrix $G$ can be used to rotate the factor loadings ${\beta}$ into ${\beta} G$ without changing the covariance matrix $\Omega$. In this case, model (ref) changes to the observationally equivalent $y_t=({\beta} G)(G^{-1} f_t)+\epsilon_t$, and therefore ${\beta}$ is not uniquely identified, which is called rotational invariance. Further restrictions are {required} to achieve unique identification of ${\beta}$, {however we work with unconstrained loading matrices} in this paper.
{As mentioned previously, we} contribute to the literature on the identification of $\Sigma_0$, called variance identification. More precisely, we consider the basic factor model as part of a sparse Bayesian factor analysis (BFA) model, where the factor loading matrix ${\beta}$ follows an unknown zero-nonzero pattern that is estimated from the data along with the other parameters. {Sparsity is allowed a priori, but not enforced, because the estimated pattern may be a fully nonzero matrix}. In such a setting, the identifiability of $\Sigma_0$ is not known a priori, because it depends on the estimated zero-nonzero pattern of ${\beta}$. Further details on the sparse BFA model are given in Section (ref).
We emphasize, that the sufficient condition that we employ to achieve variance identification is a condition on ${\beta}$. In other words, we investigate the properties of ${\beta}$ to say something about the identifiability of $\Sigma_0$. and-rub:sta provide such a condition, later modified in tum-sat:ide for overfitting factor models, called the extended row deletion property in Fruehwirth2023When.
and-rub:sta show that $\text{RD}({r},{1})$ is a sufficient condition for the identification of $\Sigma_0$. Their result holds for every ${\beta}$ with real-valued entries and thus also for a matrix with some exact zero entries; this makes it relevant for sparse BFA. A second application of $\text{RD}({r},{s})$ is variance identification in overfitting factor models, when series-specific {(spurious)} factors are allowed, which do not contribute to the off-diagonal elements of $\Omega$, as introduced by tum-sat:ide. In both and-rub:sta and tum-sat:ide, however, the theory lacks practical ways to verify $\text{RD}({r},{s})$. Later, sat:stu develop a new condition for ${\beta}$, based on counting nonzero rows in all rotations of ${\beta}$, which is shown to be necessary for $\text{RD}({r},{s})$. However, since the authors consider infinitely many rotations of ${\beta}$, {this} condition is unverifiable {in practice} and is not widely applied in the literature to our knowledge.
Recently, Fruehwirth2023When directly build on results by and-rub:sta, tum-sat:ide, and sat:stu, and introduce a framework in sparse BFA for the joint identification of ${\beta}$ and $\Sigma_0$. The framework is based on an identifying assumption\footnote{ ${\beta}$ is assumed to adhere to the so-called generalized lower triangular structure. For details, refer to Fruehwirth2023When.} on ${\beta}$ combined with the following counting rule, which is in the same spirit as the counting rule of sat:stu.
The authors show, among others, that, in their framework, an $m\timesr$ binary matrix ${\delta}$ satisfying $\text{CR}({r},{s})$ is sufficient for almost all $m\timesr$ factor loading matrices ${\beta}$ to satisfy $\text{RD}({r},{s})$, where ${\delta}$ is an indicator matrix for ${\beta}\neq0$; that is, $\beta_{ij}=0$ if ${\delta}_{ij}=0$. We generalize this result {to unconstrained loading matrices} in the next section.
We fix the notation and the terminology for the rest of the section. Denote by ${\beta}$ an $m\timesr$-dimensional matrix with real-valued elements $\beta_{ij}$, where we pay attention to zeros for their special interpretation in sparse Bayesian factor analysis. Additionally, denote by ${\delta}$ an $m\timesr$-dimensional binary matrix of zeros and ones. We say that ${\beta}$ is generated by ${\delta}$ if $\beta_{ij}=0$ whenever ${\delta}_{ij}=0$, and we denote the set of all ${\beta}$ generated by ${\delta}$ by $\mathcal{B}$. We think of $\mathcal{B}$ as being equivalent to a continuous probability space on {$\mathbb{R}^d$, where $d={\sum_{i,j}{\delta}_{ij}}$ is equal to the total number of non-zero elements in ${\beta}$. This definition} is motivated by the slab elements of spike-and-slab priors on ${\beta}$ in Bayesian variable selection (tadesse_handbook_2021). The probability space brings us to an expression that we extensively use below: “for all ${\beta}$ generated by ${\delta}$, except for a set of Lebesgue-measure zero” is formally understood as “the set of exceptional ${\beta}$ matrices form a Lebesgue-nullset in $\mathcal{B}$”.
This section is mainly concerned with {proving} the sufficiency of $\text{CR}({r},{s})$ for $\text{RD}({r},{s})$, formally stated in Theorem (ref) {at the end of this section. The} proof is presented in multiple pieces through Propositions (ref), (ref), and (ref), and surrounding lemmas. Since both $\text{CR}({r},{s})$ and $\text{RD}({r},{s})$ imply $m\ge2r+s$ for $m\timesr$ matrices, we assume that $m\ge2r+s$ throughout this section. First, {we show in Proposition (ref) that the problem can be} reduced to the case of $s=0$. Later, {based on Proposition (ref),} we only need to prove that $\text{CR}({r},{0})$ implies $\text{RD}({r},{0})$, except for a set of Lebesgue-measure zero.
{Before we exploit the simplification provided by Proposition (ref) below in Propositions (ref) and (ref)}, we take a {detour} into elementary graph theory to prove our results in a structured way. The centerpiece of the {final} proof of Theorem (ref) is the classical duality theorem by K\H{o}nig kon:gra and Egerv\'ary ege:mat in graph theory, which is stated later in this section, where we put $\text{CR}({r},{0})$ and $\text{RD}({r},{0})$ on the two sides of the duality. Graph theory also provides us with a convenient representation of the problem through a specific mapping of a binary matrix ${\delta}$ to its corresponding bipartite graph, {which will be} introduced in Definition (ref). Notably, these bipartite graphs are an equivalent representation of ${\delta}$ up to reordering of its rows and columns. This is a helpful framework, since both $\text{CR}({r},{0})$ and $\text{RD}({r},{0})$ are invariant to row and column permutations. In the following, we think of a nonzero entry ${\delta}_{ij}=1$ in ${\delta}$ as a line connecting a point $z_i$ that represents row $i$ with another point $c_j$ that represents column $j$. This construction is formally defined next in order to provide the language for our proofs. For a more detailed introduction to graph theoretic notions, see chapters 3.1 and 3.2 of wes:gra.
Figure (ref) shows a {$3 \times 3$} binary matrix ${\delta}$ and its bipartite graph $B=(V_\text{row},V_\text{col},E_B)$ with $V_\text{row}=\{z_1,z_2,z_3\}$ and $V_\text{col}=\{c_1,c_2,c_3\}$. The edges are $E_B=\{\{z_1,c_1\},\{z_1,c_2\},\{z_1,c_3\},\{z_2,c_1\},\{z_3,c_1\}\}$. One way we use graph theory is {displayed in the same figure:} thick edges correspond to a so-called maximal matching in $B$, formally introduced below in Definition (ref). {Loosely speaking, we are trying to reorganize the rows of ${\delta}$ such that the main diagonal exhibits as many ones as possible. } {For the ${\delta}$ considered in Figure (ref),} swapping rows $z_1$ and $z_2$ results in a reorganized ${\delta}$ with a diagonal of two ones, shown in bold face in the original ${\delta}$, which is the {maximum number that can be achieved in this case}. {However, if element ${\delta}_{3,3}$ were also a one, then ${\delta}$ could be reorganized such that the diagonal was full and contained ones, only. In this case,} the maximum matching in $B$ would be of size three.
Let us zoom out of this particular example: the existence of a {submatrix with a full diagonal of ones} will be the proof of full-rankedness of all factor loading matrices ${\beta}$ generated by ${\delta}$, except for a set of Lebesgue-measure zero. {This is stated in Lemma (ref).} {To generalize Lemma (ref), we will then introduce the notion of a matching in bipartite graphs. It will be shown that graph matching allows us to simplify the search for a submatrix ${\delta}^{\scalebox{0.5}{$\,\square$}}$ that satisfies the conditions of Lemma (ref). Instead of rearranging the rows of ${\delta}$ till a suitable ${\delta}^{\scalebox{0.5}{$\,\square$}}$ is found, a dual optimization task is performed on the bipartite graph $B$ corresponding to $\delta$.}
A matching in a bipartite graph implicitly encodes both a subset of rows and their permutation in the bi-adjacency matrix ${\delta}$ {in such a way} that the permuted rows form a submatrix in ${\delta}$ with a full diagonal of ones. The key is the unique correspondence between row indices and column indices that take part in the matching. A small example {has been provided} in Figure (ref), and larger examples are the thick edges in the two graphs in Figure (ref). In the latter, for instance, by reordering vertices $z_1,\ldots,z_7$ such that thick edges do not cross, we obtain a full diagonal of ones in the reordered ${\delta}$.
The following statement generalizes Lemma (ref) {where ${\delta}^{\scalebox{0.5}{$\,\square$}}$ was supposed to exhibit a full diagonal of ones.} In Lemma (ref), we implicitly find ${\delta}^{\scalebox{0.5}{$\,\square$}}$ even if the {rows of $\delta$ are not ordered such that an approriate submatrix with a full diagonal of ones exists. We show that we can find ${\delta}^{\scalebox{0.5}{$\,\square$}}$ with the help of a $V_\text{col}$-saturating matching in the bipartite graph $B$ of $\delta$}.
In the proof, we use a correspondence between a full-column-rank submatrix in ${\beta}$, generated by ${\delta}$, and a $V_\text{col}$-saturating matching in the bipartite graph $B$ of ${\delta}$. {For illustration, consider the matrix} ${\delta}$ and its bipartite graph $B=(V_\text{row},V_\text{col},E_B)$ in {the top of} Figure (ref), where $V_\text{row}=\{z_1,\ldots,z_7\}$, $V_\text{col}=\{c_1,c_2,c_3\}$, and $E_B$ is the set of edges connecting $V_\text{col}$ and $V_\text{row}$. The $V_\text{col}$-saturating matching $\{\{z_1,c_2\},\{z_2,c_3\},\{z_4,c_1\}\}$ is shown using thick edges in $B$, and the same positions in ${\delta}$ are typed in bold face. These elements form a full diagonal of three ones in ${\delta}$ after reordering its rows beginning with $z_4$, $z_1$, and $z_2$.
{So far, we} have discussed a condition that ensures the almost sure existence of a full-column-rank submatrix in ${\beta}$ generated by ${\delta}$. To satisfy $\text{RD}({r},{0})$, however, we look for two {distinct} submatrices of ${\beta}$ that are of rank $r$, i.e. two disjoint sets of $r$ rows $({\beta}_{i_1,\cdot}, \ldots, {\beta}_{i_{r},\cdot})$ and $({\beta}_{i_{r+1},\cdot}, \ldots, {\beta}_{i_{2r},\cdot})$ in ${\beta}$ such that the corresponding submatrices $({\beta}_{i_1,\cdot}^\top \ldots {\beta}_{i_{r},\cdot}^\top)^\top$ and $({\beta}_{i_{r+1},\cdot}^\top \ldots {\beta}_{i_{2r},\cdot}^\top)^\top$ are both of rank $r$. According to Lemma (ref), assuming for now that the rows of ${\beta}$ are ordered {appropriately}, this amounts to finding two disjoint $r\timesr$-dimensional submatrices of ${\delta}$ with a diagonal of ones, up to a Lebesgue-nullset. In Lemma (ref), we replace this search task with a simpler one {by introducing the notion of {\em duplicated binary (DB) matrices}.} An example {of a binary matrix and the corresponding DB matrix} is {provided in} Figure (ref) for $m=7$ and $r=3$.
{To utilize Lemma (ref) for our final goal, we resort to graph theory to verify the existence of a $(2r)\times(2r)$-dimensional submatrix with a diagonal ones in the DB matrix ${\delta}^{||}$ of $\delta$. To this aim we introduce the notation of a duplicated bipartite (DB) graph of a binary matrix $\delta$ in Definition (ref).} {For illustration, the} DB matrix ${\delta}^{||}$ and {the corresponding} DB graph, denoted by $B^{||}$, are shown in Figure (ref) for a $7\times3$-dimensional binary matrix ${\delta}$. {Lemma (ref) then relates the existence of a $V_\text{col}^{||}$-saturating matching in $B^{||}$ to the existence of two “disjoint” $V_\text{col}$-saturating matchings in $B$.}
An instructive way to think of Lemma (ref) is it being the same as Lemma (ref), but with a permutation $\rho$ applied to the rows of ${\delta}$. Since ${\delta}^{||}$ has the same row labels as ${\delta}$, $\rho$ can be applied to the rows of ${\delta}^{||}$ as well. This way, $\rho$ maps task (i) in Lemma (ref) to task (i) in Lemma (ref) and does the same for task (ii). The inverse $\rho^{-1}$ maps the tasks back from Lemma (ref) to Lemma (ref), thus closing the loop. Here, we present a direct mechanical proof.
One final statement considers the size of a maximum matching and thus concludes one side of the aforementioned duality theorem by K\H{o}ning and Egerv\'ary, which is formally introduced later with its necessary terminology.
Now, we place $\text{CR}({r},{0})$ onto the other side of the duality. $\text{CR}({r},{0})$ is a statement about columns of ${\delta}$ being connected to sufficient number of its rows via ones. Below, that notion of sufficiency is translated into the language of bipartite graphs. We show that $\text{CR}({r},{0})$ is sufficient for $B^{||}$ to have many edges in a specific sense. So many that one cannot do better than pick the entire $V_\text{col}^{||}$ as a set of vertices to cover all edges $E_{B^{||}}$. In the following, we define vertex covers for bipartite graphs and present the duality theorem.
{For illustration, consider Figure (ref).} On the right hand side, $\{c_1,c_2,c_3\}$ is a vertex cover in $B$ because all edges in $E_B$ touch at least one of these vertices. There is a smaller vertex cover: the set $\{z_1,c_1\}$. This vertex cover has size two, which is also the size of a matching in $B$, shown with thick edges. According to the duality theorem by K\H{o}nig and Egerv\'ary, said matching is therefore a maximum matching and $\{z_1,c_1\}$ a minimum vertex cover. With the help of the duality theorem and Proposition (ref), it only remains to show that {the DB graph $B^{||}$ of $\delta$} has a minimum vertex cover of size $2r$ if ${\delta}$ satisfies $\text{CR}({r},{0})$. In {this case,} the duality theorem implies the existence of a matching of size $2r$ in $B^{||}$ and, {consequently,} Proposition (ref) implies that all factor loading matrices ${\beta}$ generated by ${\delta}$ satisfy $\text{RD}({r},{0})$, except for a set of Lebesgue measure zero.
{The following Lemma (ref) resembles a counting rule $\text{CR}({r},{0})$ for DB matrices and is instrumental in characterizing vertex covers in $B^{||}$. Lemma (ref) is employed to prove the Proposition (ref), which is the final piece required for the proof of Theorem (ref).}
We conclude the section with a corollary that is applied in Section (ref) to identify variance identified models in sparse Bayesian factor analysis.
The link from ${\beta}$ to variance identification is the Anderson-Rubin theorem and-rub:sta stating that $\Sigma_0$ is identified if ${\beta}$ satisfies $\text{RD}({r},{1})$. Note, however, that the Anderson-Rubin theorem and thus $\text{CR}({r},{1})$ are not necessary for variance identification even in sparse BFA. Appendix (ref) provides an example of a sparse variance identified model that does not satisfy $\text{RD}({r},{1})$.
In order to apply Corollary (ref) in practice, one needs to verify $\text{CR}({r},{1})$ for a given binary matrix ${\delta}$. We establish the applicability of Corollary (ref) for large Bayesian sparse factor models in the next section.
We extend the previous section and describe an efficient algorithm that verifies $\text{CR}({r},{1})$. An initial idea might be to visit all the nonempty submatrices of ${\delta}$ that consist of $q$ columns for $1\le q\ler$ and count the number of nonzero rows. However, that approach examines $2^r-1$ matrices, which is computationally infeasible for large $r$: in {Bayesian inference,} where ${\delta}$ is sampled from {the posterior distribution and many binary matrices need} to be checked, this step may {induce considerable} computational cost. In this section, we develop a representation of the verification task in graph theory that helps us to prove our second main result: a feasible algorithm for the verification of $\text{CR}({r},{1})$ {even for large $m$}. Formally, we show in Theorem (ref) that $\text{CR}({r},{1})$ can be verified in a number of steps that is polynomial in $m$ and $r$. In this section, we provide a constructive proof of the theorem, which can be implemented in practice to verify variance identification in sparse BFA based on Corollary (ref).
{At this point, a few remarks concerning zero rows and zero columns in ${\delta}$ are in order. Obviously, } zero columns are not allowed in $\text{CR}({r},{1})$ matrices. {On the other hand, zero rows might be present and} can be removed from ${\delta}$ without loss of generality. In particular, the addition or removal of zero rows does not influence rank conditions and $\text{RD}({r},{1})$ holds for ${\beta}$ if and only if it holds for ${\beta}$ without its zero rows. Henceforth, we assume that every row and column of ${\delta}$ has at least one nonzero element.
Now, we introduce an extended notion of bipartite graphs that allows us to represent the verification of $\text{CR}({r},{1})$ in a graph-theoretical framework.
Now we present Propositions (ref) and (ref), which constitute the two pieces for the proof of Theorem (ref) In Proposition (ref), we design a weighted bipartite graph $B$ such that the verification of $\text{CR}({r},{1})$ on ${\delta}$ is equivalent to computing the total weight $M^\star$ of the MWVC on $B$. Finally, in Proposition (ref), we show that the MWVC at hand can be solved efficiently via a polynomial algorithm.
Then, the following proposition provides the basis for the polynomial algorithm. The intuition behind the vertex cover is that the submatrix formed by the rows and columns that are left out is a zero matrix in ${\delta}$.
In the next statement, we use the Big-$O$ notation to describe the computational complexity of the algorithm.
See the proof below. We do not directly work on the weighted bipartite graph to find the MWVC, but we rather reformulate the problem as a minimal network cut problem and refer to known solutions for that problem, such as Dinic's algorithm tar:dat. Therefore, in order to present the reformulation, we first introduce some notions from network theory. Even though a network is also a graph, we deliberately use different terminology for its parts to improve readability. In particular, we use “node” instead of “vertex” and “arrow” instead of “edge”. We denote arrows as tuples by round brackets, e.g., $(c_2,z_1)$, because they are directed and the order of the nodes matters in networks, in contrast to the set-notation of the curly brackets $\{z_1,c_2\}$ used for undirected edges in all graphs in this paper.
{For illustration, consider Figure (ref).} The {upper part} shows an example of the weighted bipartite graph $B$ for a $8\times3$-dimensional ${\delta}$. The bottom half of Figure (ref) shows an example of a network $N=(V,E,\kappa)$ created from $B$ as described below in the proof of Proposition (ref), where $V=\{s,c_1,c_2,c_3,z_1,z_2,z_3,z_4,\allowbreak z_5,z_6,z_7,z_8,t\}$. There are twenty-four edges in $E$, and capacities are $\kappa(s,c_i)=7$, $\kappa(c_i,z_j)=\infty$ if ${\delta}_{j,i}=1$ and otherwise $0$, and $\kappa(z_j,t)=3$. A cut can be, for example, $C=\{s,c_1,z_1,z_3,z_4,z_5\}$ with capacity $\kappa(C)=2\cdot7+0+4\cdot3=26$.
In the proof of Proposition (ref), we also find that the number of steps increases with $P(r,m)=(m+r+2)^2(m+r+mr)$. For fixed $r$, the computational complexity of our method is therefore $O(m^3)$, and, for fixed $m$, which is often the case, it is $O(r^3)$ instead of complexity $O(2^r)$ for the brute force search through all submatrices.
Although Theorem (ref) only concerns $\text{CR}({r},{1})$, the result may also be used to build an algorithm that verifies $\text{CR}({r},{s})$ in smaller settings. Namely, it is easy to see that ${\delta}$ satisfies $\text{CR}({r},{s})$ if and only if after removing any $s-1$ rows the remaining binary matrix satisfies $\text{CR}({r},{1})$, which we can verify efficiently. That realization gives rise to a recursive algorithm with complexity $O(m^{s-1} P(r,m))$, which may be practical for $s=2$ or $s=3$ for small $m$.
Finally, note that, in its current form, the proof cannot be extended to a polynomial complexity $P(r,m,s)$ algorithm also in $s$ for $\text{CR}({r},{s})$ by choosing different weights for $V_\text{row}$ or $V_\text{col}$. In particular, if $x$ denotes the ratio of vertex weights in $V_\text{col}$ and $V_\text{row}$ (i.e., $x=(2r+1)/r$ above), then it can be shown that $x>2+s-1$ and $x\le2+s/r$ are both necessary for the proof of Proposition (ref) and thus for Theorem (ref). This interval is non-empty only if $s\le1$.
{In concluding we note that implementations of this algorithm are available in R and MATLAB at \url{https://hdarjus.github.io/sparvaride/}}.
We demonstrate that missing variance identification may unnecessarily inflate the estimated number of factors during exploratory factor analysis (EFA). We choose the Bayesian paradigm, which allows us to emulate matrix sparsity using a spike-and-slab prior distribution on ${\beta}$ ({to be} introduced in Section (ref)) and to consider variance identification as a domain restriction on that prior distribution. Consequently, we can estimate the posterior distribution via a Markov chain Monte Carlo (MCMC) sampler under the unrestricted prior and apply the domain restriction as a post-processing step by discarding the unsatisfactory draws. The model, its estimation, a simulation study, and a real data study are detailed below.
To facilitate variance identification through $\text{CR}({r},{1})$, we follow the tradition of wes:bay_fac and introduce indicator variables ${\delta}_{ij}\in\{0,1\}$ for every factor loading $\beta_{ij}$ as parameters to estimate for $i=1,\ldots,m$, and $j=1,\ldots,{H}$, and collected in the $m\times{H}$ matrix ${\delta}=({\delta}_{ij})$. Following established procedures (con-etal:bay_exp,kau-sch:bay,Fruehwirth2024Sparse), Bayesian posterior sampling is applied with a conjugate prior on ${\beta}$ and $\Sigma_0$, combined with column-wise shrinkage on the indicator {variables ${\delta}_{ij}$}. For completeness, we {provide the full hierarchical model specification by combining model (ref) with a corresponding prior:}
where $IG(c_0, C_0)$ denotes the inverted gamma distribution with kernel density $x^{-c_0-1}\exp(-C_0/x)\mathbbm{1}(x>0)$, $Ber(\tau_j)$ is the Bernoulli distribution with success probability $\tau_j\in(0,1)$, $B(a_0, b_0)$ is the beta distribution with kernel density $x^{a_0-1}(1-x)^{b_0-1}\mathbbm{1}(0<x<1)$, and $t=1,\ldots,T$. The choice of $\sigma_i^2$ as the variance lets $\beta_{ij}$ capture potential scaling differences between the observation series. Moreover, two settings are considered below for the prior on $\tau_j$: following roc-geo:fas and Fruehwirth2023Cusp, the finite one-parameter beta (1PB) prior $(a_0,b_0)=(\alpha/{H},1)$ is chosen first, which we call shrinkage below, and the uniform prior $(a_0,b_0)=(1,1)$ is picked as an alternative for sensitivity analysis.
A potentially influential question is the choice of ${H}$. One solution is the use of infinite factor models, initiated by gha-etal:bay and popularized by bha-dun:spa and leg-etal:bay, where one theoretically lets ${H}$ diverge to $\infty$ while cumulatively shrinking the columns a priori {towards zero} as the column index increases. Here, we assume\footnote{ Necessarily, $2{H}+1\lem=m$, where $m$ is the number of observation series. That is essential for variance identification via $\text{RD}({r},{1})$, and therefore also via $\text{CR}({r},{1})$.} ${H}=\min(30,\lfloor(m-1)/2\rfloor)$ to both allow for parameter identification via Corollary (ref) and keep Monte Carlo simulations manageable. Notably, recently, Fruehwirth2023Cusp showed that our column-wise exchangeable prior $p({\delta}\mid{H})$ in Equation (ref) is strongly related to both the framework of gha-etal:bay and leg-etal:bay.
Model (ref) specifies a sparse Bayesian factor model with a spike-and-slab prior on ${\beta}$. The prior $p({\beta}\mid{H})$ is exchangeable both row-wise and column-wise, and the elements of $\{\sigma_i^2\}$ are independent a priori, which results in an order-independent model for the observation series. Furthermore, the choice of standard conjugate priors for $({\beta},\{\sigma^2_i\})$ and $\{\tau_j\}$ enables simple Gibbs sampling. See Appendix (ref) for the steps of the MCMC algorithm.
Throughout the demonstration, we compare three domain restrictions, which we implement via post-processing of the MCMC output. Under the unrestricted scenario, variance identification as a step is ignored, and the entire output of the MCMC procedure is retained. In the second scenario, the necessary condition for variance identification of and-rub:sta is applied as a post-processing step, similar to kau-sch:bay. Namely, if in all columns of ${\beta}$, at least three nonzero elements are present, then the MCMC draw is retained, and, otherwise, it is excluded from summaries of the posterior distribution. In the third scenario, the sufficient condition $\text{CR}({r},{1})$ from Corollary (ref) is enforced during post-processing by only keeping the MCMC draws that satisfy the condition. In both cases, the MCMC output is filtered before proceeding further: before any {subsequent} analysis, we discard the joint draws of {$({\beta},\Sigma_0,f_1, \ldots, f_T)$} when ${\beta}$ does not satisfy the necessary or, respectively, the sufficient condition.
Further steps during post-processing are estimating {the number of factors} $r$ and the covariance matrix $\Omega={\beta}{\beta}^\top+\Sigma_0$ from the filtered or unfiltered MCMC output, depending on the scenario above. Following Fruehwirth2023Cusp and Fruehwirth2023When, we assume a potentially too large number of factors ${H}$ and estimate the posterior distribution for {$r$} by counting the number of active columns in ${\beta}$ for every MCMC draw. Active columns of ${\beta}$ are those that contain at least two nonzero elements, and zero columns are deemed inactive. Columns with a single nonzero element are automatically transformed to zero columns during post-processing by moving the square of the single factor loading and adding {it} to the corresponding diagonal element of $\Sigma_0$. The reason is that these columns are actually {spurious factors and} they capture the variance of a single observation series, as explained in Fruehwirth2023When. Finally, one acquires a posterior sample for the covariance matrix by calculating $\Omega = {\beta}{\beta}^\top+\Sigma_0$ for every joint draw of $({\beta},\Sigma_0)$.
We follow leg-etal:bay and conduct a simulation study with three different combinations of $(m,r)$, namely, $(20,5)$, $(50,10)$, and $(100,15)$. For each combination, $25$ repetitions of {$T=100$} observations are generated. Following Fruehwirth2023Cusp, we examine two settings for generating ${\delta}$: in the dense setting, ${\delta}$ is a fully nonzero {binary} matrix, and in the sparse setting, random 30% of {the indicators in ${\delta}$ are} set to zero and the remaining 70% to one. We always enforce the true ${\delta}$ to satisfy $\text{CR}({r},{1})$ by re-sampling until {this condition is met}. In all scenarios, $\Sigma_0$ is the identity matrix, and $\beta_{ij}$ is standard normal {whenever} ${\delta}_{ij}$ is nonzero.
Turning to the priors, the fairly vague setting $(c_0, C_0)=(1, 0.3)$ is adopted from leg-etal:bay. Finally, contrary to Fruehwirth2023Cusp, we do not estimate $\alpha$ to keep the model simple but rather fix $\alpha=5$, which is consistent with their findings. Including the choice of shrinkage and uniform priors for $\tau_j$, 300 posterior distributions are estimated in this simulation study in total.
To facilitate MCMC convergence diagnostics, four independent posterior Markov chains are simulated with distant i\-ni\-tia\-li\-za\-tions: zero, one, ${H}-1$, and ${H}$ randomly filled columns in ${\beta}$ with standard normal draws. In the small settings $(20,5)$ and $(50,10)$, the MCMC chains are run for $50\,000$ iterations, and the first $10\,000$ are discarded as burn-in. However, we face significant computational challenges with our simple Gibbs sampler in the biggest setting $(100,15)$, where we run the MCMC chains for one million iterations on a cluster of 400 cores and one terrabyte memory for a total of 20 hours to see full convergence.
Figures (ref) and (ref) provide details on the results under the shrinkage prior on $\tau_j$ and follow a similar structure. The six facets of Figure (ref) depict the posterior probability {$p(r=r_\text{true} \mid{y})$} of the true number of factors, where ${y}$ are the observed data and {$r_\text{true}$} is the true number of factors in the data generating process (DGP). The first and the second rows correspond to the dense and, resp., the sparse setting, while the columns correspond to the true number of factors $r_\text{true}$. Within a facet, from left to right, the three boxplots summarize posterior probabilities under the unrestricted, the necessary, and, respectively, the sufficient scenario, each showing a distribution over 25 DGP repetitions. The final ingredients of the chart are the lines that connect the corresponding repetitions, i.e., posterior summaries under different scenarios but the same data set. The six facets of Figure (ref) depict the root mean squared error (RMSE) of the estimated covariance matrix and follow the same structure as the {six facets in Figure (ref).} We find that variance identification consistently reduces the estimated number of factors $r$ without affecting the quality of the estimated covariance matrix. In more than 50% of the dense cases, the posterior probability of the true number of factors is below 0.5 under the unrestricted scenario but over 0.5 under both restricted scenarios, which can be seen as an important jump. In the sparse setting, the posterior probabilities are generally lower, but the same pattern is observed. We also find that the necessary and the sufficient scenarios yield very similar results, which we read as the necessary and sufficient conditions being almost equivalent for our DGP's. In summary, we see that variance identification improves the estimate for the number of factors for all simulated data sets.
Results not reported here indicate the same conclusion under the uniform prior for {$\tau_j$}. In particular, variance identification improves the estimate for the number of factors without affecting the quality of the estimated covariance matrix. One difference is, however, that the posterior probabilities of the true number of factors are generally lower under the uniform prior than under the shrinkage prior. While the probabilities range even up to 0.7 under the shrinkage prior, as Figure (ref) shows, the largest ones are already below 0.04 in the $(50,10)$ dense setting and below 0.001 in most of the $(100,15)$ sparse settings under the uniform prior. The uniform prior {on $\tau_j$} does not provide as strong a signal for the correct number of factors as the shrinkage prior does, which is consistent with Fruehwirth2023Cusp.
Weekly returns of 17 currencies against the EUR are investigated between January, 2003, and December, 2005. The series include the currencies of big trading partners of the Eurozone (Australian Dollar, Canadian Dollar, British Pound, Hong Kong Dollar, Japanese Yen, South Korean Won, New Zealand Dollar, Russian Ruble, Turkish Lira, and US Dollar), and important local partners (Swiss Franc, Czech Koruna, Danish Krone, Norwegian Krone, Polish Zloty, Romanian Leu, and Swedish Krona). The chosen time period mostly avoids large international crises and heavy-tailed return distributions, as depicted in Figure (ref), which renders the static latent factor model (ref) appropriate for its analysis.
Estimation is done {using a moving window of} 52 weekly returns and the predictive performance is examined. In particular, the log posterior predictive likelihood $\text{LPPL}=\log\int_\theta p(y_{53}\mid\theta,y_1,\ldots,y_{52}) p(\theta\mid y_1,\ldots,y_{52})\,d\theta$ of the next weekly return is estimated as the natural logarithm of the mean of the sampled posterior predictive likelihoods $\mathcal{S}=\{p(y_{53}\mid\theta,y_1,\ldots,y_{52})\}_{\theta\in\text{posterior}}$, i.e., $\log((\sum_{s\in\mathcal{S}}s)/|\mathcal{S}|)$, where $\theta$ collects all parameters of the model, and $y_t$, $t=1,\ldots,53$, are the weekly returns for a given time window, including the next weekly return $y_{53}$. The sample means of the 52 weekly returns are subtracted from the input data before estimation and from the vector of next weekly returns before computing the LPPL. Then, the time window is shifted by one week, and estimation and prediction are repeated. The procedure is {repeated} 100 times, which covers approximately two years of weekly predictions under a moving window regime. ${H}=8$ is chosen for this exercise, which is the largest ${H}$ that satisfies ${H}\le(m-1)/2$ for $m=17$, and the same priors as in the simulation study are used. Importantly, we again consider two priors for {$\tau_j$} (shrinkage and uniform), which results in 200 posterior distributions in total for this exercise.
During post-processing, the three scenarios {regarding variance identification used in the simulation study} (unrestricted, necessary, and sufficient) are applied to the MCMC output. Two measures are computed for comparing the scenarios: the LPPL and the estimated number of factors. If model $\mathcal{M}_1$ has the same LPPL as model $\mathcal{M}_2$ but fewer factors, then $\mathcal{M}_1$ is preferred for its simplicity.
The dots in Figure (ref) show the LPPL of the sufficient scenario for a moving window of width 52 relative to the LPPL of the baseline unrestricted scenario, both under the shrinkage prior on $\tau_j$. The grey area represents the 5th to 95th percentiles of the posterior sample $\mathcal{S}$ used to estimate the LPPL under the unrestricted scenario, also relative to said baseline. We see that the LPPL of the sufficient scenario is very close to the LPPL of the unrestricted scenario as the difference stays close to zero. Moreover, the difference is considerably smaller than the width of the middle 90% region of the sampling distribution of the LPPL under the unrestricted scenario. Results not reported here show that both the uniform prior {on $\tau_j$} and the necessary scenario provide the same conclusion. In summary, restricting the prior to {variance} identified patterns does not significantly affect predictive performance of the factor model. Since this predictive measure purely depends on the estimated covariance matrix $\Omega$, this finding is consistent with the simulation study.
The top panel of Figure (ref) displays the shift in the posterior distribution $p(r\mid y_1,\ldots,y_{52})$ {across time} when switching from the unrestricted scenario to the sufficient scenario under the uniform prior {on $\tau_j$}. The bottom panel shows the same for the shrinkage prior. For instance, the blue triangle at the “2003-07/2004-06” label in the “Uniform prior” facet at $r=4$ denotes approximately 0.15, which means that the posterior probability of $r=4$ is 15 percentage points higher under the sufficient scenario than under the unrestricted scenario. Probabilities of large $r$ are generally reduced, and the probabilities of small $r$ are increased. We do not report results for the necessary scenario here, but the image is similar. The sea of downward-pointing triangles {lies} above the sea of upward-pointing triangles in both panels, which indicates that the sufficient scenario consistently reduces the estimated number of factors $r$ compared to the unrestricted scenario.
In our experience, the share of variance identified matrices increases in the posterior sample with more shrinkage, and this is reflected in Figure (ref), which shows the posterior proportion of variance identified ${\delta}$ matrices under the two prior specifications. The shrinkage prior prefers either close to empty or close to full columns a priori, separately for each column. In contrast, the uniform prior produces close to half full columns a priori. This spills over to the posterior distribution for this data set as can be seen from the proportions. The counting rule $\text{CR}({r},{1})$ is more likely satisfied with more crowded columns, which results in slightly higher acceptance rates in all time windows. Rates are mostly between 25% and 45%, and the difference between the two priors is consistent but not substantial. Further investigations not reported here show that the necessary scenario results in a similar increase in the proportion of variance identified matrices as the sufficient scenario does. Moreover, increasing shrinkage by decrasing $\alpha$ from five to three increases the distance between the two priors in the proportion of variance identified matrices, further supporting the conclusion that column shrinkage is beneficial for variance identification.
Overall, {both} the simulation study and the real world application consistently show that variance identification reduces the estimated number of factors without affecting the quality of the estimated covariance matrix. One drawback is increased computational time for the same number of draws, as parts of the MCMC output are discarded, but the {ensuing} reduction in efficiency is {small} compared to the benefits of an improved estimator.
{In this paper, we studied factor models which are a highly useful technique for dimension reduction in multivariate statistical analysis. To add to the mathematical understanding of these models, we focused on variance identification to uniquely identify the variance decomposition in the factor representation of a covariance matrix. We proved that a well-known counting rule based on the zero-nonzero pattern of the loading matrix is a sufficient condition for achieving variance identification. The proof relied on connecting factor analysis with some classical elements from graph and network theory which to our knowledge has not been exploited so far.}
{To enhance the relevance of this mathematical insight for practical factor analysis, we provide a computationally efficient algorithm for verifying the counting rule that again relies on results in graph and network theory.} Our methodology is illustrated for simulated as well as real data in the context of post-processing posterior draws in Bayesian sparse factor analysis. {As a main conclusion we find that certain inference tasks in factor analysis such as a predictive analysis are robust to whether posterior draws are variance identified, while others inference tasks such as identifying number of factors may be hugely impacted by the presence of unidentified posterior draws.}
We thank the Editor, Associate Editor and referees.