# Sparse Bayesian factor analysis when the number of factors is unknown

Sylvia Frühwirth-Schnatter111Department of Finance, Accounting, and
Statistics, WU Vienna University of Economics and Business, Austria. Email:
<EMAIL_ADDRESS>Darjus Hosszejni222Department of Finance, Accounting, and
Statistics, WU Vienna University of Economics and Business, Austria. Email:
<EMAIL_ADDRESS>Hedibert Freitas Lopes333School of Mathematical and
Statistical Sciences, Arizona State University, Tempe, USA & Insper Institute
of Education and Research, São Paulo, Brazil. Email<EMAIL_ADDRESS>
###### Abstract

There has been increased research interest in the subfield of sparse Bayesian
factor analysis with shrinkage priors, which achieve additional sparsity
beyond the natural parsimonity of factor models. In this spirit, we estimate
the number of common factors in the highly implemented sparse latent factor
model with spike-and-slab priors on the factor loadings matrix. Our framework
leads to a natural, efficient and simultaneous coupling of model estimation
and selection on one hand and model identification and rank estimation (number
of factors) on the other hand. More precisely, by embedding the unordered
generalised lower triangular loadings representation into overfitting sparse
factor modelling, we obtain posterior summaries regarding factor loadings,
common factors as well as the factor dimension via postprocessing draws from
our efficient and customised Markov chain Monte Carlo scheme.

Keywords: Hierarchical model; identifiability; point-mass mixture priors;
marginal data augmentation; reversible jump MCMC; prior distribution;
sparsity; Heywood problem; rotational invariance; ancillarity-sufficiency
interweaving strategy; fractional priors

JEL classification: C11, C38, C63

##  1 Introduction

Factor analysis aims at identifying common variation in multivariate
observations and relating it to hidden causes, the so-called common factors,
see Thurstone (1947) and, more recently, Anderson (2003). The common setup
consists of a sample
${\mathbf{y}}=\\{{\mathbf{y}}_{1},\ldots,{\mathbf{y}}_{T}\\}$ of $T$
multivariate observations
${\mathbf{y}}_{t}=(y_{1t},\ldots,y_{mt})^{{}^{\prime}}$ of dimension $m$. For
a given factor dimension $r$, the basic factor model is defined as a latent
variable model, involving the common factors ${\mathbf{f}}_{t}=(f_{1t}\cdots
f_{rt})^{{}^{\prime}}$:

$\displaystyle{\mathbf{f}}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},{{\mathbf{I}}}_{r}\right),\quad{\mathbf{y}}_{t}=\boldsymbol{\Lambda}{\mathbf{f}}_{t}+\boldsymbol{\epsilon}_{t},\quad\boldsymbol{\epsilon}_{t}\sim\mathcal{N}_{m}\left({\mathbf{0}},{\mathbf{\Sigma}}_{0}\right),\quad{\mathbf{\Sigma}}_{0}=\mbox{\rm
Diag}\,(\sigma^{2}_{1},\ldots,\sigma^{2}_{m}),$ (1.1)

where $\boldsymbol{\Lambda}$ is the $m\times r$ matrix of factor loadings
$\Lambda_{ij}$ and ${\mathbf{\Sigma}}_{0}$ is the covariance matrix of the
idiosyncratic errors $\boldsymbol{\epsilon}_{t}$. Model (1.1) implies that
conditional on ${\mathbf{f}}_{t}$ the $m$ elements of ${\mathbf{y}}_{t}$ are
independent and all dependence among these variables is explained through the
common factors. Assuming independence of ${\mathbf{f}}_{t}$ and
$\boldsymbol{\epsilon}_{t}$ implies that, marginally, ${\mathbf{y}}_{t}$
arises from a multivariate normal distribution,
${\mathbf{y}}_{t}\sim\mathcal{N}_{m}\left({\mathbf{0}},{\mathbf{\Omega}}\right)$,
with zero mean and a covariance matrix ${\mathbf{\Omega}}$ with the following
constrained structure:

$\displaystyle{\mathbf{\Omega}}=\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{{}^{\prime}}+{\mathbf{\Sigma}}_{0}.$
(1.2)

Since $r$ typically is (much) smaller than $m$, factor models yield a
parsimonious representation of ${\mathbf{\Omega}}$ with (at most) $m(r+1)$
instead of the $m(m+1)/2$ parameters of an unconstrained covariance matrix.
Hence, factor models proved to be very useful for covariance estimation,
especially if $m$ is large; see Fan et al. (2008), Forni et al. (2009),
Bhattacharya and Dunson (2011) and Kastner (2019), among others.

From the very beginning, the goal of factor analysis has been to extract the
hidden factors to understand the driving forces behind the observed
covariance. Using (1.2), this requires a decomposition of the covariance
matrix ${\mathbf{\Omega}}$ into the cross-covariance matrix
$\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{{}^{\prime}}$ and the covariance
matrix ${\mathbf{\Sigma}}_{0}$ of the uncorrelated idiosyncratic errors. With
the only source of information being the observed covariance of the data, this
is more challenging than estimating ${\mathbf{\Omega}}$ itself. A huge
literature, going back to Koopmans and Reiersøl (1950), Reiersøl (1950) and
Anderson and Rubin (1956) has addressed this problem of identification which
can be resolved by imposing additional structure on the factor model, see
Neudecker (1990), Geweke and Zhou (1996), Bai and Ng (2013), Chan et al.
(2018), and Williams (2020), among many others. Recently, a new identification
strategy based on unordered generalized lower triangular (UGLT) structures
(Frühwirth-Schnatter and Lopes, 2018; Frühwirth-Schnatter et al., 2022) was
introduced that allows to address not only the commonly known rotational
invariance problem, but also the problem of variance identification (Anderson
and Rubin, 1956) of which the literature is still less aware.

The main goal of the present paper is to foster mathematically rigorous
identification in Bayesian factor analysis. The recent years have seen many
contributions in the field of sparse Bayesian factor analysis which achieves
additional sparsity beyond the natural parsimonity of factor models through
the choice of shrinkage priors. One strand of literature considers continuous
shrinkage priors on the factor loadings, often in combination with a prior
that allows infinitely many columns in the loading matrix of which only a
finite number is non-zero such as the Indian buffet process prior (Griffiths
and Ghahramani, 2006) or the multiplicative Gamma process (Bhattacharya and
Dunson, 2011), see Zhao et al. (2016), Ročková and George (2017), Kastner
(2019), and Legramanti et al. (2020) for recent contributions. Alternatively,
following the pioneering paper by West (2003), many authors considered point
mass mixture priors on the factor loadings in basic factor models (Carvalho et
al., 2008; Frühwirth-Schnatter and Lopes, 2010), in dedicated factor models
with correlated factors (Conti et al., 2014) and dynamic factor models
(Kaufmann and Schuhmacher, 2019).

Sparse Bayesian factor analysis with point mass mixture (also called spike-
and-slab) priors allow factor loadings to be exactly zero and treat the
identification of these elements as a variable selection problem. This allows
to identify “simple structures” where in each row only a few nonzero loadings
are present (Anderson and Rubin, 1956), a long standing issue in factor
analysis addressed in Conti et al. (2014). It also allows to identify
irrelevant variables $y_{it}$ which are uncorrelated with the remaining
variables in ${\mathbf{y}}_{t}$, since the entire row of the factor loading
matrix is zero for these variables. Among other fields, this is of relevance
in economics (Kaufmann and Schuhmacher, 2017), where it is common practice to
include as many variables as possible in factor analysis (Stock and Watson,
2002; Boivin and Ng, 2006), and in bioinformatics (Lucas et al., 2006), where
typically only a few out of potentially ten thousands of genes may be related
to a certain physiological outcome.

A challenging problem in factor analysis is choosing the factor dimension $r$
which is usually unknown, see Owen and Wang (2016) for an excellent review.
Often information criteria (Bai and Ng, 2002) are used also in a Bayesian
context, see e.g. Aßmann et al. (2016) and Chan et al. (2018), other authors
employ marginal likelihoods (Lee and Song, 2002; Lopes and West, 2004).
Learning about the factor dimension $r$ is intrinsic in Bayesian approaches
that allow infinitely many columns in the loading matrix (Griffiths and
Ghahramani, 2006; Bhattacharya and Dunson, 2011; Ročková and George, 2017;
Legramanti et al., 2020), even if this approach requires careful tuning of
hyperparameters (Durante, 2017). A number of authors exploit variable
selection in a finite-dimensional overfitting factor model (Frühwirth-
Schnatter and Lopes, 2010; Conti et al., 2014; Kaufmann and Schuhmacher, 2017)
and we will follow their lead in the present paper. Our approach relies on the
following sparse Bayesian exploratory factor analysis (EFA) model:

$\displaystyle{\mathbf{y}}_{t}=\boldsymbol{\beta}_{k}{\mathbf{f}}^{k}_{t}+\boldsymbol{\epsilon}_{t},\qquad\boldsymbol{\epsilon}_{t}\sim\mathcal{N}_{m}\left({\mathbf{0}},{\mathbf{\Sigma}}_{k}\right),\quad{\mathbf{f}}^{k}_{t}\sim\mathcal{N}_{k}\left({\mathbf{0}},{{\mathbf{I}}}_{k}\right),$
(1.3)

where $\boldsymbol{\beta}_{k}$ is an $m\times k$ loading matrix with elements
$\beta_{ij}$ and ${\mathbf{\Sigma}}_{k}$ is a diagonal matrix with strictly
positive diagonal elements. Model (1.3) is potentially overfitting, since we
assume a finite factor dimension $k$ which is larger than true number of
factors $r$. We employ spike-and-slab priors, where the elements $\beta_{ij}$
of $\boldsymbol{\beta}_{k}$ are allowed to be exactly zero, with the
corresponding $m\times k$ sparsity matrix with elements $\delta_{ij}$ being
denoted by $\boldsymbol{\delta}_{k}$. To learn the factor dimension $r$, we
exploit a finite version of the two-parameter Beta prior (Ghahramani et al.,
2007; Frühwirth-Schnatter, 2022) to define a shrinkage process prior on
$\boldsymbol{\delta}_{k}$ that induces increasing shrinkage of the factor
loadings toward zero as the column index increases, extending recent work by
Legramanti et al. (2020).

To achieve identification, we impose a UGLT structure (Frühwirth-Schnatter et
al., 2022) on the loading matrix $\boldsymbol{\beta}_{k}$ and the
corresponding sparsity matrix $\boldsymbol{\delta}_{k}$ in the EFA model
(1.3). Commonly, in machine learning and statistics, no constraints are
imposed on $\boldsymbol{\delta}_{k}$; however, leaving the sparsity pattern
unconstrained makes it more difficult to recover the true number of factors
and reconstruct $\boldsymbol{\Lambda}$ from $\boldsymbol{\beta}_{k}$. In
econometrics as well as in statistics, $\boldsymbol{\delta}_{k}$ is often
constrained to a lower triangular matrix, however such a constraint is too
restrictive (Jöreskog, 1969; Carvalho et al., 2008). The UGLT structure that
we impose is a much weaker constraint which only requires the top non-zero
elements in each column of the loading matrix to lie in different rows and
still leaves the model unidentified. As shown in Frühwirth-Schnatter et al.
(2022), on the one hand it is weak enough to ensure that any loading matrix
can be rotated into a GLT representation, on the other hand it strong enough
to ensure “controlled unidentifiability” which can be easily resolved.

For practical Bayesian inference, we develop a new and efficient Markov chain
Monte Carlo (MCMC) procedure that delivers posterior draws from the EFA model
(1.3) under point mass mixture priors, which is known to be particularly
challenging (Pati et al., 2014). As part of our algorithm, we design a
(simple) reversible jump MCMC sampler to navigate through the space of UGLT
loading matrices of varying factor dimension. We achieve mathematically
rigorous identification in the spirit of Anderson and Rubin (1956) through
post-processing the posterior draws
$(\boldsymbol{\beta}_{k},{\mathbf{\Sigma}}_{k})$ from this unidentified EFA
model. By ensuring variance identification through the algorithm of Hosszejni
and Frühwirth-Schnatter (2022) and resolving rotational invariance during
postprocessing, we are able recover the factor dimension $r$, the
idiosyncratic variances ${\mathbf{\Sigma}}_{0}$ and an ordered GLT
representation $\boldsymbol{\Lambda}$ of the loading matrix from the posterior
draws. In this way, we add a mathematically rigorous approach to a growing
literature where commonly more heuristic post-processing procedures are
applied for this purpose; see (Aßmann et al., 2016; Kaufmann and Schuhmacher,
2019; Poworoznek et al., 2021; Papastamoulis and Ntzoufras, 2022).

Our sampling as well as our identification strategy work under arbitrary
choices regarding the slab distribution, including fractional priors
(Frühwirth-Schnatter and Lopes, 2010), the horseshoe prior (Zhao et al., 2016)
and the Lasso prior (Ročková and George, 2017) . In high-dimensional models,
we work with structured priors with column-specific shrinkage (Legramanti et
al., 2020) and employ the triple gamma prior (Cadonna et al., 2020) to achieve
good separation of signal and noise.

The rest of the paper is organized as follows. Section 2 introduces sparse
Bayesian EFA models with UGLT structures, while prior choices are discussed in
Section 3. Section 4 introduces our MCMC sampler for this model class and
Section 5 discusses post-processing MCMC draws to achieve identification.
Section 6 considers applications to exchange rate data and NYSE100 returns,
while Section 7 concludes.

##  2 Sparse Bayesian EFA models with UGLT structures

###  2.1 Model definition

Throughout the paper, we work with the exploratory factor analysis (EFA) model
(1.3), with $k$ potential common factors. Factor analysis based on this EFA
model yields the extended variance decomposition

$\displaystyle{\mathbf{\Omega}}=\boldsymbol{\beta}_{k}\boldsymbol{\beta}_{k}^{{}^{\prime}}+{\mathbf{\Sigma}}_{k},$
(2.1)

instead of the true variance decomposition (1.2) and the question arises how
to recover $(r,\boldsymbol{\Lambda},{\mathbf{\Sigma}}_{0})$ from
$(\boldsymbol{\beta}_{k},{\mathbf{\Sigma}}_{k})$. Without imposing constraints
on $\boldsymbol{\beta}_{k}$, this question is not easily answered, due to the
many identifiability issues in factor analysis. In recent work, Frühwirth-
Schnatter et al. (2022) discuss econometric identification of
$(r,\boldsymbol{\Lambda},{\mathbf{\Sigma}}_{0})$ from the sparse overfitting
Bayesian factor model (1.3), if the loading matrices are restricted to
unordered generalized lower triangular (UGLT) structures.

######  Definition 1.

UGLT structures. Let $\boldsymbol{\delta}_{r}$ be an binary matrix of $r$ non-
zero columns. Let $l_{j}$ denote the row index (also called pivot) of the top
non-zero entry in the $j$th column of $\boldsymbol{\delta}_{r}$ (i.e.
${\delta}_{ij}=0,\forall\,i<l_{j}$). $\boldsymbol{\delta}_{r}$ is said to be a
UGLT structure, if the pivot elements ${\mathbf{l}}_{r}=(l_{1},\ldots,l_{r})$
lie in different rows. More generally, a binary matrix
$\boldsymbol{\delta}_{k}$ with zero columns has an unordered GLT structure, if
the submatrix $\boldsymbol{\delta}_{r}$ containing the $r$ non-zero column of
$\boldsymbol{\delta}_{k}$ is a UGLT structure.

To facilitate identification in sparse Bayesian factor analysis, we assume in
the present paper that the factor loading matrix $\boldsymbol{\beta}_{k}$ and
the corresponding sparsity matrix $\boldsymbol{\delta}_{k}$ in the EFA model
(1.3) exhibit a UGLT structure. Compared to the common literature, where all
elements of $\boldsymbol{\delta}_{k}$ are left unspecified, this imposes the
constraint on $\boldsymbol{\delta}_{k}$ that the top non-zero element in the
various (non-zero) columns lie in different rows. As discussed in Frühwirth-
Schnatter et al. (2022), the (weak) UGLT constraint on
$\boldsymbol{\delta}_{k}$ is sufficient for a mathematically rigourous
identification even in an overfitting factor model, see Section 2.2 for more
details. These insights will be exploited in Section 5, where the posterior
draws from a sparse EFA model with UGLT structure are screened in a post-
processing manner to ensure identification of the posterior draws and to learn
about the unknown factor dimension $r$, the loading matrix
$\boldsymbol{\Lambda}$ as well as ${\mathbf{\Sigma}}_{0}$ from the data.

###  2.2 Econometric identification in factor models

Consider first an EFA model (1.3) that is not overfitting, i.e. $k=r$. A
rigorous approach toward identification of factor models was first offered by
Anderson and Rubin (1956) who consider identification as a two-step procedure.
The first step is variance identification, i.e. identification of
${\mathbf{\Sigma}}_{0}$ from the variance decomposition (1.2). Variance
identification is easily violated for sparse Bayesian factor models,
regardless whether $\boldsymbol{\delta}_{k}$ is unconstrained or exhibits a
UGLT structure. Frühwirth-Schnatter et al. (2022) prove that for UGLT
structures the so-called 3579 counting rule for the elements in the $r$ non-
zero columns $\boldsymbol{\delta}_{r}$ of $\boldsymbol{\delta}_{k}$ is
sufficient for variance identification. Hosszejni and Frühwirth-Schnatter
(2022) provide an efficient algorithm to verify this rule.

######  Definition 2.

3579 counting rule. A binary matrix $\boldsymbol{\delta}_{r}$ satisfies the
3579 counting rule, if the following condition is satisfied: for each
$q=1,\ldots,r$ and for each submatrix consisting of $q$ column of
$\boldsymbol{\delta}_{r}$, the number of nonzero rows in this sub-matrix is at
least equal to $2q+1$.

The 3579 counting rule states that every column of $\boldsymbol{\delta}_{r}$
should have at least 3, every pair of non-zero columns at least 5, every
subset of 3 columns at least 7 elements and so forth. If an indicator matrix
$\boldsymbol{\delta}_{r}$ obeys the 3579 counting rule, this implies that
${\mathbf{\Sigma}}_{0}={\mathbf{\Sigma}}_{r}$ and hence,
$\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{{}^{\prime}}=\boldsymbol{\beta}_{r}\boldsymbol{\beta}_{r}^{{}^{\prime}}$
is identified. A necessary condition for $\boldsymbol{\delta}_{r}$ to satisfy
the 3579 counting rule is the following upper bound for $r$:

$\displaystyle r\leq\frac{m-1}{2}.$ (2.2)

The second step of identification is solving the rotational invariance problem
provided that variance identification holds for
$(\boldsymbol{\beta}_{r},{\mathbf{\Sigma}}_{r})$. Since variance
identification implies that
$\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{{}^{\prime}}=\boldsymbol{\beta}_{r}\boldsymbol{\beta}_{r}^{{}^{\prime}}$,
it follows that $\boldsymbol{\beta}_{r}=\boldsymbol{\Lambda}{\mathbf{P}}$ for
some orthogonal matrix ${\mathbf{P}}$ (Anderson and Rubin, 1956, Lemma 5.1).
The rotational identification problem is usually solved by imposing a
structure on the loading matrices $\boldsymbol{\beta}_{r}$ that ensures
identification of $\boldsymbol{\Lambda}$ from
$\boldsymbol{\beta}_{r}\boldsymbol{\beta}_{r}^{{}^{\prime}}$ in the sense that
$\boldsymbol{\beta}_{r}=\boldsymbol{\Lambda}{\mathbf{P}}$ iff ${\mathbf{P}}$
is equal to the identity. A weaker condition is rotational identification up
to signed permutations
$\boldsymbol{\beta}_{r}=\boldsymbol{\Lambda}{\mathbf{P}}_{\pm}{\mathbf{P}}_{\rho}$,
where the permutation matrix ${\mathbf{P}}_{\rho}$ corresponds to one of the
$r$! permutations and the reflection matrix ${\mathbf{P}}_{\pm}=\mbox{\rm
Diag}\,(\pm 1,\ldots,\pm 1)$ corresponds to one of the $2^{r}$ ways to switch
the signs of the $r$ columns of $\boldsymbol{\Lambda}$. Frühwirth-Schnatter et
al. (2022) show that imposing a UGLT structure on $\boldsymbol{\beta}_{r}$ and
$\boldsymbol{\Lambda}$ leads to rotational identification up to signed
permutations. Provided that $\boldsymbol{\beta}_{r}$ is variance identified,
$\boldsymbol{\Lambda}$ is recovered from $\boldsymbol{\beta}_{r}$ by
reordering the columns of $\boldsymbol{\beta}_{r}$ such that the pivots
$l_{1}<\ldots<l_{r}$ are increasing.

In applied factor analysis, the EFA model (1.3) is typically overfitting with
$k>r$ and additional identifiability issues arise. In an overfitting EFA
model, identifiability of ${\mathbf{\Sigma}}_{k}$ and
$\boldsymbol{\beta}_{k}\boldsymbol{\beta}_{k}^{{}^{\prime}}$ from (2.1) is
lost and infinitely many representations
$(\boldsymbol{\beta}_{k},{\mathbf{\Sigma}}_{k})$ with
${\mathbf{\Sigma}}_{k}\neq{\mathbf{\Sigma}}_{0}$ exist that imply the same
covariance matrix ${\mathbf{\Omega}}$ as
$(\boldsymbol{\Lambda},{\mathbf{\Sigma}}_{0})$ (Geweke and Singleton, 1980;
Tumura and Sato, 1980). Frühwirth-Schnatter et al. (2022) show that imposing a
UGLT structure on $\boldsymbol{\beta}_{k}$ allows to reveal
$\boldsymbol{\Lambda}$ from $\boldsymbol{\beta}_{k}$ even in this case.

Consider, for illustration, an overfitting EFA model with $k=r+1$. Factor
analysis typically yields loading matrices $\boldsymbol{\beta}_{r+1}$ with
$r+1$ non-zero columns rather than matrices with $r$ non-zero and a single
zero column. Frühwirth-Schnatter et al. (2022, Theorem 7) prove that by
imposing a UGLT structure on $\boldsymbol{\beta}_{r+1}$, the additional column
in the overfitting model is equal to a so-called _spurious factor_
$\boldsymbol{\Xi}$, e.g.:

$\displaystyle{\small\boldsymbol{\beta}_{r+1}=\left(\begin{array}[]{cc}\boldsymbol{\Lambda}&\boldsymbol{\Xi}\end{array}\right),\quad\boldsymbol{\Xi}=\left(\begin{array}[]{c}{\mathbf{0}}\\\
\Xi_{l_{\footnotesize\rm sp}}\\\
{\mathbf{0}}\end{array}\right),\quad{\mathbf{\Sigma}}_{r+1}=\mbox{\rm
Diag}\,(\sigma^{2}_{1},\ldots,\sigma^{2}_{l_{\footnotesize\rm
sp}}-\Xi_{l_{\footnotesize\rm sp}}^{2},\ldots,\sigma^{2}_{m}),}$ (2.7)

with a single non-zero factor loading $\Xi_{l_{\footnotesize\rm sp}}$
satisfying $0<\Xi_{l_{\footnotesize\rm
sp}}^{2}<\sigma^{2}_{l_{\footnotesize\rm sp}}$ which lies in a pivot row
$l_{\footnotesize\rm sp}$ different from the pivot rows $l_{1},\ldots,l_{r}$
in $\boldsymbol{\Lambda}$.

This result allows to identify spurious columns $\boldsymbol{\Xi}$ in the
loading matrix $\boldsymbol{\beta}_{r+1}$ of the overfitting EFA model (1.3)
and to recover $\boldsymbol{\Lambda}$ from the remaining (active) columns up
to a signed permutation,
$\boldsymbol{\beta}_{r}=\boldsymbol{\Lambda}{\mathbf{P}}_{\pm}{\mathbf{P}}_{\rho}$.
The covariance matrix
${\mathbf{\Sigma}}_{0}={\mathbf{\Sigma}}_{r+1}+\boldsymbol{\Xi}\boldsymbol{\Xi}^{{}^{\prime}}$
is identified by moving the spurious column $\boldsymbol{\Xi}$ to the
idiosyncratic errors. This results also holds for a higher degree of
overfitting and is the cornerstone of our MCMC sampler which relates
exploratory and confirmatory Bayesian factor analysis.

###  2.3 Relating exploratory to confirmatory Bayesian factor analysis

The sparsity matrix $\boldsymbol{\delta}_{k}$ in the EFA model (1.3) allows to
identify which factors are active (the corresponding column of
$\boldsymbol{\delta}_{k}$ has at least two non-zero loading), which factors
are spurious (the corresponding column of $\boldsymbol{\delta}_{k}$ has a
single non-zero loading), and which ones are inactive (the corresponding
column of $\boldsymbol{\delta}_{k}$ is zero). This allows to split the
sparsity matrix $\boldsymbol{\delta}_{k}$ into an $m\times r$ submatrix
$\boldsymbol{\delta}_{r}$ with $r$ active columns, an $m\times
r_{\footnotesize sp}$ submatrix $\boldsymbol{\delta}_{\Xi}$ with
$r_{\footnotesize sp}$ spurious columns, and a submatrix with
$j_{\footnotesize 0}=k-r-r_{\footnotesize sp}$ zero columns. The loading
matrix $\boldsymbol{\beta}_{k}$ is split accordingly into the $m\times r$
submatrix $\boldsymbol{\beta}_{r}$, the $m\times r_{\footnotesize sp}$
submatrix $\boldsymbol{\Xi}$, and $j_{\footnotesize 0}$ zero columns, while
the factors ${\mathbf{f}}_{t}^{k}$ are split into ${\mathbf{f}}_{t}^{r}$,
${\mathbf{f}}_{t}^{\Xi}$ and ${\mathbf{f}}^{0}_{t}$.

Exploiting representation (2.7), we extract the following model of factor
dimension $r$ which is embedded in any EFA model with UGLT structure,

$\displaystyle{\mathbf{f}}^{r}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},{{\mathbf{I}}}_{r}\right),\quad{\mathbf{y}}_{t}=\boldsymbol{\beta}_{r}{\mathbf{f}}^{r}_{t}+\boldsymbol{\epsilon}_{t},\quad\boldsymbol{\epsilon}_{t}\sim\mathcal{N}_{m}\left({\mathbf{0}},{\mathbf{\Sigma}}_{r}\right),\quad{\mathbf{\Sigma}}_{r}={\mathbf{\Sigma}}_{k}+\boldsymbol{\Xi}\boldsymbol{\Xi}^{{}^{\prime}},$
(2.8)

by moving the $r_{\footnotesize sp}$ spurious columns $\boldsymbol{\Xi}$ to
the idiosyncratic variances ${\mathbf{\Sigma}}_{r}$. We call (2.8) the
confirmatory factor analysis (CFA) model induced by the active columns of the
sparsity matrix $\boldsymbol{\delta}_{k}$ in the EFA model. The likelihood
function is invariant to moving from the EFA model to the CFA model (2.8),
since the implied covariance matrix
${\mathbf{\Omega}}=\boldsymbol{\beta}_{k}\boldsymbol{\beta}_{k}^{{}^{\prime}}+{\mathbf{\Sigma}}_{k}=\boldsymbol{\beta}_{r}\boldsymbol{\beta}_{r}^{{}^{\prime}}+{\mathbf{\Sigma}}_{r}$
remains the same. On the other hand, taking the CFA model (2.8) as a starting
point, we can move to the EFA model (1.3) without changing the likelihood
function by adding $r_{\footnotesize sp}$ spurious columns
$\boldsymbol{\delta}_{\Xi}$ to $\boldsymbol{\delta}_{r}$. The only relevant
information needed for expanding the CFA model in this way is that
$r_{\footnotesize sp}$ columns outside of $\boldsymbol{\delta}_{r}$ have
column-size one, whereas all remaining columns are zero.

For $r_{\footnotesize sp}=1$, for instance, a single spurious column
$\boldsymbol{\delta}_{\Xi}$ is added to $\boldsymbol{\delta}_{r}$ to define an
EFA model with $k-r-1$ zero columns. The position of the only non-zero
indicator in column $\boldsymbol{\delta}_{\Xi}$, denoted by
$\delta^{\Xi}_{l_{\footnotesize\rm sp}}$, is not identified and can lie in any
row $l_{\footnotesize\rm sp}$ that is different from the pivots
${\mathbf{l}}_{r}$ in $\boldsymbol{\delta}_{r}$. A spurious column
$\boldsymbol{\Xi}$ is added to $\boldsymbol{\beta}_{r}$ to define
$\boldsymbol{\beta}_{k}$, while the covariance matrix of the idiosyncratic
errors in the EFA model is defined as
${\mathbf{\Sigma}}_{k}={\mathbf{\Sigma}}_{r}-\boldsymbol{\Xi}\boldsymbol{\Xi}^{{}^{\prime}}$.
The only non-zero loading $\Xi_{l_{\footnotesize\rm sp}}$ in
$\boldsymbol{\Xi}$ can take any value such that the $l_{\footnotesize\rm
sp}$-th diagonal element of ${\mathbf{\Sigma}}_{k}$ remains positive, i.e.
${\mathbf{\Sigma}}_{k,l_{\footnotesize\rm sp},l_{\footnotesize\rm
sp}}=\sigma^{2}_{l_{\footnotesize\rm sp}}-(\Xi_{l_{\footnotesize\rm
sp}})^{2}>0$. This entire move only affects the $l_{\footnotesize\rm sp}$-th
row $\boldsymbol{\beta}_{r,l_{\footnotesize\rm sp},\cdot}$ of
$\boldsymbol{\beta}_{r}$. More specifically, for $t=1,\ldots,T$:

$\displaystyle y_{l_{\footnotesize\rm
sp},t}=\boldsymbol{\beta}_{r,l_{\footnotesize\rm
sp},\cdot}{\mathbf{f}}^{r}_{t}+\epsilon_{l_{\footnotesize\rm
sp},t},\,\epsilon_{l_{\footnotesize\rm
sp},t}\sim\mathcal{N}\left(0,\sigma^{2}_{l_{\footnotesize\rm sp}}\right),$
(2.9) $\displaystyle y_{l_{\footnotesize\rm
sp},t}=\boldsymbol{\beta}_{r,l_{\footnotesize\rm
sp},\cdot}{\mathbf{f}}^{r}_{t}+\Xi_{l_{\footnotesize\rm
sp}}f_{t}^{\Xi}+\tilde{\epsilon}_{l_{\footnotesize\rm
sp},t},\,\tilde{\epsilon}_{l_{\footnotesize\rm
sp},t}\sim\mathcal{N}\left(0,(\sigma^{2}_{l_{\footnotesize\rm
sp}})^{\mbox{\rm\tiny new}}\right),\,(\sigma^{2}_{l_{\footnotesize\rm
sp}})^{\mbox{\rm\tiny new}}=\sigma^{2}_{l_{\footnotesize\rm
sp}}-(\Xi_{l_{\footnotesize\rm sp}})^{2}.$

By integrating model (2.9) with respect to the spurious factor $f_{t}^{\Xi}$,
it can once more be verified that both models imply the same distribution
$p(y_{l_{\footnotesize\rm sp},t}|\boldsymbol{\beta}_{r,l_{\footnotesize\rm
sp},\cdot},{\mathbf{f}}^{r}_{t},\sigma^{2}_{l_{\footnotesize\rm sp}})$,
independently of the specific value of $\Xi_{l_{\footnotesize\rm sp}}$. Also
for arbitrary $r_{\footnotesize sp}\in\\{1,\ldots,k-r\\}$, neither the
position of the pivots ${\mathbf{l}}_{\Xi}$ of $\boldsymbol{\delta}_{\Xi}$
(which are the only non-zero indicator in each column) nor the non-zero
elements in the corresponding spurious loading matrix $\boldsymbol{\Xi}$ are
identified. However, the pivots in ${\mathbf{l}}_{\Xi}$ must lie in different
rows.

Moving forth and back between the EFA model (1.3) and the CFA model (2.8) as
described above is the cornerstone of an efficient MCMC algorithm developed in
Section 4 that operates in the space of UGLT matrices with varying dimension
without imposing further constraints. In Section 3, suitable priors are
defined for a sparse EFA model with UGLT structure that are (largely)
invariant to these moves.

##  3 Prior specifications

###  3.1 Column sparsity through exchangeable shrinkage process priors

Bayesian inference is performed in the EFA model (1.3) with a finite number
$k$ of potential factors. Dirac-spike-and-slab prior for the factor loadings
are assumed,

$\displaystyle\beta_{ij}|\tau_{j}\sim(1-\tau_{j})\Delta_{0}+\tau_{j}P_{\tiny\rm
slab}(\beta_{ij}),$ (3.1)

where the columns of the loading matrix are increasingly pulled toward 0 as
the column index increases. This cumulative shrinkage is achieved indirectly
by placing an exchangeable shrinkage process (ESP) prior on the slab
probabilities $\tau_{1},\ldots,\tau_{k}$:

$\displaystyle\tau_{j}|k\sim\mathcal{B}\left(a_{k},b_{k}\right),\quad
j=1,\ldots,k,$ (3.2)

see (Frühwirth-Schnatter, 2022). The ESP prior turns model (1.3) into a sparse
EFA model, where the number $r$ of active columns in $\boldsymbol{\delta}_{k}$
with at least two non-zero elements is a random variable which takes values
smaller than $k$ with positive probability. As recently shown by Frühwirth-
Schnatter (2022), prior (3.2) has a representation as a finite cumulative
shrinkage process (CUSP) prior (Legramanti et al., 2020). Prominent examples
of such ESP priors are the finite one-parameter-beta (1PB) prior:

$\displaystyle\tau_{j}|k\sim\mathcal{B}\left(\frac{\alpha}{k},1\right),\quad
j=1,\ldots,k,$ (3.3)

and the finite two-parameter-beta (2PB) prior,

$\displaystyle\tau_{j}|k\sim\mathcal{B}\left(\gamma\frac{\alpha}{k},\gamma\right),\quad
j=1,\ldots,k.$ (3.4)

Alternative choices are $a_{k}=\alpha/k$ and $b_{k}=\gamma(k-1)/k$ (Paisley
and Carin, 2009). For $k\rightarrow\infty$, prior (3.3) converges to the
Indian buffet process prior (Teh et al., 2007) and has been employed by
Ročková and George (2017) in sparse Bayesian factor analysis. For
$k\rightarrow\infty$, prior (3.4) converges to the two-parameter-beta prior
introduced by Ghahramani et al. (2007) in Bayesian nonparametric latent
feature models which can be regarded as a factor model with infinitely many
columns of which only a finite number is non-zero.

As opposed to this literature, we stay within the framework of factor models
with finitely many columns. While obeying the upper bound (2.2), $k$ is
selected large enough to encourage spurious columns and zero columns in
$\boldsymbol{\delta}_{k}$. Spurious columns are essential for our strategy of
recovering the factor dimension $r$ through adding and deleting spurious
columns in the overfitting EFA model (1.3). The hyperparameters $\alpha$ and
$\gamma$ are instrumental in controlling prior column sparsity. A prior with
$\alpha<k$ and $\gamma=1$, e.g., induces sparsity, since the largest slab
probability $\tau_{(k)}\sim\mathcal{B}\left(\alpha,1\right)$, while the
smallest slab probability $\tau_{(1)}\sim\mathcal{B}\left(\alpha/k,1\right)$
(Frühwirth-Schnatter, 2022). To adapt the hyperparameters to the data at hand,
$\alpha$ is assumed to be a random with prior
$\alpha\sim\mathcal{G}\left(a^{\alpha},b^{\alpha}\right)$,while $\gamma$ is
random with prior $\gamma\sim\mathcal{G}\left(a^{\gamma},b^{\gamma}\right)$
for the 2PB-prior (3.4).

#### Imposing a UGLT structure

For given numbers $r$ and $r_{\footnotesize sp}$ of, respectively, active and
spurious columns, we define a prior on the pivots
${\mathbf{l}}_{r}=(l_{1},\ldots,l_{r})$ and
${\mathbf{l}}_{\Xi}=(l_{{\Xi},1},\ldots,l_{{\Xi},r_{\footnotesize sp}})$ such
that the non-zero columns of the sparsity matrix $\boldsymbol{\delta}_{k}$
exhibit a UGLT structure. The prior $p({\mathbf{l}}_{r})$ in the CFA model is
defined as follows. Let ${\cal
L}({\mathbf{l}})=\\{i\in\\{1,2,\ldots,m\\}:i\notin{\mathbf{l}}\\}$. Given the
pivots ${\mathbf{l}}_{r,-j}$ outside of any column $j$, condition UGLT implies
that $l_{j}$ has to be different from ${\mathbf{l}}_{r,-j}$, and we assume a
uniform prior distribution over all admissible pivots $l_{j}\in{\cal
L}({\mathbf{l}}_{r,-j})$:

$\displaystyle p(l_{j}|{\mathbf{l}}_{r,-j})=\frac{1}{|{\cal
L}({\mathbf{l}}_{r,-j})|}=\frac{1}{m-r+1}.$ (3.5)

The pivots ${\mathbf{l}}_{\Xi}$ of the spurious columns
$\boldsymbol{\delta}_{\Xi}$ are assigned a uniform prior over all admissible
values, given the pivots ${\mathbf{l}}_{r}$, using the following order-
independent construction. Given ${\mathbf{l}}_{r}$, $l_{{\Xi},1}$ is uniform
over ${\cal L}({\mathbf{l}}_{r})$; given $l_{{\Xi},1}$, $l_{{\Xi},2}$ is
uniform over ${\cal L}({\mathbf{l}}_{r}\cup\\{l_{{\Xi},1}\\})$, and so forth.

Given the pivots $l_{j}$ in all non-zero columns of $\boldsymbol{\delta}_{k}$,
by definition $\delta_{l_{j},j}=1$ and $\delta_{ij}=0$ for $i<l_{j}$, while
the $m-l_{j}$ indicators $\delta_{ij}$ below $l_{j}$ are subject to variable
selection with column-specific probabilities $\tau_{j}$ following the ESP
prior (3.2):

$\displaystyle\mbox{\rm
Pr}(\delta_{ij}=1|l_{j},\tau_{j})=\left\\{\begin{array}[]{ll}0,&i<l_{j},\\\
1,&i=l_{j},\\\ \tau_{j},&i=l_{j}+1,\ldots,m.\\\ \end{array}\right.$ (3.9)

Let $d_{j}=\sum_{i=1}^{m}{\delta_{ij}}$ be the number of non-zero indicators
in columns $j$. With $d_{j}-1$ successes and $m-l_{j}-d_{j}+1$ failures in the
experiment defined in (3.9), we obtain the following prior for column
$\boldsymbol{\delta}^{r}_{\cdot,j}$:

$\displaystyle\mbox{\rm
Pr}(\boldsymbol{\delta}^{r}_{\cdot,j}|l_{j},\tau_{j})=\tau_{j}^{d_{j}-1}(1-\tau_{j})^{m-l_{j}-d_{j}+1}.$
(3.10)

If we integrate over $\tau_{j}$, then we obtain:

$\displaystyle\mbox{\rm
Pr}(\boldsymbol{\delta}^{r}_{\cdot,j}|l_{j})=\frac{B(a_{k}+d_{j}-1,b_{k}+m-l_{j}-d_{j}+1)}{B(a_{k},b_{k})}.$
(3.11)

###  3.2 Choosing the slab distribution

To define a prior on the loading matrix $\boldsymbol{\beta}_{k}$, we split
$\boldsymbol{\beta}_{k}$ as discussed in Section 2.3 and define a prior
$p(\boldsymbol{\beta}_{r}|{\mathbf{\Sigma}}_{r},\boldsymbol{\delta}_{r})$ on
the factor loading matrix $\boldsymbol{\beta}_{r}$ in the CFA model (2.8),
conditional on ${\mathbf{\Sigma}}_{r}=\mbox{\rm
Diag}\,(\sigma^{2}_{1},\ldots,\sigma^{2}_{m})$ and $\boldsymbol{\delta}_{r}$.
When expanding the CFA model to an EFA model with $r_{\footnotesize sp}$
columns, we define a prior
$p(\boldsymbol{\Xi}|\boldsymbol{\beta}_{r},{\mathbf{\Sigma}}_{r},{\mathbf{l}}_{\Xi})$
on the spurious loadings conditional on $\boldsymbol{\beta}_{r}$,
${\mathbf{\Sigma}}_{r}$, and ${\mathbf{l}}_{\Xi}$. The spurious factor
loadings are assigned a uniform prior over all values that lead to a positive
definite matrix
${\mathbf{\Sigma}}_{k}={\mathbf{\Sigma}}_{r}-\boldsymbol{\Xi}\boldsymbol{\Xi}^{{}^{\prime}}$
in the EFA model:

$\displaystyle\Xi_{l_{\footnotesize\rm
sp}}^{2}|\sigma^{2}_{l_{\footnotesize\rm
sp}}\sim\mathcal{U}\left[0,\sigma^{2}_{l_{\footnotesize\rm sp}}\right].$
(3.12)

This ensures for all $l_{\footnotesize\rm sp}\in{\mathbf{l}}_{\Xi}$ that
${\mathbf{\Sigma}}_{k,l_{\footnotesize\rm sp},l_{\footnotesize\rm
sp}}=\sigma^{2}_{l_{\footnotesize\rm sp}}-\Xi_{l_{\footnotesize\rm
sp}}^{2}>0$. By this definition, both the likelihood and the prior are
invariant to moving between the EFA and the CFA model for a given number of
spurious columns $r_{\footnotesize sp}$, regardless of the chosen slab
distribution.444Note that this is a major improvement compared to Frühwirth-
Schnatter and Lopes (2018).

The spike-and-slab prior (3.1) is formulated for the factor loading matrix
$\boldsymbol{\beta}_{r}$ in the CFA model. A broad range of slab distributions
has been considered for sparse Bayesian factor analysis and can be combined
with the RJMCMC sampler we introduce in Section 4. Since the conditional
likelihood function factors into a product over the rows of the loading
matrix, prior independence across the rows is assumed. For each row $i$ of
$\boldsymbol{\delta}_{r}$ with $q_{i}=\sum_{j}\delta_{ij}>0$ nonzero elements,
we consider two different prior families for the vector
$\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}$ of unconstrained elements
in this row. First, hierarchical Gaussian priors are considered, taking the
form

$\displaystyle\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i}\sim\mathcal{N}_{q_{i}}\left({\mathbf{0}},{\mathbf{B}}_{i0}^{\boldsymbol{\delta}}\sigma^{2}_{i}\right),$
(3.13)

where ${\mathbf{B}}_{i0}^{\boldsymbol{\delta}}$ is a diagonal matrix. The
variance of prior (3.13) is assumed to depend on $\sigma^{2}_{i}$, because
this allows joint drawing of $\boldsymbol{\beta}_{r}$ and
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$ and, even more importantly, sampling
the sparsity matrix $\boldsymbol{\delta}_{r}$ without conditioning on the
model parameters during MCMC estimation, see Algorithm 1. For models with a
moderate number of features $m$, the choice
${\mathbf{B}}_{i0}^{\boldsymbol{\delta}}=A_{0}{{\mathbf{I}}}_{q_{i}}$ with a
fixed hyperparameter $A_{0}$ is common (Lopes and West, 2004; Ghosh and
Dunson, 2009; Conti et al., 2014). For high-dimensional models with larger
$m$, often a structured hierarchical prior is assumed to achieve better
shrinkage of the factor loadings. One example are spike-and-slab priors with a
column specific shrinkage parameter $\theta_{j}$ in the slab distribution of
$\beta_{ij}$, see e.g. Kaufmann and Schuhmacher (2019), Legramanti et al.
(2020) and Frühwirth-Schnatter (2022):

$\displaystyle\beta_{ij}|\delta_{ij}=1,\theta_{j},\sigma^{2}_{i},\kappa\sim\mathcal{N}\left(0,\kappa\theta_{j}\sigma^{2}_{i}\right).$
(3.14)

The column specific shrinkage parameters follow an exchangeable prior, either
an inverse gamma prior,
$\theta_{j}\sim\mathcal{G}^{-1}\left(c^{\theta},b^{\theta}\right)$ (Legramanti
et al., 2020), or a triple gamma prior, $\theta_{j}\sim\mbox{\rm
F}\left(2a^{\theta},2c^{\theta}\right)$ (Cadonna et al., 2020). The triple
gamma prior acts as a variance selection prior and puts considerable mass on
small values of $\theta_{j}$ which simultaneously pulls the factors $f_{jt}$
toward 0 for all time points $t$. For $a^{\theta}=c^{\theta}=0.5$, this yields
a grouped version of the horseshoe prior (Zhao et al., 2016). The global
shrinkage parameter follows
$\kappa\sim\mathcal{G}^{-1}\left(c^{\kappa},b^{\kappa}\right)$ or
$\kappa\sim\mbox{\rm F}\left(2a^{\kappa},2c^{\kappa}\right)$. Prior (3.14) can
be extended by assuming local shrinkage parameters $\omega_{ij}$ arising from
an F-distribution for each non-zero factor loading:

$\displaystyle\beta_{ij}|\delta_{ij}=1,\omega_{ij},\theta_{j},\sigma^{2}_{i},\kappa\sim\mathcal{N}\left(0,\kappa\theta_{j}\sigma^{2}_{i}\omega_{ij}\right),\quad\omega_{ij}\sim\mbox{\rm
F}\left(2a^{\omega},2c^{\omega}\right).$ (3.15)

For $a^{\omega}=c^{\omega}=0.5$, reduces to the horseshoe prior employed by
Zhao et al. (2016).

As an alternative shrinkage prior, Frühwirth-Schnatter and Lopes (2010)
introduced a conditionally conjugate fractional prior
$p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r})\propto\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{b}$
in the spirit of O’Hagan (1995). It can be interpreted as the posterior of a
non-informative prior and a small fraction $b>0$ the conditional likelihood is
derived from regression model

$\displaystyle\tilde{{\mathbf{y}}}_{i}={\mathbf{X}}_{i}^{\boldsymbol{\delta}}\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}+\tilde{\boldsymbol{\epsilon}}_{i},$
(3.16)

where $\tilde{{\mathbf{y}}}_{i}=(y_{i1}\cdots y_{iT})^{{}^{\prime}}$,
$\tilde{\boldsymbol{\epsilon}}_{i}=(\epsilon_{i1}\cdots\epsilon_{iT})^{{}^{\prime}}$
and ${\mathbf{X}}_{i}^{\boldsymbol{\delta}}$ is a regressor matrix constructed
from the latent factors
${\mathbf{f}}_{r}=\\{{\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T}\\}$ (see
Appendix A.1 for details). This yields:

$\displaystyle\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r}\sim\mathcal{N}_{q_{i}}\left({\mathbf{b}}_{iT}^{\boldsymbol{\delta}},{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}\sigma^{2}_{i}/b_{N}\right),$
(3.17)

where
${\mathbf{b}}_{iT}^{\boldsymbol{\delta}}={\mathbf{B}}_{iT}^{\boldsymbol{\delta}}({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i}$
and
${\mathbf{B}}_{iT}^{\boldsymbol{\delta}}=(({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{-1}$.
Similar fractional priors have been applied for variable selection in other
latent variable models (Frühwirth-Schnatter and Tüchler, 2008; Frühwirth-
Schnatter and Wagner, 2010). In the spirit of Foster and George (1994),
$b_{N}=1/N$, since $N=mT$ observations are available to estimate
$\boldsymbol{\beta}_{r}$. The fractional prior often works better than
hierarchical Gaussian priors for moderate $m$, however is not flexible enough
to handle larger dimensions $m$.

###  3.3 The prior on the idiosyncratic variances

Finally, we define a prior on the idiosyncratic variances
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$, independently of all other parameters.
When estimating factor models using classical statistical methods such as ML
estimation, it frequently happens that the optimal solution lies outside the
admissible parameter space with one or more of the idiosyncratic variances
$\sigma^{2}_{i}$s being negative, see e.g. (Jöreskog, 1967; Bartholomew,
1987). This difficulty became known as the Heywood problem. Using a prior on
the idiosyncratic variances $\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$ within a
Bayesian framework, typically chosen from the inverted Gamma family,

$\displaystyle\sigma^{2}_{i}\sim\mathcal{G}^{-1}\left(c_{0},C_{i0}\right),$
(3.18)

naturally avoids negative values for $\sigma^{2}_{i}$. Nevertheless, there
exists a Bayesian analogue of the Heywood problem which takes the form of
multi-modality of the posterior of $\sigma^{2}_{i}$ with one mode lying at 0.
Heywood problems typically occur, if the constraint

$\displaystyle\frac{1}{\sigma^{2}_{i}}\geq({\mathbf{\Omega}}^{-1})_{ii}\quad\Leftrightarrow\quad\sigma^{2}_{i}\leq\frac{1}{({\mathbf{\Omega}}^{-1})_{ii}}$
(3.19)

is violated, where ${\mathbf{\Omega}}$ is the covariance matrix of
${\mathbf{y}}_{t}$ defined in (1.2), see e.g. Bartholomew (1987, p. 54). It is
clear from inequality (3.19) that $1/\sigma^{2}_{i}$ has to be bounded away
from 0. For this reason, improper priors on the idiosyncratic variances such
as $p(\sigma^{2}_{i})\propto 1/\sigma^{2}_{i}$ (Martin and McDonald, 1975;
Akaike, 1987) are not able to prevent Heywood problems. Similarly, proper
inverted Gamma priors with small degrees of freedom such as $c_{0}=1.1$ (Lopes
and West, 2004) allow values too close to 0. Hence, the hyperparameters
$c_{0}$ and $C_{i0}$ in (3.18) have to be selected in such a way that Heywood
problems are avoided. In particular, $c_{0}$ should be chosen large enough to
bound the prior away from 0, typically $c_{0}=2.5$. Regarding $C_{i0}$, the
most common choice is using a fixed value $C_{i0}=C_{0}$.

An alternative choice has been suggested in Frühwirth-Schnatter and Lopes
(2010). They reduce the occurrence probability of a Heywood problem which is
equal to $\mbox{\rm Pr}(X\leq C_{i0}({\mathbf{\Omega}}^{-1})_{ii})$ where
$X\sim\mathcal{G}\left(c_{0},1\right)$ through the choice of individual
scalings $C_{i0}$. The smaller $C_{i0}$, the smaller is this probability. On
the other hand, a downward bias is introduced, if $C_{i0}$ is too small, since
$\mbox{\rm E}(\sigma^{2}_{i})=C_{i0}/(c_{0}-1)$. Choosing
$C_{i0}=(c_{0}-1)/(\widehat{{\mathbf{\Omega}}^{-1}})_{ii}$ as the largest
value for which inequality (3.19) is fulfilled by the prior expectation
$\mbox{\rm E}(\sigma^{2}_{i})$ yields the following prior:

$\displaystyle\sigma^{2}_{i}\sim\mathcal{G}^{-1}\left(c_{0},(c_{0}-1)/(\widehat{{\mathbf{\Omega}}^{-1}})_{ii}\right),$
(3.20)

based on an estimator $\widehat{{\mathbf{\Omega}}^{-1}}$ of
${\mathbf{\Omega}}^{-1}$. In Frühwirth-Schnatter and Lopes (2018), a Bayesian
estimator is proposed which combines the sample information with an inverted
Wishart prior
${\mathbf{\Omega}}^{-1}\sim\mathcal{W}_{m}\left(\nu_{o},\nu_{o}{\mathbf{S}}_{o}\right)$:

$\displaystyle\widehat{{\mathbf{\Omega}}^{-1}}=(\nu_{o}+T/2)(\nu_{o}{\mathbf{S}}_{o}+0.5\sum_{t=1}^{T}{\mathbf{y}}_{t}{\mathbf{y}}_{t}^{{}^{\prime}})^{-1}.$
(3.21)

##  4 MCMC estimation

We use MCMC techniques to sample from the posterior distribution of the EFA
model (1.3), given the priors introduced in Section 3. As noted by many
authors, e.g. Pati et al. (2014), MCMC sampling for sparse Bayesian factor
models is notoriously difficult, since sampling the indicator matrix
$\boldsymbol{\delta}_{k}$ corresponds to navigating through an extremely high
dimensional model space. This is even more challenging, if the sparse factor
model is overfitting. In the present paper, an MCMC scheme summarized in
Algorithm 1 is developed where several steps are designed specifically for
sparse Bayesian factor models with UGLT structure where the factor dimension
is unknown. Updating $\boldsymbol{\delta}_{k}$ for sparse exploratory Bayesian
factor analysis without imposing identification constraints on
$\boldsymbol{\delta}_{k}$ is fairly straightforward, see e.g. Carvalho et al.
(2008) and Kaufmann and Schuhmacher (2019), among many others. A more refined
approach is implemented in the present paper under sparse UGLT structures
which will allow us to address econometric identification in a post-processing
manner in Section 5.

######  Algorithm 1 (MCMC estimation for sparse Bayesian factor models with
unordered GLT structures).

Choose initial values for $\boldsymbol{\delta}_{r}$, ${\mathbf{f}}_{r}$ and
$r_{\footnotesize sp}$.555See Appendix H for details. Iterate $M$ times
through the following steps and discard the first $M_{0}$ draws as burn-in:

  * (CFA)

Update all unknowns in the CFA model (2.8) corresponding to
$\boldsymbol{\delta}_{r}$:

  * (H)

Update any unknown hyperparameters in the ESP prior (3.2) without conditioning
on the slab probabilities $\tau_{1},\ldots,\tau_{r}$. For $j=1,\ldots,r$,
sample
$\tau_{j}|l_{j},d_{j}\sim\mathcal{B}\left(a_{k}+d_{j}-1,b_{k}+m-l_{j}-d_{j}+1\right)$,
where $d_{j}=\sum_{i=1}^{m}{\delta_{ij}}$.

  * (D)

Loop over all columns of the sparsity matrix $\boldsymbol{\delta}_{r}$ in a
random order and perform variable selection by updating the $j$th column
$\boldsymbol{\delta}^{r}_{\cdot,j}$ below the pivot, conditional on the
remaining columns $\boldsymbol{\delta}^{r}_{\cdot,-j}$, the factors
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$ and the
hyperparameter $\tau_{j}$ without conditioning on the model parameters
$\boldsymbol{\beta}_{r}$ and $\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$:

    * (a)

Sample all indicators $\delta_{ij}$ below the pivot $l_{j}$ from
$p(\delta_{ij}|l_{j},\boldsymbol{\delta}^{r}_{\cdot,-j},{\mathbf{f}}_{r},\tau_{j},{\mathbf{y}})$.

    * (b)

If column $\boldsymbol{\delta}^{r}_{\cdot,j}$ is spurious after this update,
increase $r_{\footnotesize sp}$ by one. Remove the $j$th column from
$\boldsymbol{\delta}_{r}$, the factors $f_{jt},t=1,\ldots,T$ from
${\mathbf{f}}_{r}$ and $\tau_{j}$ from $\boldsymbol{\tau}_{r}$ to define,
respectively, $\boldsymbol{\delta}_{r-1}$, ${\mathbf{f}}_{r-1}$ and
$\boldsymbol{\tau}_{r-1}$ and decrease $r$ by one.

  * (L)

Loop over all columns of $\boldsymbol{\delta}_{r}$ in a random order and
update the pivots without conditioning on $\boldsymbol{\beta}_{r}$,
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$ and the slab probabilities
$\boldsymbol{\tau}_{r}$:

    * (a)

Sample a new pivot row $l_{j}$ in each column $j$ from
$p(l_{j}|\boldsymbol{\delta}^{r}_{\cdot,-j},{\mathbf{f}}_{r},{\mathbf{y}})$.

    * (b)

If column $\boldsymbol{\delta}^{r}_{\cdot,j}$ is spurious after this update,
increase $r_{\footnotesize sp}$ by one. Remove the $j$th column from
$\boldsymbol{\delta}_{r}$, the factors $f_{jt},t=1,\ldots,T$ from
${\mathbf{f}}_{r}$ and $\tau_{j}$ from $\boldsymbol{\tau}_{r}$ to define,
respectively, $\boldsymbol{\delta}_{r-1}$, ${\mathbf{f}}_{r-1}$ and
$\boldsymbol{\tau}_{r-1}$ and decrease $r$ by one.

  * (P)

Sample the model parameters $\boldsymbol{\beta}_{r}$ and
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$ jointly conditional on the sparsity
matrix $\boldsymbol{\delta}_{r}$ and the factors
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$ from
$p(\boldsymbol{\beta}_{r},\sigma^{2}_{1},\ldots,\sigma^{2}_{m}|\boldsymbol{\delta}_{r},{\mathbf{f}}_{r},{\mathbf{y}})$.

  * (F)

Sample the latent factors
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$
conditional on the model parameters $\boldsymbol{\beta}_{r}$ and
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$ from
$p({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T}|\boldsymbol{\beta}_{r},\sigma^{2}_{1},\ldots,\sigma^{2}_{m},{\mathbf{y}})$.

  * (S)

For hierarchical Gaussian priors, update all unknown scaling factors, i.e. the
global shrinkage parameter $\kappa$, the column-specific shrinkage parameters
$\theta_{j}$, and the local shrinkage parameters $\omega_{ij}$.

  * (A)

Perform a boosting step to enhance mixing.

  * (EFA)

Move from the current CFA model to an EFA model (1.3) with $r_{\footnotesize
sp}$ spurious columns:

    * (R)

Use Algorithm 2 to change $r_{\footnotesize sp}$, while holding the number of
active factors $r$ fixed. Loop over all columns of the spurious factors
$\boldsymbol{\delta}_{\Xi}$ and try to turn them into active factors.

Move from the current EFA model back to the CFA model (2.8) and preserve the
updated number $r_{\footnotesize sp}$ of spurious columns.

Algorithm 1 consists of two main blocks. Block (CFA) operates in the
confirmatory factor analysis model (2.8) corresponding to
$\boldsymbol{\delta}_{r}$. Due to the prior specification in Section 3, the
number $r_{\footnotesize sp}$ of spurious columns is a sufficient statistic
for the remaining columns in $\boldsymbol{\delta}_{k}$ and no further
information, such as the position of the pivots ${\mathbf{l}}_{\Xi}$ in
$\boldsymbol{\delta}_{\Xi}$ or the spurious loading matrix $\boldsymbol{\Xi}$,
is needed to update the parameters in the CFA model, namely the inclusion
probabilities $\boldsymbol{\tau}_{r}$, the sparsity matrix
$\boldsymbol{\delta}_{r}$, the loading matrix $\boldsymbol{\beta}_{r}$, the
idiosyncratic variances $\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$ in
${\mathbf{\Sigma}}_{r}$, the latent factors
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$, and all
unknown shrinkage factors.

In Block (EFA), the sampler moves from the current CFA model to an EFA model
with $r_{\footnotesize sp}$ spurious columns. In Step (R), dimension changing
moves are performed in the much larger space underlying the EFA model.
Exploiting the results of Section 2.2, spurious factors are added and deleted
with the help of a reversible jump MCMC (RJMCMC) step described in more detail
in Section 4.1. The sampler finally returns to a CFA model with a potentially
larger number of active factors $r$.

To ensure that the loading matrix in the CFA model exhibits a UGLT structure,
Step (L) performs MH steps that navigate through the space of all admissible
$\boldsymbol{\delta}_{r}$ where the pivots
${\mathbf{l}}_{r}=(l_{1},\ldots,l_{r})$ lie in different rows, see Section
4.2. These steps are performed marginalized w.r.t.
$\boldsymbol{\tau}_{r}=(\tau_{1},\ldots,\tau_{r})$. Given the pivots
${\mathbf{l}}_{r}$, the hyperparameters $a_{k}$ and $b_{k}$ in the ESP prior
(3.2) are updated in Step (H) using an MH step without conditioning on the
slab probabilities $\tau_{1},\ldots,\tau_{k}$. Details are provided in
Appendix F. The posterior $\tau_{j}|l_{j},d_{j}$ is updated by combining the
likelihood (3.10) with the prior
$\tau_{j}\sim\mathcal{B}\left(a_{k},b_{k}\right)$ for all columns $j$. In Step
(D), variable selection is performed in each column $j$ for all indicators
$\delta_{ij}$ below the pivot element $l_{j}$. This step potentially turns an
active factor into a spurious one and in this way decreases the number of
active factors $r$, while increasing $r_{\footnotesize sp}$. All moves in Step
(D) are implemented conditionally on $\tau_{j}$ (and all shrinkage parameters
for hierarchical Gaussian priors), as this allows efficient multimove sampling
of all indicators $\\{\delta_{ij},i\in\\{l_{j}+1,\ldots,m\\}\\}$, using
Algorithm 4 in Appendix B.2.

The remaining steps are quite standard in Bayesian factor analysis, see Geweke
and Singleton (1980) and Lopes and West (2004). In Step (F), the conditional
joint posterior
$p({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T}|\boldsymbol{\beta}_{r},\sigma^{2}_{1},\ldots,\sigma^{2}_{m},{\mathbf{y}})$
factors into $T$ independent normal distributions given by:

$\displaystyle{\mathbf{f}}^{r}_{t}|{\mathbf{y}}_{t},\boldsymbol{\beta}_{r},{\mathbf{\Sigma}}_{r}\sim\mathcal{N}_{r}\left(({{\mathbf{I}}}_{r}+\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r})^{-1}\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}{\mathbf{y}}_{t},({{\mathbf{I}}}_{r}+\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r})^{-1}\right).$
(4.1)

In Step (P), we use an efficient algorithm for multi-move sampling of all
unknown model parameters $\boldsymbol{\beta}_{r}$, and
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$, see Appendix A.2. For the hierarchical
Gaussian priors (3.14) and (3.15), all unknown shrinkage parameters are
updated in Step (S), see Appendix C. Finally, the boosting Step (A) is added
to improve the mixing of the MCMC scheme, see Section 4.3 and Appendix G.

###  4.1 Split and merge moves for overfitting models

Step (R) in Algorithm 1 is based on moving from the CFA model (2.8) to an EFA
model with $r_{\footnotesize sp}$ spurious factors in
$\boldsymbol{\delta}_{k}$. As discussed in Section 2.2, spurious columns
$\boldsymbol{\Xi}$ in an EFA model can be substituted by zero columns without
changing the likelihood function, by adding $\boldsymbol{\Xi}$ to the
covariance matrix ${\mathbf{\Sigma}}_{r}$ of the idiosyncratic errors in the
CFA model. On the other hand, for any row $l_{\footnotesize\rm sp}$ that is
not a pivot row in the sparsity matrix $\boldsymbol{\delta}_{r}$ of the CFA
model, a fraction of the idiosyncratic variance
$\sigma^{2}_{l_{\footnotesize\rm sp}}$ can be used to turn a zero column in
$\boldsymbol{\beta}_{k}$ into a spurious column $\boldsymbol{\Xi}$ with a non-
zero loading $\Xi_{l_{\footnotesize\rm sp}}$ without changing the likelihood
function either. This is the cornerstone of the dimension changing procedure
in Step (R), outlined in detail in Algorithm 2.

######  Algorithm 2 (Dimension changing move in an EFA model).

Step (R) in Algorithm 1 is implemented in the following way:

  * (R-S)

Perform an RJMCMC step to change the number $r_{\footnotesize sp}$ of spurious
columns through a split move on a zero column or a merge move on a spurious
column in $\boldsymbol{\delta}_{k}$.

  * (R-L)

Given $r_{\footnotesize sp}$, sample the pivot rows
${\mathbf{l}}_{\Xi}|{\mathbf{l}}_{r}$ of all $r_{\footnotesize sp}$ spurious
columns sequentially from the set ${\cal L}({\mathbf{l}}_{r})$, where
${\mathbf{l}}_{r}$ are the pivot rows of the active factors
$\boldsymbol{\delta}_{r}$. Order the spurious columns such that
$l_{{\Xi},1}<\ldots<l_{{\Xi},r_{\footnotesize sp}}$.

  * (R-F)

Loop over all spurious columns $j_{\footnotesize\rm sp}$ and sample the
factors ${\mathbf{f}}_{j_{\footnotesize\rm sp}}=(f_{j_{\footnotesize\rm
sp},1},\ldots,f_{j_{\footnotesize\rm sp},T})$ independently for all
$t=1,\ldots,T$ from $f_{j_{\footnotesize\rm
sp},t}|{\mathbf{f}}^{r}_{t},\boldsymbol{\beta}_{r},\sigma^{2}_{l_{\footnotesize\rm
sp}},y_{l_{\footnotesize\rm sp},t}\sim\mathcal{N}\left(E_{j_{\footnotesize\rm
sp},t},V_{j_{\footnotesize\rm sp}}\right)$, where $U_{j_{\footnotesize\rm
sp}}$ is a draw from a uniform distribution on [-1,1] and

$\displaystyle V_{j_{\footnotesize\rm sp}}=1-U_{j_{\footnotesize\rm
sp}}^{2},\quad\displaystyle E_{j_{\footnotesize\rm
sp},t}=U_{j_{\footnotesize\rm sp}}\frac{y_{l_{\footnotesize\rm
sp},t}-\boldsymbol{\beta}_{r,l_{\footnotesize\rm
sp},\cdot}{\mathbf{f}}^{r}_{t}}{\sqrt{\sigma}^{2}_{l_{\footnotesize\rm sp}}}.$
(4.2)

  * (R-H)

Sample $\tau_{j_{\footnotesize\rm sp}}|l_{\footnotesize\rm
sp}\sim\mathcal{B}\left(a_{k},b_{k}+m-l_{\footnotesize\rm sp}\right)$ for all
spurious columns $j_{\footnotesize\rm sp}$.

  * (R-D)

Loop over all spurious columns from the last (with the largest pivot) to the
first (with the smallest pivot) and try to turn spurious columns into active
ones:

    * (R-Da)

sample all indicators $\delta_{i,j_{\footnotesize\rm sp}}$ with $i\in
I_{j_{\footnotesize\rm sp}}=\\{l_{\footnotesize\rm sp}+1,\ldots,m\\}$ below
the pivot $l_{\footnotesize\rm sp}$ conditional on $\tau_{j_{\footnotesize\rm
sp}}$, $\boldsymbol{\delta}_{r}$, ${\mathbf{f}}_{r}$ and the spurious factors
${\mathbf{f}}_{j_{\footnotesize\rm sp}}$ without conditioning on
$\boldsymbol{\beta}_{r}$, $\boldsymbol{\Xi}$ and
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$.

    * (R-Db)

If column $j_{\footnotesize\rm sp}$ remains spurious, $r$ and
$r_{\footnotesize sp}$ are unchanged and the EFA model is integrated over
factor $j_{\footnotesize\rm sp}$ by removing
$\boldsymbol{\delta}_{\cdot,j_{\footnotesize\rm sp}}$ from
$\boldsymbol{\delta}_{\Xi}$ and ${\mathbf{f}}_{j_{\footnotesize\rm sp}}$ from
${\mathbf{f}}_{\Xi}$. Otherwise, decrease $r_{\footnotesize sp}$ by 1,
increase $r$ by 1, add $\boldsymbol{\delta}_{\cdot,j_{\footnotesize\rm sp}}$
to $\boldsymbol{\delta}_{r}$ and ${\mathbf{f}}_{j_{\footnotesize\rm sp}}$ to
${\mathbf{f}}_{r}$.

Step (R-S) in Algorithm 2 changes $r_{\footnotesize sp}$ by adding and
deleting spurious columns in $\boldsymbol{\delta}_{k}$. For a given
$r_{\footnotesize sp}$, the conditional prior
$p(\boldsymbol{\delta}_{k},\boldsymbol{\beta}_{k},{\mathbf{\Sigma}}_{k}|r_{\footnotesize
sp})$ is invariant to the specific choice of $l_{\footnotesize\rm sp}$ and
$\boldsymbol{\Xi}$. However, the prior odds that a zero column in
$\boldsymbol{\delta}_{k}$ can be turned into an additional spurious column
depends both on $r$ and $r_{\footnotesize sp}$ (see Appendix D for a proof):

$\displaystyle O^{\mbox{\rm\tiny sp}}(r,r_{\footnotesize sp})=\frac{a_{k}(m-r-
r_{\footnotesize sp})}{b_{k}+m-r-r_{\footnotesize sp}-1}.$ (4.3)

Hence, simply adding or deleting spurious columns would lead to an invalid
MCMC procedure and an RJMCMC step that incorporates $O^{\mbox{\rm\tiny
sp}}(r,r_{\footnotesize sp})$ is performed in Step (R-S). As opposed to other
applications of RJMCMC, the acceptance rate is extremely easy to compute, see
(4.4) and (4.5).

At each sweep of the sampler, a split or a merge move is performed with,
respectively, probability $p_{\footnotesize split}(r,r_{\footnotesize sp})$ or
$p_{\footnotesize merge}(r,r_{\footnotesize sp})$. In a split move (which
requires $r<k$) a randomly chosen zero column in $\boldsymbol{\delta}_{k}$ is
turned into a spurious column. The corresponding proposal density reads
$q_{\mbox{\rm\footnotesize split}}(\boldsymbol{\delta}_{k}^{\mbox{\rm\tiny
new}}|\boldsymbol{\delta}_{k})=p_{\footnotesize split}(r,r_{\footnotesize
sp})/(k-r-r_{\footnotesize sp})$. The merge move (which requires
$r_{\footnotesize sp}>0$) is obtained by reversing the split move and turns
one of the spurious columns into a zero column. The corresponding proposal
density reads $q_{\mbox{\rm\footnotesize
merge}}(\boldsymbol{\delta}_{k}^{\mbox{\rm\tiny
new}}|\boldsymbol{\delta}_{k})=p_{\footnotesize merge}(r,r_{\footnotesize
sp})/r_{\footnotesize sp}$. A symmetric proposal is selected for all $0\leq
r_{\footnotesize sp}<k-r$ such that $p_{\footnotesize
split}(r,r_{\footnotesize sp})=p_{\footnotesize merge}(r,r_{\footnotesize
sp}+1)=p_{r_{\footnotesize sp}}$, where $p_{r_{\footnotesize sp}}\leq 0.5$ is
a tuning parameter, while $p_{\footnotesize merge}(r,r_{\footnotesize sp})=0$
for $r_{\footnotesize sp}=0$ and $p_{\footnotesize split}(r,r_{\footnotesize
sp})=0$ for $r_{\footnotesize sp}=k-r$. We found it useful to choose a fixed
probability $p_{r_{\footnotesize sp}}=p_{s}$, although other choices are
possible. A split move is accepted with probability
$\min(1,A_{\mbox{\rm\footnotesize split}}(r,r_{\footnotesize sp}))$, where:

$\displaystyle A_{\mbox{\rm\footnotesize split}}(r,r_{\footnotesize
sp})=\frac{q_{\mbox{\rm\footnotesize
merge}}(\boldsymbol{\delta}_{k}|\boldsymbol{\delta}_{k}^{\mbox{\rm\tiny
new}})}{q_{\mbox{\rm\footnotesize
split}}(\boldsymbol{\delta}_{k}^{\mbox{\rm\tiny
new}}|\boldsymbol{\delta}_{k})}O^{\mbox{\rm\tiny sp}}(r,r_{\footnotesize
sp})=\frac{a_{k}(m-r-r_{\footnotesize sp})(k-r-r_{\footnotesize
sp})}{(r_{\footnotesize sp}+1)(b_{k}+m-r-r_{\footnotesize sp}-1)},$ (4.4)

whereas a merge move is accepted with probability
$\min(1,A_{\mbox{\rm\footnotesize merge}}(r,r_{\footnotesize sp}))$, where

$\displaystyle A_{\mbox{\rm\footnotesize merge}}(r,r_{\footnotesize
sp})=\frac{1}{A_{\mbox{\rm\footnotesize split}}(r,r_{\footnotesize
sp}-1)}=\frac{r_{\footnotesize sp}(b_{k}+m-r-r_{\footnotesize sp})}{a_{k}(m-r-
r_{\footnotesize sp}+1)(k-r-r_{\footnotesize sp}+1)}.$ (4.5)

Very conveniently, $A_{\mbox{\rm\footnotesize split}}(r,r_{\footnotesize sp})$
and $A_{\mbox{\rm\footnotesize merge}}(r,r_{\footnotesize sp})$ are
independent of the pivots ${\mathbf{l}}_{\Xi}$ in the spurious columns. Note
that there is a dynamic feature underlying this RJMCMC algorithm, with
acceptance depending on the number of spurious columns $r_{\footnotesize sp}$.
For $b_{k}=1$, for instance, $A_{\mbox{\rm\footnotesize
split}}(r,r_{\footnotesize sp})$ is monotonically decreasing and
$A_{\mbox{\rm\footnotesize merge}}(r,r_{\footnotesize sp})$ is monotonically
increasing in $r_{\footnotesize sp}$.

Once $r_{\footnotesize sp}$ has been updated, Step (R-L) is trying to turn
each spurious column into an active one. Since the likelihood is non-
informative about spurious columns, pivots $l_{\footnotesize\rm sp}$ are
sampled from the prior ${\mathbf{l}}_{\Xi}|{\mathbf{l}}_{r}$, while the
spurious factor loadings $\Xi_{l_{\footnotesize\rm sp}}$ are sampled from the
prior (3.12). Given $l_{\footnotesize\rm sp}$, the idiosyncratic variance
$\sigma^{2}_{l_{\footnotesize\rm sp}}$ in the CFA model is split, with the
help of a random variable $U_{j_{\footnotesize\rm
sp}}\sim\mathcal{U}\left[-1,1\right]$, between $\Xi_{l_{\footnotesize\rm sp}}$
and an updated idiosyncratic variance $(\sigma^{2}_{l_{\footnotesize\rm
sp}})^{\mbox{\rm\tiny new}}$. More specifically:

$\displaystyle\Xi_{l_{\footnotesize\rm sp}}=U_{j_{\footnotesize\rm
sp}}\sqrt{\sigma^{2}_{l_{\footnotesize\rm
sp}}},\qquad(\sigma^{2}_{l_{\footnotesize\rm sp}})^{\mbox{\rm\tiny
new}}=(1-U_{j_{\footnotesize\rm sp}}^{2})\sigma^{2}_{l_{\footnotesize\rm
sp}}.$ (4.6)

Given $\Xi_{l_{\footnotesize\rm sp}}$ and $(\sigma^{2}_{l_{\footnotesize\rm
sp}})^{\mbox{\rm\tiny new}}$, factors $f_{j_{\footnotesize\rm sp},t}$ are
proposed in Step (R-F) independently for $t=1,\ldots,T$, from the conditional
density $p(f_{j_{\footnotesize\rm
sp},t}|{\mathbf{f}}^{r}_{t},\boldsymbol{\beta}_{r},\sigma^{2}_{l_{\footnotesize\rm
sp}},y_{l_{\footnotesize\rm sp},t})$ given in (4.2), see Appendix D for a
proof. The slab probabilities $\tau_{j_{\footnotesize\rm sp}}$ are sampled in
Step (R-H) as in Algorithm 1, Step (H) using that $d_{j_{\footnotesize\rm
sp}}=1$. Finally, in Step (R-Da), variable selection is performed in each
spurious column on all indicators below $l_{\footnotesize\rm sp}$ as in Step
(D) of Algorithm 1, conditional on $f_{j_{\footnotesize\rm sp},t}$, but
marginalizing w.r.t. the idiosyncratic variances and the factor loading
matrices $\boldsymbol{\beta}_{r}$ and $\boldsymbol{\Xi}$. More details about
this step are provided in Appendix D. Any spurious column that is active after
this step is integrated into the CFA model in Step (R-Db), increasing in this
way the number of active columns $r$.

Figure 1: MCMC moves to change the leading indices of an unordered GLT
structure; from left to right: shifting the leading index, adding a new
leading index, deleting a leading index and switching the leading elements

###  4.2 Special MCMC moves for unordered GLT structures

Step (L) in Algorithm 1 implements moves that explicitly change the position
of the pivots ${\mathbf{l}}_{r}=(l_{1},\ldots,l_{r})$ in the $r$ columns of
the UGLT indicator matrix $\boldsymbol{\delta}_{r}$. We scan all columns of
$\boldsymbol{\delta}_{r}$ in a random order and propose to change
$l_{j}|{\mathbf{l}}_{r,-j}$ given the pivots ${\mathbf{l}}_{r,-j}$ in the
other columns. To this aim, we use one of four MH moves, namely shifting the
pivot, adding a new pivot, deleting a pivot and switching the pivots (and
additional indicators) between column $j$ and a randomly selected column
$j^{\prime}$; see Figure 1 for illustration. All moves are performed
marginalized w.r.t. $\boldsymbol{\tau}_{r}$. Changing the pivot from $l_{j}$
to $l_{j}^{\mbox{\rm\tiny new}}$ changes the number of unconstrained
indicators, whereas the prior ratio

$\frac{p(l^{\mbox{\rm\tiny
new}}_{j}|{\mathbf{l}}_{r,-j})}{p(l_{j}|{\mathbf{l}}_{r,-j})}=1,$

since $l_{j}$ has a uniform prior over ${\cal L}({\mathbf{l}}_{r,-j})$, the
set of admissible pivot rows in column $j$. With $d_{j}^{\mbox{\rm\tiny new}}$
being the new number of non-zero elements in column $j$, the prior ratio
$R_{\mbox{\rm\footnotesize move}}$ can be derived from (3.11):

$\displaystyle R_{\mbox{\rm\footnotesize move}}=\frac{\mbox{\rm
Pr}(\boldsymbol{\delta}_{\cdot,j}^{\mbox{\rm\tiny new}}|l_{j}^{\mbox{\rm\tiny
new}})}{\mbox{\rm
Pr}(\boldsymbol{\delta}_{\cdot,j}|l_{j})}=\frac{B(a_{k}+d_{j}^{\mbox{\rm\tiny
new}}-1,b_{k}+m-l_{j}^{\mbox{\rm\tiny new}}-d_{j}^{\mbox{\rm\tiny
new}}+1)}{B(a_{k}+d_{j}-1,b_{k}+m-l_{j}-d_{j}+1)}.$ (4.7)

Further details are provided in Appendix E.

###  4.3 Boosting MCMC

Step (F) and Step (P) in Algorithm 1 perform full conditional Gibbs sampling
for a confirmatory factor model with sparsity matrix
$\boldsymbol{\delta}_{r}$, by sampling the factors ${\mathbf{f}}^{r}_{t}$
conditionally on the loading matrix $\boldsymbol{\beta}_{r}$ and the
idiosyncratic covariance matrix ${\mathbf{\Sigma}}_{r}$ and sampling
$\boldsymbol{\beta}_{r}$ and ${\mathbf{\Sigma}}_{r}$ conditionally on
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$.
Depending on the signal-to-noise ratio of the latent variable representation,
such full conditional Gibbs sampling tends to be poorly mixing. In a CFA
model, where
${\mathbf{f}}^{r}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},{{\mathbf{I}}}_{r}\right)$,
the information in the data (the “signal”) can be quantified by the matrix
$\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r}$
in comparison to the identity matrix ${{\mathbf{I}}}_{r}$ (the “noise”) in the
filter for
${\mathbf{f}}^{r}_{t}|{\mathbf{y}}_{t},\boldsymbol{\beta}_{r},{\mathbf{\Sigma}}_{r}$,
see (4.1). In particular for large factor models with many measurements, one
would expect that the data contain ample information to estimate the factors
${\mathbf{f}}^{r}_{t}$. However, this is true only, if the information matrix
$\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r}$
increases with $m$, hence if most of the factor loadings are nonzero. For
sparse factor models many columns with quite a few zero loadings are present,
leading to a low signal-to-noise ratio and ,consequently, to a poorly mixing
Gibbs sampler, as illustrated in the left-hand panel in Figure 2 showing
posterior draws of $\mbox{\rm
tr}(\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r})$
without boosting Step (A) for the exchange data to be discussed in Section
6.1.

|  |
---|---|---
Figure 2: Exchange rate data; fractional prior. Posterior draws of $\mbox{\rm
tr}(\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r})$
without boosting (left-hand side), boosting through ASIS with the largest
loading (in absolute values) in each nonzero column serving as
$\sqrt{\Psi_{j}}$ (middle) and boosting through MDA based on the inverted
Gamma working prior $\Psi_{j}\sim\mathcal{G}^{-1}\left(1.5,1.5\right)$ (right-
hand side).

Hence, boosting steps are essential to obtain efficient MCMC schemes for
sparse factor models. Several papers (Ghosh and Dunson, 2009; Frühwirth-
Schnatter and Lopes, 2010; Conti et al., 2014; Piatek and Papaspiliopoulos,
2018) apply marginal data augmentation (MDA) in the spirit of van Dyk and Meng
(2001); others (Kastner et al., 2017; Frühwirth-Schnatter and Lopes, 2018)
exploit the ancillarity-suffiency interweaving strategy (ASIS) introduced by
Yu and Meng (2011). Some boosting strategies enhances mixing at the cost of
changing the prior of the factor loading matrix $\boldsymbol{\beta}_{r}$
(Ghosh and Dunson, 2009), however, this appears undesirable in a variable
selection context and is avoided by the boosting strategies applied in the
present paper.

Boosting is based on moving from the CFA model (2.8) where
${\mathbf{f}}^{r}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},{{\mathbf{I}}}_{r}\right)$
to an expanded model with a more general prior:

$\displaystyle{\mathbf{y}}_{t}=\tilde{\boldsymbol{\beta}}_{r}\tilde{\mathbf{f}}^{r}_{t}+\boldsymbol{\epsilon}_{t},\quad\boldsymbol{\epsilon}_{t}\sim\mathcal{N}_{m}\left({\mathbf{0}},{\mathbf{\Sigma}}_{r}\right),\qquad\tilde{\mathbf{f}}^{r}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},\boldsymbol{\Psi}\right),$
(4.8)

where $\boldsymbol{\Psi}=\mbox{\rm Diag}\,(\Psi_{1},\ldots,\Psi_{r})$ is a
diagonal matrix. The relation between the two systems is given by the
transformations
$\tilde{\mathbf{f}}^{r}_{t}=(\boldsymbol{\Psi})^{1/2}{\mathbf{f}}^{r}_{t}$ and
$\tilde{\boldsymbol{\beta}}_{r}=\boldsymbol{\beta}_{r}(\boldsymbol{\Psi})^{-1/2}$,
where the nonzero elements in $\tilde{\boldsymbol{\beta}}_{r}$ have the same
position as the nonzero elements in $\boldsymbol{\beta}_{r}$ and the sparsity
matrix $\boldsymbol{\delta}_{r}$ is not affected by the transformation. The
main difference between MDA and ASIS lies in the choice of
$\boldsymbol{\Psi}$. While $\Psi_{j}$ is sampled from a working prior for MDA,
$\Psi_{j}$ is chosen in a deterministic fashion for ASIS, see Appendix G for
details. For illustration, Figure 2 shows considerable efficiency gain in the
posterior draws of $\mbox{\rm
tr}(\boldsymbol{\beta}_{r}^{{}^{\prime}}{\mathbf{\Sigma}}_{r}^{-1}\boldsymbol{\beta}_{r})$,
when a boosting strategy such as ASIS (middle panel) or MDA (right-hand panel)
is applied.

For hierarchical Gaussian priors with column specific shrinkage parameters,
such as (3.14) and (3.15), we found it more useful to apply column boosting
and interweave $\theta_{j}$ into the state equation by choosing
$\Psi_{j}=\theta_{j}$. Column boosting and an additional boosting step which
interweaves the global shrinkage parameter $\kappa$ into the prior of the
shrinkage parameter $\theta_{j}$ is discussed in detail in Appendix G.

##  5 Post-processing posterior draws

Algorithm 1 delivers posterior draws
$(\boldsymbol{\delta}_{r},\boldsymbol{\beta}_{r},{\mathbf{\Sigma}}_{r})$ in a
CFA model with a varying number $r$ of active columns. Instead of sampling
these parameters without any constraints, our sampler imposes the (mild)
constraint that the pivots (the first non-zero loading in each column) lie in
different rows, ensuring that posterior draws of the loading matrix exhibit a
UGLT structure. While these draws are not identified in the rigorous sense
discussed in Section 2.2, the UGLT structure allows identification during
post-processing.

We use the 3579 counting rule and the algorithm of Hosszejni and Frühwirth-
Schnatter (2022) to check for each draw $\boldsymbol{\delta}_{r}$, if the
variance decomposition is unique, and remove all posterior draws that are not
variance identified. In our experience, the fraction $M_{V}$ of variance
identified draws is relatively high for reasonably chosen priors. A low
fraction of variance identified draws can be the result of poorly chosen
shrinkage priors and is a hint to consider another prior family. The
fractional prior on the loadings $\beta_{ij}$, for instance, tends to deliver
a low fraction of variance identified draws in high-dimensional factor models,
see Section 6.

Many quantities can be inferred from the variance identified posteriors draws
with varying factor dimension $r$, such as the marginal covariance matrix
${\mathbf{\Omega}}=\boldsymbol{\beta}_{r}\boldsymbol{\beta}_{r}^{T}+{\mathbf{\Sigma}}_{r}$
and the idiosyncratic variances $\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$
appearing in the diagonal of ${\mathbf{\Sigma}}_{r}$. In addition, overall
sparsity in terms of the number $d=\sum_{j=1}^{r}\sum_{i=1}^{m}\delta_{ij}$ of
nonzero elements in $\boldsymbol{\delta}_{r}$, can be evaluated. Posterior
draws of $d$ are particularly useful to check convergence and assessing
efficiency of the MCMC sampler, as $d$ captures the ability of the sampler to
move across (variance identified) UGLT factor models of different dimensions.
Functionals of ${\mathbf{\Omega}}$ and ${\mathbf{\Sigma}}_{r}$, such as the
trace, the log determinant, and
${\mathbf{1}}_{m}^{{}^{\prime}}{\mathbf{\Omega}}^{-1}{\mathbf{1}}_{m}$ are
further useful means to assess MCMC convergence. Furthermore, for each
variable $y_{it}$, inference with respect to the proportion of the variance
explained by the common factors (also known as communalities $R^{2}_{i}$) is
possible:

$\displaystyle R^{2}_{i}=\sum_{j=1}^{r}R_{ij}^{2},\qquad
R_{ij}^{2}=\frac{\beta_{ij}^{2}}{\sum_{l=1}^{r}\beta_{il}^{2}+\sigma^{2}_{i}}.$
(5.1)

Most importantly, variance identified draws are instrumental for estimating
the number of factors from an overfitting factor model with cumulative
shrinkage process priors. As shown by Frühwirth-Schnatter et al. (2022,
Theorem 7), the number $r$ of nonzero columns in $\boldsymbol{\delta}_{r}$ is
equal to the factor dimension, if the variance decomposition is unique for
$r$. In this way, posterior draws of $r$ can be used for variance identified
factor loading matrices to draw posterior inference on the factor dimension
$r$. The mode $\tilde{r}$ of the posterior distribution $p(r|{\mathbf{y}})$,
e.g., can be used to estimate the factor dimension.

Due to the UGLT structure imposed on $\boldsymbol{\beta}_{r}$, rotational
invariance reduces for variance identified draws to sign and column switching
(Frühwirth-Schnatter et al., 2022, Theorem 1) and $\boldsymbol{\Lambda}$ and
$\boldsymbol{\delta}^{\Lambda}$ are easily recovered. First, the columns of
$\boldsymbol{\delta}_{r}$ and $\boldsymbol{\beta}_{r}$ are ordered such that
the pivots ${\mathbf{l}}_{r}=(l_{1},\ldots,l_{r})$ obey $l_{1}<\ldots<l_{r}$;
i.e.
$\boldsymbol{\delta}^{\Lambda}=\boldsymbol{\delta}_{r}{\mathbf{P}}_{\rho}$.
Then, the sign of all columns in $\boldsymbol{\beta}_{r}{\mathbf{P}}_{\rho}$
is switched if the leading element is negative; i.e.
$\boldsymbol{\Lambda}=\boldsymbol{\beta}_{r}{\mathbf{P}}_{\rho}{\mathbf{P}}_{\pm}$.
In addition, the factors ${\mathbf{f}}^{r}_{t}$ are reordered through
${\mathbf{P}}_{\pm}^{{}^{\prime}}{\mathbf{P}}_{\rho}^{{}^{\prime}}{\mathbf{f}}^{r}_{t}$
for $t=1,\ldots,T$. Finally, ${\mathbf{P}}_{\rho}$ is used to reorder the
draws of $\boldsymbol{\tau}_{r}=(\tau_{1},\ldots,\tau_{r})$ and any local and
column-specific shrinkage parameters in a hierarchical Gaussian shrinkage
prior.

The posterior draws of $\boldsymbol{\delta}^{\Lambda}$ are exploited in
various ways. The highest probability model (HPM), i.e. the indicator matrix
$\boldsymbol{\delta}^{\Lambda}_{H}$ visited most often, its frequency $p_{H}$
(an estimator of the posterior probability of the HPM), its factor dimension
$r_{H}$, its model size $d_{H}$, and its sequence of pivots ${\mathbf{l}}^{H}$
are of interest. For each draw of $\boldsymbol{\delta}^{\Lambda}$, the ordered
pivots $l_{1}<\ldots<l_{r}$ are draws from the marginal posterior distribution
$p(l_{1},\ldots,l_{r}|{\mathbf{y}})$, allowing additional posterior inference
w.r.t. ${\mathbf{l}}_{r}$. First of all, the posterior probability that a
specific feature (row) serves as a pivot, i.e. $\mbox{\rm
Pr}(i\in{\mathbf{l}}_{r}|{\mathbf{y}})$ for $i=1,\ldots,m$ is of interest.
Second, the sequence of pivots
${\mathbf{l}}^{\star}=(l_{1}^{\star},\ldots,l_{r^{\star}}^{\star})$ visited
most often is determined together with its frequency $p_{L}$ which reflects
posterior uncertainty with respect to choosing the pivots. The number
$r^{\star}$ of elements in ${\mathbf{l}}^{\star}$ provides yet another
estimator of the number of factors. If the frequencies $p_{H}$ and $p_{L}$ are
small, the estimator $r^{\star}$ and $r_{H}$ might not coincide with the
posterior mode $\tilde{r}$ and ${\mathbf{l}}^{H}$ will different from
${\mathbf{l}}^{\star}$.

To estimate the factor loading matrix for a chosen number of factors $r$,
Bayesian model averaging is performed conditionally on an estimator
$\hat{{\mathbf{l}}}_{r}$ of the pivots. Averaging over all posterior draws
with these specific pivots avoids column switching and provides an estimate of
the factor loading matrix $\boldsymbol{\Lambda}$ and the marginal inclusion
probabilities $\mbox{\rm
Pr}(\delta^{\Lambda}_{ij}=1|{\mathbf{y}},\hat{{\mathbf{l}}}_{r})$ for all
elements of the corresponding sparsity matrix. The median probability model
(MPM) $\boldsymbol{\delta}^{\Lambda}_{M}$ is obtained by setting each
indicator to one whenever $\mbox{\rm
Pr}(\delta_{ij}^{\Lambda}=1|{\mathbf{y}},\hat{{\mathbf{l}}}_{r})\geq 0.5$. As
for any Bayesian approach, the reliability of any of these estimators might
depend on how informative the data are. In settings where the data are not
very informative, the chosen degree of sparsity in the prior might exhibit
considerable influence on the result.

##  6 Applications

###  6.1 Sparse factor analysis for exchange rate data

As a first application, we analyze log returns from $m=22$ exchange rates with
respect to the Euro, observed for $T=96$ months.666The data was obtained from
the European Central Bank’s Statistical Data Warehouse and ranges from January
3, 2000 to December 3, 2007. It contains the 22 exchange rates listed in Table
I.1 in Appendix I from which we derived monthly returns based on the first
trading day in a month. The data are demeaned and standardized.A sparse EFA
with UGLT structure is fitted with $k=10$ equals the upper bound (2.2).
Regarding the prior on the sparsity matrix $\boldsymbol{\delta}_{k}$, we
consider the 1PB prior (3.3) and the 2PB prior (3.4), where
$\alpha\sim\mathcal{G}\left(6,3\right)$ and
$\gamma\sim\mathcal{G}\left(6,6\right)$ for the 2PB prior. We combine the
fractional prior (3.17) as slab distribution for the non-zero loadings in
$\boldsymbol{\beta}_{r}$ with the prior (3.20) on the idiosyncratic variances
$\sigma^{2}_{1},\ldots,\sigma^{2}_{m}$, where $c_{0}=2.5$ and
$\widehat{{\mathbf{\Omega}}^{-1}}$ is estimated from (3.21) with $\nu_{o}=3$
and ${\mathbf{S}}_{o}={{\mathbf{I}}}_{m}$.

Figure 3: Exchange rate data; posterior draws of the factor dimension $r$
(left-hand side) and model size $d$ (right-hand side) including burn-in,
starting from $r=2$ and $r_{\footnotesize sp}=3$. Table 1: Exchange rate data;
posterior distribution $p(r|{\mathbf{y}})$ of the number of factors and
posterior mean $\mbox{\rm E}(d|{\mathbf{y}})$ of model size (based on
$100p_{V}$ percent variance identified draws).

|  | $p(r|{\mathbf{y}})$ |  |  |
---|---|---|---|---|---
| $\mbox{\rm E}(\alpha|{\mathbf{y}})$ | $\mbox{\rm E}(\gamma|{\mathbf{y}})$ | 0-2 | 3 | 4 | 5 | 6 | 7-10 | $100\cdot p_{V}$ | $\mbox{\rm E}(d|{\mathbf{y}})$
1PB | 2.0 | - | 0 | 0.128 | 0.867 | 0.005 | $\approx 0$ | 0 | 93.3 | 28
2PB | 2.0 | 1.1 | 0 | 0 | 0.988 | 0.012 | 0 | 0 | 95.2 | 28

Note: non-zero probabilities smaller than $10^{-3}$ are indicated by $\approx
0$.

Table 2: Exchange rate data; posterior probability of the event $\mbox{\rm Pr}(q_{i}=0|{\mathbf{y}})$, where $q_{i}$ is the row sum of $\boldsymbol{\delta}_{r}$, for various currencies. | $\mbox{\rm Pr}(q_{i}=0|{\mathbf{y}})$
---|---
Currency | CHF | CZK | MXN | NZD | RON | RUB | remaining
1PB prior | 0.87 | 0.73 | 0.81 | 0.48 | 0.62 | 0.62 | 0
2PB prior | 0.88 | 0.76 | 0.82 | 0.49 | 0.62 | 0.62 | 0

Algorithm 1 is run for $M=40{,}000$ after a burn-in of $20{,}000$ draws. As
discussed in Section 4, this sampler navigates in the space of all unordered
GLT structures with an unknown number of nonzero columns and unknown pivots
without forcing any further constraint. To verify convergence, independent
MCMC chains are started with, respectively, $r=2$ and $r=7$, and
$r_{\footnotesize sp}=3$ spurious columns. The sampler shows good mixing
across models of different dimension, in particular for the 2PB prior, with an
inefficiency factor of roughly 5. For illustration, Figure 3 shows all
posterior draws of $r$ and the model size $d$, including burn-in, for the
first run under the 2PB prior.

Table 3: Exchange rate data; total number of visited models $N_{v}$; frequency $p_{H}$ (in percent), pivots ${\mathbf{l}}_{r}^{H}$ and model size $d_{H}$ of the HPM; pivots ${\mathbf{l}}^{\star}$ visited most often, corresponding frequency $p_{L}$ (in percent) and corresponding number of factors $r^{\star}$; pivots ${\mathbf{l}}_{H}$, model size $d_{H}$ and frequency $p_{H}$ (in percent) of the HPM; model size $d_{M}$ of the MPM. Prior | $N_{v}$ | $100p_{H}$ | ${\mathbf{l}}_{H}$ | $d_{H}$ | ${\mathbf{l}}^{\star}$ | $100p_{L}$ | $r^{\star}$ | $d_{M}$
---|---|---|---|---|---|---|---|---
1PB prior | 16508 | 4.5 | (1,2,5,7) | 26 | (1,2,5,7) | 84.8 | 4 | 26
2PB prior | 11933 | 5.0 | (1,2,5,7) | 26 | (1,2,5,7) | 91.0 | 4 | 26

As outlined in Section 5, we resolve identification during post-processing.
First, we screen for variance identified draws. The fraction $p_{V}$ of
variance identified draws, reported in Table 1, is very high for both EPS
priors. For the variance identified draws, the number $r$ of columns of the
sparsity matrix $\boldsymbol{\delta}_{r}$ in the CFA model is regarded as a
draw of the number $r$ of factors under the specific prior choice. Table 1
reports the posterior distribution $p(r|{\mathbf{y}})$ for both ESP priors.
For the 2PB prior, this posterior is highly concentrated at a four factors.
Also under the 1PB prior, the posterior mode equals 4, but also three factors
receive some posterior probability. The indicator matrix
$\boldsymbol{\delta}_{r}$ is pretty sparse, with an average posterior model
size of 28.

Furthermore, the variance identified draws are used to explore if some
measurements are uncorrelated with the remaining measurements. This is
investigated in Table 2 through the posterior probability $\mbox{\rm
Pr}(q_{i}=0|{\mathbf{y}})$, where $q_{i}$ is the row sum of
$\boldsymbol{\delta}_{r}$. Regardless of the chosen prior, the Swiss franc
(CHF), the Mexican peso (MXN) and the Czech koruna (CZK) have considerable
probability to be uncorrelated with the rest, while the situation is less
clear for the New Zealand dollar (NZD), the Romania fourth leu (RON), and the
Russian ruble (RUB). The remaining currencies are clearly correlated.

Figure 4: Exchange rate data; sparsity matrix $\boldsymbol{\delta}_{4}$
corresponding both to the HPM and the MPM which are identical for the 1PB and
the 2PB prior.

To proceed with identification for all variance identified draws, rotation
indeterminacy is resolved by ordering the pivots (and all column-specific
variables) such that $l_{1}<\ldots<l_{r}$ and an ordered GLT structure is
imposed. Several posterior summaries for the ordered GLT draws are reported in
Table 3. For both ESP priors,
${\mathbf{l}}^{\star}={\mathbf{l}}^{H}=(1,2,5,7)$ is the most likely sequence
of pivots. As a final step, all variance identified, ordered GLT draws where
the pivots coincide with ${\mathbf{l}}^{\star}={\mathbf{l}}^{H}=(1,2,5,7)$ are
used for both ESP priors to identify the marginal inclusion probabilities
$\mbox{\rm Pr}(\delta_{ij}=1|{\mathbf{y}},{\mathbf{l}}^{\star})$, the
corresponding MPM, its model size $d_{M}$ and the factor loading matrix
$\boldsymbol{\Lambda}$ for a 4-factor model. Sign switching in the posterior
draws is resolved by imposing the constraint $\Lambda_{11}>0$,
$\Lambda_{22}>0$, $\Lambda_{53}>0$, and $\Lambda_{74}>0$ on
$\boldsymbol{\Lambda}$. Both ESP priors yield the same MPM which coincides
with the sparsity matrix of both HPMs, see Figure 4 for illustration. Further
results are reported in Appendix I. The resulting model indicates considerable
sparsity, with many factor loadings being shrunk toward zero. Factor 2 is a
common factor among the correlated currencies, while the remaining factors are
three group specific, for the most part dedicated factors.

Figure 5: NYSE data; posterior draws of the total number of non-zero columns
$r+r_{\footnotesize sp}$ (top left), the number of spurious columns
$r_{\footnotesize sp}$ (top right), the extracted number of factors $r$
(bottom left), and the model dimension $d$ (bottom right).

###  6.2 Sparse factor analysis for NYSE stock returns

As a second application, we consider monthly log returns from $m=63$ firms
from the NYSE observed for $T=247$ months from February 1999 till August
2019.777The top 150 companies (as of September 13, 2019) listed on the NYSE
were downloaded from Bloomberg on September 13, 2019. Since many of the
companies entered the NYSE after 1990, we picked February 1999 as a starting
date and, after removing all companies that were founded later, we use monthly
return (determined on the last trading day in each month) for the 103
remaining companies listed on the NYSE for all $T=247$ observation periods
till August 2019. For our case study, we consider the 63 firms belonging to
the following five sectors: basic industries (1-9), non-durable consumer goods
(8-17), energy (18-27), finance (28-45) and health care (46-63). The data are
ordered according to industry and are The data are demeaned and standardized.
An EFA model with UGLT structure is fitted with $k=31$ being equal to the
upper bound (2.2). Regarding the ESP prior on the sparsity matrix
$\boldsymbol{\delta}_{k}$, we consider the 2PB prior (3.4) with
$\alpha\sim\mathcal{G}\left(6,12\right)$ and
$\gamma\sim\mathcal{G}\left(6,6\right)$. As in Section 6.1, we tried to apply
a fractional prior as slab distribution, however the fraction of variance
identified posterior draws was extremely low (less than 1%). Hence, the
hierarchically structured Gaussian shrinkage prior (3.15) is chosen in the
slab. To introduce aggressive shrinkage for the factor loadings, the local
scaling parameters are assumed to follow a triple gamma prior with
$a^{\omega}=c^{\omega}=0.2$, i.e. $\omega_{ij}\sim\mbox{\rm
F}\left(0.4,0.4\right)$. Also the column specific shrinkage parameters follow
a triple gamma prior, $\theta_{j}\sim\mbox{\rm F}\left(5,5\right)$, whereas
the global shrinkage parameter follows an inverse gamma distribution,
$\kappa\sim\mathcal{G}^{-1}\left(10,50\right)$. Regarding the idiosyncratic
variances, we choose the prior
$\sigma^{2}_{i}\sim\mathcal{G}^{-1}\left(2.5,1.5\right)$. The fraction of
variance identified MCMC draws under this prior is roughly 20%.

Table 4: NYSE data; posterior distribution $p(r|{\mathbf{y}})$ of the number
of factors and the mean and the quartiles (in parenthesis) of the posterior
distribution $p(d|{\mathbf{y}})$ of model size ($100p_{V}$ percent variance
identified draws).

| $p(r|{\mathbf{y}})$ |  |
---|---|---|---
$\mbox{\rm E}(\alpha|{\mathbf{y}})$ | $\mbox{\rm E}(\gamma|{\mathbf{y}})$ | 0-12 | 13 | 14 | 15 | 16 | 17 | 18-31 | $p(d|{\mathbf{y}})$
2.5 | 0.8 | 0 | $\approx 0$ | 0.413 | 0.488 | 0.098 | 0.001 | 0 | 203 (210,218)

Note: non-zero probabilities smaller than $10^{-3}$ are indicated by $\approx
0$.

Figure 6: NYSE data; posterior probability $\mbox{\rm
Pr}(i\in{\mathbf{l}}_{r}|{\mathbf{y}})$ that a specific firm serves as pivot;
the red dots indicate the first firm in each sector.

Figure 7: NYSE data; estimated marginal correlation matrix $\mbox{\rm
E}({\mathbf{\Omega}}^{\star}|{\mathbf{y}})$, where
$\Omega^{\star}_{i\ell}=\mbox{\rm Corr}((y_{it}-\Lambda_{i1}f_{1t})(y_{\ell
t}-\Lambda_{\ell 1}f_{1t}))$.

Algorithm 1 was applied to obtain $M=50{,}000$ posterior draws after a burn-in
of $40{,}000$ draws, starting with $r=10$ factors and $r_{\footnotesize sp}=3$
spurious columns. The MCMC scheme shows relatively good mixing, despite the
high dimensionality, as illustrated by Figure 5 showing posterior draws of the
total number of non-zero columns, $r+r_{\footnotesize sp}$, the number of
spurious columns $r_{\footnotesize sp}$, the extracted number of factors $r$,
and the model dimension $d$. As shown in Table 4, the posterior distribution
$p(r|{\mathbf{y}})$ derived from the variance identified draws yields a
posterior mode of $\tilde{r}=15$, but also 14 factors receive considerable
posterior evidence.

Figure 8: NYSE data; factors with pivots equal to (1, 8, 18, 19, 28, 30, 46,
53)

For further inference, an ordered GLT structure is imposed on all variance
identified draws with $r=15$. Since each draw of $\boldsymbol{\Lambda}$ turned
out to have a different shrinkage matrix $\boldsymbol{\delta}^{\Lambda}$,
choosing a unique set of pivots is challenging. In Figure 6, the posterior
probability $\mbox{\rm Pr}(i\in{\mathbf{l}}_{r}|{\mathbf{y}})$ that a specific
firm serves as a pivot is displayed. While the first three factors use the
pivots 1, 2, and 3 (as a PLT structure would do), we find that the remaining
factors clearly exhibit a GLT structure and the first firm listed in a
specific sectors (indicated by a red dot) typically serves as a pivot for new
factor.

Since the pivot of the first factor is equal to 1 for all posterior draws, we
can estimate the first column of the indicator matrix and the corresponding
factor loadings as the average of all posterior draws. This factor is a market
factor that loads on all 63 firms, see also Figure 8. The remaining 14 factors
mainly capture industry specific correlations as well as cross-sectional
correlations between specific firms, see the estimated marginal correlation
matrix ${\mathbf{\Omega}}^{\star}$ that remains after extracting the first
factor in Figure 7. For further illustration, factors with pivots equal to (1,
8, 18, 19, 28, 30, 46, 53) are extracted from all GLT draws with $r=15$, where
the sequence of pivots ${\mathbf{l}}_{r}$ contains these pivots. Beyond the
market factor 1, the sector specific factors 3, 6 and 7 mainly load on firms
in, respectively, the energy, the finance and the health care sector. The
remaining factors are weak factors with very sparse loadings.

##  7 Concluding remarks

We have estimated (from a Bayesian viewpoint) a fairly important and highly
implemented class of sparse factor models when the number of common factors is
unknown. Our framework leads to a natural, efficient and simultaneous coupling
of model estimation and selection on one hand and model identification and
rank estimation (number of factors) on the other hand. More precisely, by
combining point-mass mixture priors with overfitting sparse factor modelling,
in an unordered generalised lower triangular loadings representation
(Frühwirth-Schnatter et al., 2022), we obtain posterior summaries regarding
factor loadings, common factors as well as the factor dimension via
postprocessing draws from our highly efficient and customised MCMC scheme.

The new framework is readily available for some straightforward extensions.
Relatively immediate extensions are (i) idiosyncratic errors following
Student’s $t$-distributions or more general Gaussian mixtures and (ii) dynamic
sparse factor models with stationary common factors; both extensions commonly
found in econometrics applications, see e.g. the recent papers by Piatek and
Papaspiliopoulos (2018) and Kaufmann and Schuhmacher (2019). A further
interesting extension would be to design a prior on the sparsity matrix that a
priori distinguishes between pervasive factors that loads on (nearly) all
measurements, group specific factors that load on selected measurements and
factors that mainly capture (weak) cross-sectional heterogeneity which is not
built into the basic factor model. Such approximate factor models are very
popular in non-Bayesian factor analysis, see e.g. Chamberlain and Rothschild
(1983) and Bai and Ng (2002) and would deserve more attention from the
Bayesian community. However, we leave this interesting idea for future
research.

## References

  * Akaike (1987) Akaike, H. (1987).  Factor analysis and AIC.  Psychometrika 52, 317–332.
  * Anderson (2003) Anderson, T. W. (2003).  An Introduction to Multivariate Statistical Analysis (3 ed.).  Chichester: Wiley.
  * Anderson and Rubin (1956) Anderson, T. W. and H. Rubin (1956).  Statistical inference in factor analysis.  In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume V, pp. 111–150.
  * Aßmann et al. (2016) Aßmann, C., J. Boysen-Hogrefe, and M. Pape (2016).  Bayesian analysis of static and dynamic factor models: An ex-post approach toward the rotation problem.  Journal of Econometrics 192, 190–206.
  * Bai and Ng (2002) Bai, J. and S. Ng (2002).  Determining the number of factors in approximate factor models.  Econometrica 70, 191–221.
  * Bai and Ng (2013) Bai, J. and S. Ng (2013).  Principal components estimation and identification of static factors.  Journal of Econometrics 176, 18–29.
  * Bartholomew (1987) Bartholomew, D. J. (1987).  Latent Variable Models and Factor Analysis.  London: Charles Griffin.
  * Bhattacharya and Dunson (2011) Bhattacharya, A. and D. Dunson (2011).  Sparse Bayesian infinite factor models.  Biometrika 98, 291–306.
  * Boivin and Ng (2006) Boivin, J. and S. Ng (2006).  Are more data always better for factor analysis?  Journal of Econometrics 132, 169–194.
  * Cadonna et al. (2020) Cadonna, A., S. Frühwirth-Schnatter, and P. Knaus (2020).  Triple the gamma – A unifying shrinkage prior for variance and variable selection in sparse state space and TVP models.  Econometrics 8, 20.
  * Carvalho et al. (2008) Carvalho, C. M., J. Chang, J. E. Lucas, J. Nevins, Q. Wang, and M. West (2008).  High-dimensional sparse factor modeling: Applications in gene expression genomics.  Journal of the American Statistical Association 103, 1438–1456.
  * Chamberlain and Rothschild (1983) Chamberlain, G. and M. Rothschild (1983).  Arbitrage, factor structure, and mean-variance analysis on large asset markets.  Econometrica 51, 1281–1304.
  * Chan et al. (2018) Chan, J., R. Leon-Gonzalez, and R. W. Strachan (2018).  Invariant inference and efficient computation in the static factor model.  Journal of the American Statistical Association 113, 819–828.
  * Conti et al. (2014) Conti, G., S. Frühwirth-Schnatter, J. J. Heckman, and R. Piatek (2014).  Bayesian exploratory factor analysis.  Journal of Econometrics 183, 31–57.
  * Durante (2017) Durante, D. (2017).  A note on the multiplicative gamma process.  Statistics and Probability Letters 122, 198–204.
  * Fan et al. (2008) Fan, J., Y. Fan, and J. Lv (2008).  High dimensional covariance matrix estimation using a factor model.  Journal of Econometrics 147, 186–197.
  * Forni et al. (2009) Forni, M., D. Giannone, M. Lippi, and L. Reichlin (2009).  Opening the black box: Structural factor models with large cross sections.  Econometric Theory 25, 1319–1347.
  * Foster and George (1994) Foster, D. P. and E. I. George (1994).  The risk inflation criterion for multiple regression.  The Annals of Statistics 22, 1947–1975.
  * Frühwirth-Schnatter (2022) Frühwirth-Schnatter, S. (2022).  Generalized cumulative shrinkage process priors with applications to sparse bayesian factor analysis.
  * Frühwirth-Schnatter et al. (2022) Frühwirth-Schnatter, S., D. Hosszejni, and H. Lopes (2022).  When it counts - Econometric identification of factor models based on GLT structures.  ArXiv forthcoming soon.
  * Frühwirth-Schnatter and Lopes (2010) Frühwirth-Schnatter, S. and H. Lopes (2010).  Parsimonious Bayesian Factor Analysis when the Number of Factors is Unknown.  Research report, Booth School of Business, University of Chicago.
  * Frühwirth-Schnatter and Lopes (2018) Frühwirth-Schnatter, S. and H. Lopes (2018).  Sparse Bayesian Factor Analysis when the Number of Factors is Unknown.  arXiv 1804.04231.
  * Frühwirth-Schnatter and Tüchler (2008) Frühwirth-Schnatter, S. and R. Tüchler (2008).  Bayesian parsimonious covariance estimation for hierarchical linear mixed models.  Statistics and Computing 18, 1–13.
  * Frühwirth-Schnatter and Wagner (2010) Frühwirth-Schnatter, S. and H. Wagner (2010).  Stochastic model specification search for Gaussian and partially non-Gaussian state space models.  Journal of Econometrics 154, 85–100.
  * Geweke and Singleton (1980) Geweke, J. F. and K. J. Singleton (1980).  Interpreting the likelihood ratio statistic in factor models when sample size is small.  Journal of the American Statistical Association 75, 133–137.
  * Geweke and Zhou (1996) Geweke, J. F. and G. Zhou (1996).  Measuring the pricing error of the arbitrage pricing theory.  Review of Financial Studies 9, 557–587.
  * Ghahramani et al. (2007) Ghahramani, Z., T. L. Griffiths, and P. Sollich (2007).  Bayesian nonparametric latent feature models (with discussion and rejoinder).  In J. M. Bernardo, M. J. Bayarri, J. O. Berger, A. P. Dawid, D. Heckerman, A. F. M. Smith, and M. West (Eds.), Bayesian Statistics 8, pp. ADD–ADD. Oxford: Oxford University Press.
  * Ghosh and Dunson (2009) Ghosh, J. and D. B. Dunson (2009).  Default prior distributions and efficient posterior computation in Bayesian factor analysis.  Journal of Computational and Graphical Statistics 18, 306–320.
  * Griffiths and Ghahramani (2006) Griffiths, T. L. and Z. Ghahramani (2006).  Infinite latent feature models and the Indian buffet process.  In Y. Weiss, B. Schölkopf, and J. Platt (Eds.), Advances in Neural Information Processing Systems, Volume 18, pp. 475–482. Cambridge, MA: MIT Press.
  * Hosszejni and Frühwirth-Schnatter (2022) Hosszejni, D. and S. Frühwirth-Schnatter (2022).  Cover it up! Bipartite graphs uncover identifiability in sparse factor analysis.  arXiv 2211.00671.
  * Jöreskog (1967) Jöreskog, K. G. (1967).  Some contributions to maximum likelihood factor analysis.  Psychometrika 32, 443–482.
  * Jöreskog (1969) Jöreskog, K. G. (1969).  A general approach to confirmatory maximum likelihood factor analysis.  Psychometrika 34, 183–202.
  * Kastner (2019) Kastner, G. (2019).  Sparse Bayesian time-varying covariance estimation in many dimensions.  Journal of Econometrics 210, 98–115.
  * Kastner et al. (2017) Kastner, G., S. Frühwirth-Schnatter, and H. F. Lopes (2017).  Efficient Bayesian inference for multivariate factor stochastic volatility models.  Journal of Computational and Graphical Statistics 26, 905–917.
  * Kaufmann and Schuhmacher (2017) Kaufmann, S. and C. Schuhmacher (2017).  Identifying relevant and irrelevant variables in sparse factor models.  Journal of Applied Econometrics 32, 1123–1144.
  * Kaufmann and Schuhmacher (2019) Kaufmann, S. and C. Schuhmacher (2019).  Bayesian estimation of sparse dynamic factor models with order-independent and ex-post identification.  Journal of Econometrics 210, 116–134.
  * Koopmans and Reiersøl (1950) Koopmans, T. C. and O. Reiersøl (1950).  The identification of structural characteristics.  The Annals of Mathematical Statistics 21, 165–181.
  * Lee and Song (2002) Lee, S.-Y. and X.-Y. Song (2002).  Bayesian selection on the number of factors in a factor analysis model.  Behaviormetrika 29, 23–39.
  * Legramanti et al. (2020) Legramanti, S., D. Durante, and D. B. Dunson (2020).  Bayesian cumulative shrinkage for infinite factorizations.  Biometrika 107, 745–752.
  * Lopes and West (2004) Lopes, H. F. and M. West (2004).  Bayesian model assessment in factor analysis.  Statistica Sinica 14, 41–67.
  * Lucas et al. (2006) Lucas, J., C. Carvalho, Q. Wang, A. Bild, J. R. Nevins, and M. West (2006).  Sparse statistical modelling in gene expression genomics.  In K. Do, P. Müller, and M. Vannucci (Eds.), Bayesian Inference for Gene Expression and Proteomics, pp. 155–176. Cambridge, UK: Cambridge University Press.
  * Martin and McDonald (1975) Martin, J. K. and R. P. McDonald (1975).  Bayesian estimation in unrestricted factor analysis: a treatment for Heywood cases.  Psychometrika 40, 505–517.
  * Neudecker (1990) Neudecker, H. (1990).  On the identification of restricted factor loading matrices: An alternative condition.  Journal of Mathematical Psychology 34, 237–241.
  * O’Hagan (1995) O’Hagan, A. (1995).  Fractional Bayes factors for model comparison.  Journal of the Royal Statistical Society, Ser. B 57, 99–138.
  * Owen and Wang (2016) Owen, A. B. and J. Wang (2016).  Bi-Cross-Validation for Factor Analysis.  Statistical Science 31, 119–139.
  * Paisley and Carin (2009) Paisley, J. W. and L. Carin (2009).  Nonparametric factor analysis with beta process priors.  In Proceedings of the 26th International Conference on Machine Learning (ICML ’09), pp. 777–784.
  * Papastamoulis and Ntzoufras (2022) Papastamoulis, P. and I. Ntzoufras (2022).  On the identifiability of Bayesian factor analytic models.  Statistics and Computing 32.
  * Pati et al. (2014) Pati, D., A. Bhattacharya, N. S. Pillai, and D. B. Dunson (2014).  Posterior contraction in sparse Bayesian factor models for massive covariance matrices.  Annals of Statistics 42, 1102–1130.
  * Piatek and Papaspiliopoulos (2018) Piatek, R. and O. Papaspiliopoulos (2018).  A bayesian nonparametric approach to factor analysis.  Submitted.
  * Poworoznek et al. (2021) Poworoznek, E., F. Ferrari, and D. Dunson (2021).  Efficiently resolving rotational ambiguity in Bayesian matrix sampling with matching.  ArXiv 2107.13783.
  * Reiersøl (1950) Reiersøl, O. (1950).  On the identifiability of parameters in Thurstone’s multiple factor analysis.  Psychometrika 15, 121–149.
  * Ročková and George (2017) Ročková, V. and E. I. George (2017).  Fast Bayesian factor analysis via automatic rotation to sparsity.  Journal of the American Statistical Association 111, 1608–1622.
  * Rue and Held (2005) Rue, H. and L. Held (2005).  Gaussian Markov Random Fields: Theory and Applications, Volume 104 of Monographs on Statistics and Applied Probability.  London: Chapman & Hall/CRC.
  * Stock and Watson (2002) Stock, J. H. and M. W. Watson (2002).  Macroeconomic forecasting using diffusion indexes.  Journal of Business & Economic Statistics 20, 147–162.
  * Teh et al. (2007) Teh, Y. W., D. Görür, and Z. Ghahramani (2007).  Stick-breaking construction for the Indian buffet process.  In M. Meila and X. Shen (Eds.), Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, Volume 2 of Proceedings of Machine Learning Research, San Juan, Puerto Rico, pp. 556–563. PMLR.
  * Thurstone (1947) Thurstone, L. L. (1947).  Multiple factor analysis.  Chicago: University of Chicago.
  * Tumura and Sato (1980) Tumura, Y. and M. Sato (1980).  On the identification in factor analysis.  TRU Mathematics 16, 121–131.
  * van Dyk and Meng (2001) van Dyk, D. and X.-L. Meng (2001).  The art of data augmentation.  Journal of Computational and Graphical Statistics 10, 1–50.
  * West (2003) West, M. (2003).  Bayesian factor regression models in the “large p, small n” paradigm.  In J. M. Bernardo, M. J. Bayarri, J. O. Berger, A. P. Dawid, D. Heckerman, A. F. M. Smith, and M. West (Eds.), Bayesian Statistics 7, pp. 733–742. Oxford: Oxford University Press.
  * Williams (2020) Williams, B. (2020).  Identification of the linear factor model.  Econometric Reviews 39, 92–109.
  * Yu and Meng (2011) Yu, Y. and X.-L. Meng (2011).  To center or not to center: That is not the question - An ancillarity-suffiency interweaving strategy (ASIS) for boosting MCMC efficiency.  Journal of Computational and Graphical Statistics 20, 531–615.
  * Zhao et al. (2016) Zhao, S., C. Gao, S. Mukherjee, and B. E. Engelhardt (2016).  Bayesian group factor analysis with structured sparsity.  Journal of Machine Learning Research 17, 1–47.

Supplementary material for:
“Sparse Bayesian Factor Analysis when the Number of Factors is Unknown”

##  Appendix A Details on Step (P)

###  A.1 Posterior distributions in a confirmatory factor model

Step (P) of Algorithm 1 updates the parameters in the confirmatory sparse
factor model factor model

$\displaystyle{\mathbf{f}}^{r}_{t}\sim\mathcal{N}_{r}\left({\mathbf{0}},{{\mathbf{I}}}_{r}\right),\quad{\mathbf{y}}_{t}=\boldsymbol{\beta}_{r}{\mathbf{f}}^{r}_{t}+\boldsymbol{\epsilon}_{t},\quad\boldsymbol{\epsilon}_{t}\sim\mathcal{N}_{m}\left({\mathbf{0}},{\mathbf{\Sigma}}_{r}\right),\quad{\mathbf{\Sigma}}_{r}=\mbox{\rm
Diag}\,(\sigma^{2}_{1},\ldots,\sigma^{2}_{m}),$ (A.1)

where the indicator matrix $\boldsymbol{\delta}_{r}$ imposes a certain zero
structure on the loading matrix $\boldsymbol{\beta}_{r}$. The joint posterior
distribution
$p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i}|{\mathbf{y}},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r})$
of the nonzero factor loadings
$\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}$ and the idiosyncratic
variance $\sigma^{2}_{i}$ is derived for each row $i$ ($i=1,\ldots,m$)
conditional on the factors
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$ and the
indicator matrix $\boldsymbol{\delta}_{r}$ from the following regression
model:

$\displaystyle\tilde{{\mathbf{y}}}_{i}={\mathbf{X}}_{i}^{\boldsymbol{\delta}}\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}+\tilde{\boldsymbol{\epsilon}}_{i},$
(A.2)

where $\tilde{{\mathbf{y}}}_{i}=(y_{i1}\cdots y_{iT})^{{}^{\prime}}$ and
$\tilde{\boldsymbol{\epsilon}}_{i}=(\epsilon_{i1}\cdots\epsilon_{iT})^{{}^{\prime}}\sim\mathcal{N}_{T}\left(0,\sigma^{2}_{i}{\mathbf{I}}\right)$.
${\mathbf{X}}_{i}^{\boldsymbol{\delta}}$ is a regressor matrix for
$\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}$ constructed from the
$T\times r$ dimensional latent factor matrix
${\mathbf{F}}=({\mathbf{f}}^{r}_{1}\cdots{\mathbf{f}}^{r}_{T})^{{}^{\prime}}$
in the following way. If no element in row $i$ of $\boldsymbol{\beta}_{r}$ is
restricted to 0, then ${\mathbf{X}}_{i}^{\boldsymbol{\delta}}={\mathbf{F}}$.
If some elements are restricted to 0, then
${\mathbf{X}}_{i}^{\boldsymbol{\delta}}$ is obtained from ${\mathbf{F}}$ by
deleting all columns $j$ where $\delta_{ij}=0$, i.e.
${\mathbf{X}}_{i}^{\boldsymbol{\delta}}={\mathbf{F}}\boldsymbol{\Pi}_{i}^{\boldsymbol{\delta}}$,
where $\boldsymbol{\Pi}_{i}^{\boldsymbol{\delta}}$ is a
$r\times\sum_{j=1}^{r}\delta_{ij}$ selection matrix, selecting those columns
$j$ of ${\mathbf{F}}$ where $\delta_{ij}\neq 0$. The likelihood derived from
(A.2) is combined with the inverted Gamma prior (3.18) on $\sigma^{2}_{i}$
and, respectively, the hierarchical Gaussian prior (3.13) or the fractional
prior (3.17) for
$\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i}$. In a
sparse factor model, the dimension of this posterior depends on the number of
nonzero elements in the $i$th row of $\boldsymbol{\beta}_{r}$, i.e.
$q_{i}=\sum_{j=1}^{r}\delta_{ij}$. There are basically three types of rows,
when it comes to updating the parameters: zero rows, dedicated rows and rows
with multiple loadings.

For zero rows (i.e. $q_{i}=0$), (A.2) reduces to a “null” model without
regressors ${\mathbf{X}}_{i}^{\boldsymbol{\delta}}$, that is
$\tilde{{\mathbf{y}}}_{i}=\tilde{\boldsymbol{\epsilon}}_{i}$. Hence, the
posterior of $\sigma^{2}_{i}$ is simply given by

$\displaystyle\sigma^{2}_{i}|\tilde{{\mathbf{y}}}_{i},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r}\sim\mathcal{G}^{-1}\left(c_{T}^{n},C_{iT}^{n}\right),\qquad
c_{T}^{n}=c_{0}+\frac{T}{2},\qquad
C_{iT}^{n}=C_{i0}+\frac{1}{2}\sum_{t=1}^{T}y_{it}^{2}.$ (A.3)

For all nonzero rows (i.e. $q_{i}>0)$, the posterior
$(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})$ for a
specific row $i$ is given by:

$\displaystyle\sigma^{2}_{i}|\tilde{{\mathbf{y}}}_{i},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r}\sim\mathcal{G}^{-1}\left(c_{T},C_{iT}^{\boldsymbol{\delta}}\right),\qquad\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},\tilde{{\mathbf{y}}}_{i},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r}\sim\mathcal{N}_{q_{i}}\left({\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}},{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}\sigma^{2}_{i}\right).$
(A.4)

For the hierarchical Gaussian prior (3.13), the moments are given by:

$\displaystyle({\mathbf{B}}_{iT}^{\boldsymbol{\delta}})^{-1}=({\mathbf{B}}_{i0}^{\boldsymbol{\delta}})^{-1}+({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{X}}_{i}^{\boldsymbol{\delta}},\qquad{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}=({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i},$
(A.5) $\displaystyle c_{T}=c_{0}+\frac{T}{2},\qquad
C_{iT}^{\boldsymbol{\delta}}=C_{i0}+\frac{1}{2}{\mbox{\rm
SSR}}_{i},\quad{\mbox{\rm SSR}}_{i}=\
\tilde{{\mathbf{y}}}_{i}^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i}-({\mathbf{m}}_{iT}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}.$

For the fractional prior (3.17), the moments are given by:

$\displaystyle({\mathbf{B}}_{iT}^{\boldsymbol{\delta}})^{-1}=({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{X}}_{i}^{\boldsymbol{\delta}},\quad{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}=({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i},$
(A.6) $\displaystyle c_{T}=c_{0}+\frac{(1-b)T}{2},\quad
C_{iT}^{\boldsymbol{\delta}}=C_{i0}+\frac{(1-b)}{2}{\mbox{\rm SSR}}_{i},$

where ${\mbox{\rm SSR}}_{i}$ is the same as in (A.5) and, for the fractional
prior, identical to the residual sum of squares errors.888If the residual
$\boldsymbol{\epsilon}_{i}=\tilde{{\mathbf{y}}}_{i}-{\mathbf{X}}_{i}^{\boldsymbol{\delta}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}$
is defined in the usual way, then:
$\displaystyle\boldsymbol{\epsilon}^{{}^{\prime}}_{i}\boldsymbol{\epsilon}_{i}$
$\displaystyle=$
$\displaystyle\tilde{{\mathbf{y}}}_{i}^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i}-({\mathbf{m}}_{iT}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i}-\tilde{{\mathbf{y}}}_{i}^{{}^{\prime}}{\mathbf{X}}_{i}^{\boldsymbol{\delta}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}+({\mathbf{m}}_{iT}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}({\mathbf{X}}_{i}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{X}}_{i}^{\boldsymbol{\delta}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}$
$\displaystyle=$
$\displaystyle\tilde{{\mathbf{y}}}_{i}^{{}^{\prime}}\tilde{{\mathbf{y}}}_{i}-({\mathbf{m}}_{iT}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}-({\mathbf{m}}_{iT}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}+({\mathbf{m}}_{iT}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}{\mathbf{m}}_{iT}^{\boldsymbol{\delta}}={\mbox{\rm
SSR}}_{i}.$ For dedicated rows (i.e. $q_{i}=1$) only a single nonzero factor
loading $\beta_{i,{j_{i}}}$ is present in a particular column ${j_{i}}$ and
the posterior given in (A.4) simplifies considerably:

$\displaystyle\sigma^{2}_{i}|\tilde{{\mathbf{y}}}_{i},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r}\sim\mathcal{G}^{-1}\left(c_{T},C_{iT}\right),\qquad\beta_{i,{j_{i}}}|\sigma^{2}_{i},\tilde{{\mathbf{y}}}_{i},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r}\sim\mathcal{N}\left(B_{iT}m_{iT},B_{iT}\sigma^{2}_{i}\right).$
(A.7)

For a hierarchical Gaussian prior, the posterior moments are given by:

$\displaystyle
B_{iT}=1/(B_{i0,{j_{i}}{j_{i}}}^{-1}+\sum_{t=1}^{T}f_{{j_{i}},t}^{2}),\quad
m_{iT}=\sum_{t=1}^{T}f_{{j_{i}},t}y_{it},$ (A.8) $\displaystyle
c_{T}=c_{0}+\frac{T}{2},\qquad
C_{iT}=C_{i0}+\frac{1}{2}\left(\sum_{t=1}^{T}y_{it}^{2}-m_{iT}^{2}B_{iT}\right),$

where $B_{i0,{j_{i}}{j_{i}}}$ is the prior variance of the dedicated factor
loading in row $i$, whereas for a fractional prior:

$\displaystyle B_{iT}=1/\left(\sum_{t=1}^{T}f_{{j_{i}},t}^{2}\right),\quad
m_{iT}=\sum_{t=1}^{T}f_{{j_{i}},t}y_{it},$ (A.9) $\displaystyle
c_{T}=c_{0}+\frac{(1-b)T}{2},\quad
C_{iT}=C_{i0}+\frac{(1-b)}{2}\left(\sum_{t=1}^{T}y_{it}^{2}-m_{iT}^{2}B_{iT}\right).$

###  A.2 Block sampling of idiosyncratic variances and factor loadings

Step (P) in Algorithm 1 could be implemented as in Lopes and West (2004), by
sampling $\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}$ and
$\sigma^{2}_{i}$ from the posterior distribution
$p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i}|{\mathbf{y}},{\mathbf{f}}_{r},\boldsymbol{\delta}_{r})$
derived in Section A.1 row by row. However, an important improvement is
feasible through block sampling of all idiosyncratic variances and all nonzero
factor loadings, summarized in Algorithm 3.999This algorithm has been
implemented for the first time in the unpublished research report by
Frühwirth-Schnatter and Lopes (2010). The use of the Cholesky decomposition of
the information matrix (instead of the covariance matrix) to sample from a
high-dimensional density is fashioned after Rue and Held (2005, Theorem 2.5
and Algorithm 2.5) who consider Gaussian random fields.

######  Algorithm 3.

Sampling parameters for a sparse Bayesian factor model

  * (P-a)

For all zero rows, sample $\sigma^{2}_{i}$ from (A.3), which can be trivially
vectorized.

  * (P-b)

If the remaining rows are all dedicated with a single nonzero loading in
column $j_{i}$ (which can be different for different rows), then sampling from
(A.7) is easily vectorized, since all posterior moments are univariate.

  * (P-c)

Even if some of the nonzero rows are not dedicated, joint sampling of all
idiosyncratic variances and all factor loadings is feasible for all nonzero
rows. Let $i_{1},\ldots,i_{n}$ be the indices of all $n=m-m_{0}$ nonzero rows
of $\boldsymbol{\beta}_{r}$, i.e. $q_{i_{l}}>0$ for $l=1,\ldots,n$. Let
$d=\sum_{i}q_{i}$ be the total number of nonzero elements in
$\boldsymbol{\beta}_{r}$ and let
$\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}=(\boldsymbol{\beta}_{i_{1}\cdot}^{\boldsymbol{\delta}}\cdots\boldsymbol{\beta}_{i_{n}\cdot}^{\boldsymbol{\delta}})^{\prime}$
be the $d$-dimensional vector obtained by stacking row by row all nonzero
elements in each row of $\boldsymbol{\beta}_{r}$. To sample the idiosyncratic
variances $\sigma^{2}_{i_{1}},\ldots,\sigma^{2}_{i_{n}}$ and the nonzero
factor loadings $\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}$ in
$\boldsymbol{\beta}_{r}$ jointly, proceed in the following way:

    1. (P-c1)

Construct the information matrix ${\mathbf{P}}$ and the covector
${\mathbf{m}}$ of the joint posterior

$\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}|\sigma^{2}_{i_{1}},\ldots,\sigma^{2}_{i_{n}},{\mathbf{f}}_{r},{\mathbf{y}}\sim\mathcal{N}_{d}\left({\mathbf{P}}^{-1}{\mathbf{m}},{\mathbf{P}}^{-1}\mathbf{D}\right).$

The matrix $\mathbf{D}=\mbox{\rm
Diag}\,(\sigma^{2}_{i_{1}}{\mathbf{1}}_{1\times
q_{i_{1}}}\cdots\sigma^{2}_{i_{n}}{\mathbf{1}}_{1\times q_{i_{n}}})$, with
${\mathbf{1}}_{1\times l}$ being a $1\times l$ row vector of ones, is a
$d\times d$ diagonal matrix containing the idiosyncratic variances, while the
$d\times d$ matrix ${\mathbf{P}}$ and the $d\times 1$ vector ${\mathbf{m}}$
are given by:

$\displaystyle{\mathbf{P}}=\left(\begin{array}[]{cccc}({\mathbf{B}}_{i_{1},T}^{\boldsymbol{\delta}})^{-1}&{\mathbf{O}}&\cdots&{\mathbf{O}}\\\
{\mathbf{O}}&({\mathbf{B}}_{i_{2},T}^{\boldsymbol{\delta}})^{-1}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&{\mathbf{O}}\\\
{\mathbf{O}}&\cdots&{\mathbf{O}}&({\mathbf{B}}_{i_{n},T}^{\boldsymbol{\delta}})^{-1}\\\
\end{array}\right),\qquad{\mathbf{m}}=\left(\begin{array}[]{c}{\mathbf{m}}_{i_{1},T}^{\boldsymbol{\delta}}\\\
\vdots\\\ {\mathbf{m}}_{i_{n},T}^{\boldsymbol{\delta}}\\\ \end{array}\right),$

where $({\mathbf{B}}_{i_{l},T}^{\boldsymbol{\delta}})^{-1}$ and
${\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}}$ are the information matrix and
the covector appearing in the posterior (A.4) of the nonzero elements in row
$i_{l}$. ${\mathbf{P}}$ is a sparse band matrix with maximal band width equal
to $\max q_{i_{l}}$.

    2. (P-c2)

Compute the Cholesky decomposition
${\mathbf{P}}={\mathbf{L}}{\mathbf{L}}^{{}^{\prime}}$, where ${\mathbf{L}}$ is
lower triangular, using a special algorithm developed for band matrices. Next,
solve ${\mathbf{L}}{\mathbf{x}}={\mathbf{m}}$ for ${\mathbf{x}}$ using an
algorithm specially designed for triangular matrices. Evidently,
${\mathbf{x}}$ is a $d\times 1$ vector.

    3. (P-c3)

Sample $\sigma^{2}_{i_{1}},\ldots,\sigma^{2}_{i_{n}}$ jointly from (A.4). The
squared sum ${\mathbf{x}}^{{}^{\prime}}{\mathbf{x}}$ can be used to vectorize
the computation of $C_{i_{l},T}^{\boldsymbol{\delta}}$ for each
$l=1,\ldots,n$, since

$\displaystyle{\mathbf{x}}_{i_{l}}^{{}^{\prime}}{\mathbf{x}}_{i_{l}}=({\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{i_{l},T}^{\boldsymbol{\delta}}{\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}},$
(A.11)

where ${\mathbf{x}}_{i_{l}}$ is the $q_{i_{l}}$-dimensional sub vector of
${\mathbf{x}}$ corresponding to
$\boldsymbol{\beta}_{i_{l},\cdot}^{\boldsymbol{\delta}}$.

    4. (P-c4)

Finally, define the diagonal matrix $\mathbf{D}$ from
$\sigma^{2}_{i_{1}},\ldots,\sigma^{2}_{i_{n}}$ as described above and draw
${\mathbf{z}}\sim\mathcal{N}_{d}\left({\mathbf{0}},\mathbf{D}\right)$. Solving
the system

$\displaystyle{\mathbf{L}}^{{}^{\prime}}\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}={\mathbf{x}}+{\mathbf{z}}$
(A.12)

for $\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}$ leads to a draw from the
joint posterior
$\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}|\sigma^{2}_{i_{1}},\ldots,\sigma^{2}_{i_{n}},{\mathbf{y}},{\mathbf{f}}_{r}$.

To derive (A.11), let ${\mathbf{L}}_{i_{l}}$ be the $q_{i_{l}}\times
q_{i_{l}}$ submatrix of ${\mathbf{L}}$ corresponding to
$\boldsymbol{\beta}_{i_{l},\cdot}^{\boldsymbol{\delta}}$. Evidently,
${\mathbf{L}}_{i_{l}}$ is equal to the Cholesky decomposition of the
individual information matrix
$({\mathbf{B}}_{i_{l},T}^{\boldsymbol{\delta}})^{-1}$. Furthermore, the
$q_{i_{l}}$-dimensional sub vector ${\mathbf{x}}_{i_{l}}$ corresponding to
$\boldsymbol{\beta}_{i_{l},\cdot}^{\boldsymbol{\delta}}$ satisfies
${\mathbf{L}}_{i_{l}}{\mathbf{x}}_{i_{l}}={\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}}$.
Therefore:

$\displaystyle{\mathbf{x}}_{i_{l}}^{{}^{\prime}}{\mathbf{x}}_{i_{l}}=({\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}})^{{}^{\prime}}({\mathbf{L}}_{i_{l}}^{{}^{\prime}})^{-1}{\mathbf{L}}_{i_{l}}^{-1}{\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}}=({\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}})^{{}^{\prime}}({\mathbf{L}}_{i_{l}}{\mathbf{L}}_{i_{l}}^{{}^{\prime}})^{-1}{\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}}=({\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}})^{{}^{\prime}}{\mathbf{B}}_{i_{l},T}^{\boldsymbol{\delta}}{\mathbf{m}}_{i_{l},T}^{\boldsymbol{\delta}}.$

It is easy to prove that the solution
$\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}$ of (A.12) is a draw from the
posterior
$p(\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}|\sigma^{2}_{i_{1}},\ldots,\sigma^{2}_{i_{n}},{\mathbf{y}},{\mathbf{f}})$.
Note that
${\mathbf{L}}{\mathbf{L}}^{{}^{\prime}}\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}={\mathbf{L}}{\mathbf{x}}+{\mathbf{L}}{\mathbf{z}}={\mathbf{m}}+{\mathbf{L}}{\mathbf{z}}$.
Therefore

$\displaystyle\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}=({\mathbf{L}}{\mathbf{L}}^{{}^{\prime}})^{-1}{\mathbf{m}}+({\mathbf{L}}{\mathbf{L}}^{{}^{\prime}})^{-1}{\mathbf{L}}{\mathbf{z}}={\mathbf{P}}^{-1}{\mathbf{m}}+({\mathbf{L}}^{{}^{\prime}})^{-1}{\mathbf{z}}.$

Evidently, $\mbox{\rm
E}(\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r})={\mathbf{P}}^{-1}{\mathbf{m}}$.
Since for each $l=1,\ldots,n$,
${\mathbf{L}}_{i_{l}}\sigma^{2}_{i_{l}}=\sigma^{2}_{i_{l}}{\mathbf{L}}_{i_{l}}$,
it holds that ${\mathbf{L}}\mathbf{D}=\mathbf{D}{\mathbf{L}}$ and therefore
$\mathbf{D}{\mathbf{L}}^{-1}={\mathbf{L}}^{-1}\mathbf{D}$. Since $\mbox{\rm
V}(\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r})=({\mathbf{L}}^{{}^{\prime}})^{-1}\mathbf{D}{\mathbf{L}}^{-1}=({\mathbf{L}}^{{}^{\prime}})^{-1}{\mathbf{L}}^{-1}\mathbf{D}={\mathbf{P}}^{-1}\mathbf{D}$,
it follows that
$\boldsymbol{\beta}^{\boldsymbol{\delta}}_{r}\sim\mathcal{N}_{d}\left({\mathbf{P}}^{-1}{\mathbf{m}},{\mathbf{P}}^{-1}\mathbf{D}\right)$.

##  Appendix B Details on Step (D)

###  B.1 Marginal likelihoods when the factors are known

Although we work throughout this paper with a factor model where the factors
${\mathbf{f}}^{r}_{t}$ are latent, several steps of Algorithm 1 perform model
selection with respect to $\boldsymbol{\delta}_{r}$ conditional on the most
recent draw of the factors
${\mathbf{f}}_{r}=({\mathbf{f}}^{r}_{1},\ldots,{\mathbf{f}}^{r}_{T})$ in the
confirmatory factor model (A.1). Hence, to sample new indicators
$\boldsymbol{\delta}_{i\cdot}$ in row $i$, the marginal likelihood
$p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\delta}_{r})$ of
regression model (A.2) is needed.

If $\boldsymbol{\delta}_{i\cdot}$ is a zero row (i.e $q_{i}=0$), then the
marginal likelihood simplifies to

$\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\delta}_{r})=p(\tilde{{\mathbf{y}}}_{i})=\frac{\Gamma(c_{T}^{n})(C_{i0})^{c_{0}}}{(2\pi)^{T/2}\Gamma(c_{0})(C_{iT}^{n})^{c_{T}^{n}}},$
(B.1)

where $c_{T}^{n}$ and $C_{iT}^{n}$ are the posterior moments of
$\sigma^{2}_{i}$ under the “null” model given by (A.3).

If at least one element of $\boldsymbol{\delta}_{i\cdot}$ is different from
zero, then the marginal likelihood computation differs between the
hierarchical Gaussian prior (3.13) and the fractional prior (3.17).

#### Marginal likelihoods for a hierarchical Gaussian prior.

For a hierarchical Gaussian prior, a well-known exercise in Bayesian
regression analysis yields:

$\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|\boldsymbol{\delta}_{r},{\mathbf{f}}_{r})=\frac{1}{(2\pi)^{T/2}}\frac{|{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}|^{1/2}}{|{\mathbf{B}}_{i0}^{\boldsymbol{\delta}}|^{1/2}}\frac{\Gamma(c_{T})(C_{i0})^{c_{0}}}{\Gamma(c_{0})(C_{iT}^{\boldsymbol{\delta}})^{c_{T}}},$
(B.2)

where ${\mathbf{B}}_{iT}^{\boldsymbol{\delta}}$, $c_{T}$ and
$C_{iT}^{\boldsymbol{\delta}}$ are the posterior moments of
$p(\boldsymbol{\beta}^{\boldsymbol{\delta}}_{i\cdot},\sigma^{2}_{i}|\boldsymbol{\delta}_{r},\tilde{{\mathbf{y}}}_{i},{\mathbf{f}}_{r})$
given by (A.5).

#### Marginal likelihoods for a fractional prior.

For a fractional prior, the derivation of the marginal likelihood if at least
one element of $\boldsymbol{\delta}_{i\cdot}$ is different from zero, is less
standard and can be obtained in a similar way as in Frühwirth-Schnatter and
Wagner (2010). A fraction $b$ of the full conditional likelihood of regression
model (A.2) is used to define the fractional prior
$p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r})$
in (3.17):

$\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})=p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{1-b}p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{b}\propto
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{1-b}p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r}).$

The remaining part of the likelihood, that is
$p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{1-b}$,
is used for model selection and is combined with the prior $p(\sigma^{2}_{i})$
defined in (3.18) and the normalized fractional prior
$p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r})$,
given by:

$\displaystyle
p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r})=\frac{p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{b}}{c_{i}(\sigma^{2}_{i},{\mathbf{f}}_{r},b)}.$

The normalizing constant $c_{i}(\sigma^{2}_{i},{\mathbf{f}}_{r},b)$ is given
by:

$\displaystyle c_{i}(\sigma^{2}_{i},{\mathbf{f}}_{r},b)=\int
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{b}\,d\,\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}=(2\pi\sigma^{2}_{i})^{\frac{q_{i}-Tb}{2}}b^{-\frac{q_{i}}{2}}|{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}|^{1/2}\exp\left(-\frac{b}{2\sigma^{2}_{i}}{\mbox{\rm
SSR}}_{i}\right),$ (B.3)

where ${\mathbf{B}}_{iT}^{\boldsymbol{\delta}}$ and ${\mbox{\rm SSR}}_{i}$ are
the posterior moments of
$p(\boldsymbol{\beta}^{\boldsymbol{\delta}}_{i\cdot},\sigma^{2}_{i}|{\mathbf{f}}_{r},\boldsymbol{\delta}_{i\cdot},\tilde{{\mathbf{y}}}_{i})$
given by (A.6). Integrating the fractional posterior

$\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{1-b}p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|\sigma^{2}_{i},b,{\mathbf{f}}_{r})p(\sigma^{2}_{i})$

over $\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}$, yields the
fractional likelihood
$p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\sigma^{2}_{i},b)$:

$\displaystyle\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\sigma^{2}_{i},b)$
$\displaystyle=$ $\displaystyle\int
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})^{1-b}p(\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}|{\mathbf{f}}_{r},\sigma^{2}_{i},b)\,d\,\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}=\frac{1}{c_{i}(\sigma^{2}_{i},{\mathbf{f}}_{r},b)}\int
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}},\sigma^{2}_{i})\,d\,\boldsymbol{\beta}_{i\cdot}^{\boldsymbol{\delta}}$
$\displaystyle=$
$\displaystyle\displaystyle\left(\frac{1}{2\pi\sigma^{2}_{i}}\right)^{\frac{(T-q_{i})-(Tb-
q_{i})}{2}}b^{\frac{q_{i}}{2}}\frac{|{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}|^{1/2}}{|{\mathbf{B}}_{iT}^{\boldsymbol{\delta}}|^{1/2}}\exp\left(-\frac{(1-b)}{2\sigma^{2}_{i}}{\mbox{\rm
SSR}}_{i}\right)$ $\displaystyle=$
$\displaystyle\displaystyle\left(\frac{1}{2\pi\sigma^{2}_{i}}\right)^{\frac{T(1-b)}{2}}b^{\frac{q_{i}}{2}}\exp\left(-\frac{(1-b)}{2\sigma^{2}_{i}}{\mbox{\rm
SSR}}_{i}\right).$

When we combine
$p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\sigma^{2}_{i},b)$ with the prior
$p(\sigma^{2}_{i})$, then we obtain:

$\displaystyle
p(\tilde{{\mathbf{y}}}_{i}|{\mathbf{f}}_{r},\sigma^{2}_{i},b)p(\sigma^{2}_{i})=\frac{C_{i0}^{c_{0}}}{\Gamma(c_{0})}\left(\frac{1}{2\pi}\right)^{\frac{T(1-b)}{2}}b^{\frac{q_{i}}{2}}\left(\frac{1}{\sigma^{2}_{i}}\right)^{\frac{c_{0}+T(1-b)}{2}}\exp\left(-\frac{C_{i0}+{\mbox{\rm
SSR}}_{i}(1-b)/2}{\sigma^{2}_{i}}\right),$

which is the kernel of the inverted Gamma distribution in (A.4). Integrating