Quasi Maximum Likelihood Estimation of High-Dimensional Factor Models: A Critical Review

Quasi Maximum Likelihood Estimation of High-Dimensional Factor Models: A Critical Review

\thispagestyle{empty}

\footnotetext{Department of Economics - Universit\`a di Bologna. Email: [email removed] \\ Sections (ref) and (ref) heavily draw on BLqml and MBQMLPCA to which the reader is referred for the original results under a more general setting, and for their proof.\\ \\ I thank Matteo Luciani and Rolf Sundberg for helpful comments.\\ \\ A shorter version of this paper is published as: M. Barigozzi, “Quasi Maximum Likelihood Estimation of High-Dimensional Factor Models”, in Oxford Encyclopedia of Economics and Finance, Oxford University Press, 2024. }

Introduction

Factor analysis is one of the earliest proposed multivariate statistical techniques: it dates back to the studies of spearman04 in experimental psychology. The idea is to decompose a vector of $N$ observed random variables into two components: (i) a common component driven by $r<N$ latent factors, and (ii) a component which is idiosyncratic, in the sense that it can be due either to measurement errors or to specific features of a single, or a sub-group of variables. As such we can retrospectively consider factor analysis as a pioneering technique in the filed of unsupervised statistical learning.

Over the years, two main strands of applications have emerged: first, in psychometrics in a low-dimensional setting thomson36,bartlett37,bartlett38,lawley1940,thomson51,joreskog69,lawleymaxwell71,moustaki2000generalized,joreskog2001factor,bartholomew2011latent, second, in econometrics, both in a low-dimensional SS77,GS81,watsonengle83,stockwatson89,stockwatson91, and in a high-dimensional setting chamberlainrothschild83,FHLR00,stockwatson02JASA,bai03. Applications in econometrics include: the analysis of financial markets connor2006common,ait2017using,kim2019factor,barigozzi2020generalized, the measurement and prediction of macroeconomic aggregates stockwatson02JBES,de2008forecasting,Nowcasting,OGAP, the study of the dynamic effects of unexpected shocks to the economy BBE05,FGLR09,fornigambettiJME,BCL,FGLS14,smokinggun,marcellino16,BLL2, and the analysis of demand systems stone1945analysis,barigozzi2016identifying. An historical overview is in BH23.

The former strand of literature deals with random samples of data, i.e., independent observations for each variable, and estimation is almost exclusively by means of Maximum Likelihood AR56,AFP87,AA88,tippingbishop99,baili12. The latter strand of literature deals with time series data, i.e., dependent observations for each variable, and estimation is either by means of Principal Component analysis FHLR00,stockwatson02JASA,bai03,FLM13 or by means of Quasi Maximum Likelihood quahsargent93,DGRqml,baili16. Principal Components and Quasi Maximum Likelihood estimation are the most popular estimation techniques, dating back at least to hotelling1933analysis and lawley1940, respectively. A third, less popular, approach is the centroid method by thomson51.

Our focus is on high-dimensional factor models for time series. Specifically, we allow for the dimension of the vector, $N$, to grow to infinity. As a consequence, it is realistic to consider an {\it approximate} factor model where we allow for the idiosyncratic components to be cross-sectionally correlated, rather than an {\it exact} factor model where the idiosyncratic components are not correlated. Furthermore, since we deal with time series, both the factors and the idiosyncratic components are allowed to be autocorrelated, and, if we specify also a model describing the time evolution of the factors, e.g., a VAR, we speak of a {\it dynamic} factor model, while, if we do not, we speak of a {\it static} factor model. The approximate dynamic factor model considered in this paper is a restricted version of the more general model proposed by scherrer1998structure and FHLR00, which is, in fact, a representation of an infinite dimensional panel of time series as derived by fornilippi01 and hallinlippi13.

The idea of an approximate factor structure dates back to chamberlain83 and chamberlainrothschild83 who showed that, if we let the number of time series, $N$, to grow, then we are able to identify the common components even when allowing for idiosyncratic correlations. Intuitively, since estimation is based on cross-sectional aggregation, the more series we use, i.e., the more we increase $N$, the more the signal, factor-driven, component emerges from the noisy idiosyncratic component because of an assumed cross-sectional weak Law of Large Numbers. Equivalently, this is seen by noticing that, under a factor model, the gap between the largest $r$ eigenvalues of the covariance matrix of the data and the remaining $N-r$ widens as $N$ grows. As a consequence of this phenomenon, an approximate factor model does not suffer of the {\it curse of dimensionality} problem, which typically affects high-dimensional parametric models, but it enjoys instead a {\it blessing of dimensionality} property.

The idea of a dynamic factor model goes back to geweke1977dynamic and SS77, who account for the possibility that the $N$ variables co-move not only contemporaneously but also with some leads or lags, in other words, they can load the factors also with lags. In particular, the setting considered here is the same as the one formulated by FGLR09, which consists in writing the model as a linear state-space system, having as latent states the factors.

In this paper, we review likelihood-based methods for estimating high-dimensional approximate static and dynamic factor models, as considered by baili16 and DGRqml, respectively. On the one hand, baili16 treat the factors as non-stochastic and follow the classical approach based on first estimating the loadings and the idiosyncratic variances by Quasi Maximum Likelihood and then estimating the factors by Weighted Least Squares. On the other hand, DGRqml explicitly model the factors as an autocorrelated stochastic process and propose to use the Expectation Maximization (EM) algorithm together with the Kalman smoother to obtain joint estimates of the loadings and the idiosyncratic variances, which approximate the Quasi Maximum Likelihood estimators, and of the factors. Both approaches are based on the maximization of a mis-specified log-likelihood where we do not take into account any of the idiosyncratic cross- and serial correlations, nevertheless, we show that as $N\to\infty$, the mis-specification error vanishes: yet another instance of the blessing of dimensionality. A summary of the rates of convergence of the loadings and factors estimators in approximate and exact factor models is in Table (ref).

The approximate factor model

Let us consider the $T\times N$ double-array $\bm X$ having elements $\{x_{it}\in ,\, i=1,\ldots, N,\, t=1,\ldots,T\}$ of $N$ observed times series over $T$ time periods. We say that $\bm X$ follows an $r$-factor model if

where $\bm\lambda_i:=(\lambda_{i1}\cdots\lambda_{ir})^\prime$ and $\mathbf F_t:=(F_{1t}\cdots F_{rt})^\prime$ are the $r$-dimensional latent vectors of loadings for series $i$ and factors, respectively. We call $\xi_{it}$ the idiosyncratic component and $\chi_{it}:=\bm \lambda_i^\prime \mathbf F_t$ the common component.

It is convenient to write model (ref), which is written is scalar notation, also in other equivalent ways which we list below.

Let $\mathbf x_t:=(x_{1t}\cdots x_{Nt})^\prime$ be the $N$-dimensional vector of observables at time $t$. Then, in vector notation (ref) reads \[ \mathbf x_t =\bm\alpha+ \bm\Lambda\mathbf F_t +\bm \xi_t,\quad t=1\ldots, T, \] where $\bm\xi_{t}:=(\xi_{1t}\cdots\xi_{Nt})^\prime$ is the corresponding $N$-dimensional vector of idiosyncratic components, $\bm\Lambda:=(\bm\lambda_1\cdots\bm\lambda_N)^\prime=(\bm\ell_1\cdots \bm\ell_r)$ is the $n\times r$ matrix of factor loadings with rows $\bm\lambda_i^\prime$, $i=1,\ldots, N$ and columns $\bm\ell_j$, $j=1,\ldots ,r$, and $\bm\alpha:=(\alpha_1\cdots \alpha_N)^\prime$. We call $\bm\chi_t:=\bm\Lambda\mathbf F_t$ the vector of common components.

Equivalently, let $\bm x_i:=(x_{i1}\cdots x_{iT})^\prime$ be the $T$-dimensional time series vector of the $i$th observable, in vector notation (ref) also reads \[ \bm x_i = \alpha_i\bm\iota_T+\bm F\bm\lambda_i + \bm\zeta_i, \quad i=1,\ldots, N, \] where $\bm\zeta_i:=(\xi_{i1}\cdots \xi_{iT})^\prime$ is the corresponding $T$-dimensional time series vector of the $i$th idiosyncratic component, $\bm F:=(\mathbf F_1\cdots\mathbf F_T)^\prime$ is the $T\times r$ matrix of factors, and $\bm\iota_T$ is a $T$-dimensional column vector of ones.

Finally, given that $\bm X=(\mathbf x_1\cdots\mathbf x_T)^\prime$, in matrix notation (ref) reads

where $\bm \Xi:=(\bm\xi_1\cdots\bm \xi_T)^\prime$ is the $T\times N$ matrix of idiosyncratic components.

Structural assumptions

In this section we provide the minimal set of assumptions required to derive the results given in the rest of the paper. Notice that, although these assumptions do not coincide with those made by bai03, baili12, and baili16, they nest all assumptions made in those papers. More precisely, they nest all necessary assumptions in those papers, since, in fact, some of the assumptions made therein turn out to be redundant. In the following, we compare our assumptions only with baili12,baili16, while for a comparison with bai03 we refer to MBPCA.

We characterize the common component by means of the following assumption.

Parts (a) and (b), are also assumed in baili16, and are similar to what assumed by baili12. They imply pervasiveness of the factors, i.e., all factors have a non-negligible effect on all observables. Indeed, the loadings matrix has asymptotically maximum column rank $r$ (part (a)), the factors have a finite full-rank covariance matrix (part (b)), and, because of the upper bound on $\Vert\bm\lambda_i\Vert$ in part (a), for any given $N\in\mathbb N$, the contribution of each factor to every series is finite.

Part (c) is analogous to what assumed by baili12,baili16 in the case of deterministic factors. It is very general, it simply says that the sample covariance matrix of the factors is a consistent estimator of its population counterpart $\bm\Gamma^F$. It is satisfied by linear processes having innovations with finite 4th order moments, or strong mixing processes with finite 4th order moments, or, more in general, by processes having summable 4th order cross-cumulants, which is a necessary and sufficient condition for consistency of the sample covariance matrix (hannan).

Part (d) implies the existence of a finite of factors for all $N,T\in\mathbb N$. In particular, the numbers of common factors, $r$, is identified only for $N\to\infty$. Here $N_0$ is the minimum number of series we need to be able to identify $r$ so that it must be such that $N_0\ge r$. Without loss of generality hereafter when we say “for all $N\in\mathbb N$” we always mean that $N>N_0$ so that $r$ can be identified. In practice, we must always work with $N$ such that $r<\min(N,T)$.

Last, part (e), assumes a $\alpha_i$ to be a fixed idiosyncratic effect. Sometimes (ref) is also written as baili12 $x_{it}=\alpha_i+\bm\lambda_i^\prime \bar{\mathbf F}+\bm\lambda_i^\prime(\mathbf F_t-\bar{\mathbf F})+\mathbf \xi_{it} = \alpha_i^*+\bm\lambda_i^\prime \mathbf F_t^*+\xi_{it}$, with $\bar {\mathbf F}:=T^{-1}\sum_{t=1}^T \mathbf F_t$, $\mathbf F_t^* := \mathbf F_t-\bar{\mathbf F}$, and $\alpha_i^*:=\alpha_i+\bm\lambda_i^\prime \bar{\mathbf F}$ which is then a random idiosyncratic effect.

Now, because of Assumption (ref)(a), we have that the eigenvalues of $\bm\Sigma_\Lambda$, denoted as $\mu_j(\bm\Sigma_\Lambda)$, $j=1,\ldots, r$ are such that $m_\Lambda^2\le \mu_j(\bm\Sigma_\Lambda)\le M_\Lambda^2$. Similarly, because of Assumptions (ref)(b) and (ref)(c), we have that the eigenvalues of $\bm\Gamma^F$, denoted as $\mu_j(\bm\Gamma^F)$, $j=1,\ldots, r$ are such that $m_F\le \mu_j(\bm\Gamma^F)\le M_F$.

Hereafter, denote the covariance matrix of the common component as $$ \bm\Gamma^\chi:=\mathbb{E}[\bm\chi_t\bm\chi_t^\prime]=\bm\Lambda\bm\Gamma^F \bm\Lambda^\prime, $$ with $j$th largest eigenvalue $\mu_j^\chi$, $j=1,\ldots, r$. Then, it follows that (MBPCA) for all $j=1,\ldots,r$,

This means we consider only strong factors. For weak factors see, e.g., onatski2012asymptotics, uematsu2021inference, bai2021approximate, and freyaldenhoven2022factor, among many others.

To characterize the idiosyncratic component, we make the following assumptions.

Part (a), is a minimum requirement for having non-degenerate idiosyncratic components. Although, in principle, the assumption on the lower bound for $\sigma_i^2$ could be removed bartholomew2011latent, here it is maintained for simplicity, so that the log-likelihood is always well-defined. Moreover, from part (a), jointly with Assumption (ref)(b) by which $\mathbb{E}[\chi_{it}]=0$, we see that we are implicitly assuming that the observables are such that $\mathbb{E}[x_{it}]=\alpha_i$, $i=1,\ldots, N$.

Part (b) has a twofold purpose. First, it limits the degree of serial correlation of the idiosyncratic components. Second, it also limits the degree of cross-sectional correlation between idiosyncratic components. In particular, part (b) implies the conditions required by baili16 (see Lemma 1(i)-1(iii) in MBPCA), i.e.,

where $\rho$ and $M_\xi$ are defined in Assumption (ref)(b). Similar results are also in FLM13, where sparsity of the idiosyncratic covariance matrix is assumed.

Define the $N\times N$ idiosyncratic covariance matrix and the diagonal matrix of idiosyncratic variances as: $$ \bm\Gamma^\xi:=\mathbb{E}[\bm\xi_t\bm\xi_t^\prime]\quad \text{and}\quad \bm\Sigma^\xi:=\text{diag}(\sigma_1^2,\ldots \sigma_N^2), $$ respectively. Then, letting $\mu_1^\xi$ be the largest eigenvalue of $\bm\Gamma^\xi$, it follows that (see Lemma 1(v) in MBPCA)

where $M_\xi$ is defined in Assumption (ref)(b). Condition (ref) was originally assumed by chamberlainrothschild83 to control the amount of idiosyncratic correlation that we can allow for. Moreover, by setting $k=0$ in Assumption (ref)(b), it follows also that for all $i=1,\ldots, N$ and all $N\in\mathbb N$, $\sigma_i^2 \le M_\xi$. Thus, all idiosyncratic components have finite variance.

Parts (c) and (d) require finite 4th order moments as in baili12 and finite sums of 4th order cross-moments. Notice that this is weaker than the requirement of finite 8th order cross-moments assumed by baili16. Part (d) also implies that the sample covariance between $\{\xi_{it}\}$ and $\{\xi_{jt}\}$ is a consistent estimator of $\mathbb{E}[\xi_{it}\xi_{jt}]$ (hannan).

The basic set of assumption is completed by the following requirement.

This assumption obviously implies that the common components are independent of the idiosyncratic components at all leads and lags and across all units. Notice that baili12,baili16 treat the factors as constant parameters, therefore, by construction, they have no relation with the idiosyncratic components. Finally, Assumption (ref) is compatible with a structural macroeconomic interpretation of factor models, according to which the factors driving the common component are independent of the idiosyncratic components representing measurement errors or local dynamics.

We conclude with an assumption of two Central Limit Theorems which are required to allow us to conduct inference on the loadings and the factors. It coincides with the assumption made by baili16.

Identification conditions

Let $\mu_i^x$, $i=1,\ldots, N$, be the $i$th largest eigenvalue of $\bm\Gamma^x:=\mathbb{E}[(\mathbf x_t-\bm\alpha)(\mathbf x_t-\bm\alpha)^\prime]$. Then, we have the following result.

The proof is omitted since it is a direct consequence of Weyl's inequality, combined with (ref) which follows from Assumption (ref), (ref) which follows from Assumption (ref), and since $\bm\Gamma^x=\bm\Gamma^\chi+\bm\Gamma^\xi$ which follows from Assumption (ref).

Lemma (ref) shows that the spectrum of the population covariance matrix of $\{\mathbf x_t\}$ has an eigen-gap which widens as $N$ grows. This is a necessary condition for identification of the model, and, indeed, virtually all existing methods for estimating the number of factors, $r$, are based on such property (baing02, ABC10, onatski10, ahnhorenstein13, trapani2018randomized). Notice, however, that even if an eigen-gap is displayed by an observed data set, this is not sufficient to say that the underlying data generating process is a factor model.

For QML estimation it is crucial that the parameters of the model and the factors are identified. However, it is well known that in general the factor loadings, as well as the factors, are not identified. Indeed, all the structures equivalent to (ref) can be obtained through an $r\times r$ invertible matrix $\mathbf R$, which defines the transformed quantities: $\mathbf F_t^o:=\mathbf R^{-1}\mathbf F_t$, $\bm{\Lambda}^{o}:=\bm{\Lambda}\mathbf R$, and $\bm\Gamma^{\xi o}:=\bm\Gamma^{\xi}$. Under such relationships, using only first- and second-moment information, as it is the case with gaussian QML estimation, the model specified by $\mathbf F_t^o$, $\bm{\Lambda}^{o}$, and $\bm\Gamma^{\xi o}$, is indistinguishable from the one given by $\mathbf F_t$, $\bm{\Lambda}$, and $\bm\Gamma^{\xi}$.

To obtain an identified model, we then need enough a priori structure to preclude any but the trivial transformation $\mathbf R=\mathbf I_r$. This can be achieved by imposing additional $r^2$ identifying constraints. In particular, we impose the following standard identifying conditions.

Under parts (a) and (b) the loadings and the factors are identified up to a sign baing13, hence, they are only locally identified. To achieve global identification, in part (c) we fix the sign indeterminacy.

Let now $\bm V$ be the $r\times r$ diagonal matrix of eigenvalues of $(\bm\Sigma_\Lambda)^{1/2}\bm\Gamma^F(\bm\Sigma_\Lambda)^{1/2}$ or equivalently of $(\bm\Gamma^F)^{1/2}\bm\Sigma_\Lambda(\bm\Gamma^F)^{1/2}$, or also of $\bm\Sigma_\Lambda\bm\Gamma^F$. Then, because of Assumption (ref):

These conditions are sometimes assumed directly in place of Assumption (ref), in which case the model is identified only as $N,T\to\infty$.

Moreover, from parts (a) and (b) it follows also that for all $N\in\mathbb N$, the $r\times r$ matrix $N^{-1}\bm\Lambda^\prime\bm\Lambda$ is the $r\times r$ diagonal matrix having as diagonal entries $N^{-1}\mu_1^\chi,\ldots, N^{-1}\mu_r^\chi$, i.e., the eigenvalues of $\bm\Gamma^\chi$ scaled by $N$, which therefore are assumed to be distinct. Note that this implies that $\overline C_j<\underline C_{j-1}$, $j=2,\ldots, r$, in (ref) and that also the diagonal elements of $\bm V$ are distinct. Requiring distinct eigenvalues is useful for identifying the eigenvectors of $\bm\Gamma^\chi$. Although, in principle, this requirement is needed only for PC estimation bai03, it is also needed to ensure consistency of the numerical maximization algorithms presented below and used for QML estimation, which are initialized with the PC estimators.

The identifying conditions in Assumption (ref) are similar to the set of constraints IC3 considered in baili12,baili16, where, however, part (a) is replaced by the requiring $N^{-1}{\bm\Lambda^\prime(\bm\Sigma^\xi)^{-1}\bm\Lambda}$ to be diagonal for all $N\in\mathbb N$.

Factor dynamics

Under Assumptions (ref)(b)-(ref)(d) the factors can be seen either as an $rT$-dimensional sequence of constant parameters (in which case in (ref)(b) the expectations would be replaced by averages and in (ref)(c) we would have an ordinary limit), or as a realization of an $r$-dimensional stochastic process $\{\mathbf F_t,\, t\in\mathbb Z\}$. The former case is commonly considered in classical factor analysis where observations are collected from random samples. The latter case is especially relevant in time series analysis. Nevertheless, if we believe the factors to be truly autocorrelated, then it would be desirable to explicitly model the dynamic evolution of the factors. This is the approach originally adopted by SS77 and quahsargent93 for dynamic factor analysis and then re-introduced by FGLR09 and DGRfilter,DGRqml.

In order to keep things simple, but without loss of generality, here we assume the VAR(1) specification:

We make the following assumption.

Assumption (ref) is standard in VAR literature. Part (a), requires the factors to have a causal autoregressive representation. If we let $\mathbf v_t:=\mathbf H\mathbf u_t$ and $\bm\Gamma^v:=\mathbb{E}[\mathbf v_t\mathbf v_t^\prime]$, then, under parts (b) and (c), $\mathbf H$ is identified as $\mathbf H=(\bm\Gamma^v)^{1/2}$. Part (e) assumes finite and summable 4th order cross-moments of the innovations $\{\mathbf u_t\}$. In part (f), we fix the initial conditions.

As a consequence, $\{\mathbf F_t\}$ admits the Wold representation $\mathbf F_t=\sum_{k=0}^{t-1}\mathbf C_k\mathbf v_{t-k}$, with coefficients $\mathbf C_k:=\mathbf A^k$ such that $\sum_{k=0}^\infty \Vert \mathbf C_k\Vert^2\le M_C$ for some finite positive real $M_C$ and with innovations process $\{\mathbf v_t\}$ which is a zero-mean white noise with finite and summable 4th order cross-moments, and with covariance matrix $\bm\Gamma^v$ finite and positive definite. In fact, parts (a) through (e) imply Assumption (ref)(b) and, jointly with part (f), they also imply Assumption (ref)(c), and both (ref)(b) and (ref)(c) become redundant in this setting.\footnote{Specifically, it follows that $\mathbb{E}[\mathbf F_t]=\sum_{k=0}^{t-1} \mathbf C_k\mathbb{E}[\mathbf v_{t-k}]=\mathbf 0_r$, $\bm\Gamma^F= \sum_{k=0}^{t-1}\mathbf C_k\mathbf H\mathbf H^\prime \mathbf C_k^\prime\le M_C^2 M_H^2$ and it is positive definite since $\mu_r(\bm\Gamma^F)\ge\mu_{r}(\mathbf H\mathbf H^\prime)>0$, $ \mathbb{E}[\left \Vert\mathbf F_t \right\Vert^4] \le r^{10}M_C^4M_H^4 K_u, $, and, last, $\mathbb{E}[\vert T^{-1/2}\sum_{t=1}^T F_{it}F_{jt}-\mathbb{E}[F_{it}F_{jt}]\vert^2]\le r^{10}M_C^4M_H^4 K_u$, which implies Assumption (ref)(c) by Chebychev's inequality. }

A review of the Principal Component estimators

In this section we briefly review Principal Components (PC) estimation of the static approximate factor model (ref). This is arguably the most popular estimation approach and it is a reference estimator also for QML estimation. The PC estimators are non-parametric and, in particular, their implementation does not require to specify any second order structure of the idiosyncratic components as long as Assumption (ref)(b) is satisfied. There are a few different, although asymptotically equivalent, definitions of PCs. We refer to MBPCA for a discussion of the various approaches.

The idea of PCs was introduced by pearson1901 as a way to find the the $r$ directions of best fit in an $n$-dimensional space. It was then proposed as a technique to estimate factor models by hotelling1933analysis, whose definition we follow in this paper. Specifically, the loadings are estimated as eigenvectors of the sample covariance matrix of the observables and the factors are estimated as the orthonormal PCs of $\bm X$ (see also lawleymaxwell71, mardia1979multivariate, jolliffe2002principal). This is also the definition adopted by FGLR09. In particular, the PC estimator of the loadings is defined as:

where $\widehat{\mathbf v}_{i}^{x\prime}$ is the $r$-dimensional $i$th row of the $N\times r$ matrix $\widehat{\mathbf V}^x$ having has columns the $r$ normalized eigenvectors of the $N\times N$ sample covariance matrix $T^{-1}(\bm X-\bm\iota_T\bar{\mathbf x}^\prime)^\prime(\bm X-\bm\iota_T\bar{\mathbf x}^\prime)$ corresponding to the $r$ largest eigenvalues collected in the diagonal matrix $\widehat{\mathbf M}^x$. Notice that, by construction, the $N\times r$ matrix of PC estimated loadings $\widehat{\bm\Lambda}^{\text{\tiny PC}}:=(\widehat{\bm\lambda}_1^{\text{\tiny PC}}\cdots \widehat{\bm\lambda}_N^{\text{\tiny PC}})^\prime$ is such that $N^{-1}\widehat{\bm\Lambda}^{\text{\tiny PC}\prime}\widehat{\bm\Lambda}^{\text{\tiny PC}}=N^{-1}\widehat{\mathbf M}^x$, which is diagonal. Hence, $\widehat{\bm\Lambda}^{\text{\tiny PC}}$ satisfies Assumption (ref)(a).

Under the identifying Assumptions (ref), we have that the PC estimator is asymptotically equivalent to the unfeasible Ordinary Least Squares (OLS) estimator we would get if the factors were observed:

In particular, it can be shown that (MBPCA)

from which we prove the following result.

Since the PC estimator is non-parametric the asymptotic covariance matrix $\bm{\mathcal V}_{i}^{\text{\tiny \upshape OLS}}$ has a sandwich form reflecting the fact that we do not take into account possible idiosyncratic serial correlations. The final expression of $\bm{\mathcal V}_{i}^{\text{\tiny \upshape OLS}}$ comes from the law of iterated expectations, and since factors and idiosyncratic components are independent by Assumption (ref).

The PC estimator of the factors is then obtained by linear projection of $\bm X$ onto $\widehat{\bm\Lambda}^{\text{\tiny PC}}$:

It follows that the $T\times r$ matrix of PC estimated factors $\widehat{\bm F}^{\text{\tiny PC}}:=(\widehat{\mathbf F}_1^{\text{\tiny PC}}\cdots \widehat{\mathbf F}_T^{\text{\tiny PC}})^\prime$ is such that $T^{-1}\widehat{\bm F}^{\text{\tiny PC}\prime}\widehat{\bm F}^{\text{\tiny PC}}=\mathbf I_r$. Hence $\widehat{\bm F}^{\text{\tiny PC}}$ satisfies Assumption (ref)(b).

Under the identifying Assumptions (ref), we have that the PC estimator is asymptotically equivalent to the unfeasible OLS estimator we would get if the loadings were observed:

In particular, it can be shown that (MBPCA)

from which we prove the following result.

Since the PC estimator is non-parametric the asymptotic covariance matrix $\bm{\mathcal W}_{t}^{\text{\tiny \upshape OLS}}$ has a sandwich form reflecting the fact that we do not take into account possible idiosyncratic cross-sectional correlations and heteroskedasticity. Notice that the asymptotic covariance could be time dependent if we had serially heteroskedastic idiosyncratic components, since it would then depend on $\bm\Gamma_t^\xi:=\mathbb{E}[\bm \xi_t\bm \xi_t^\prime]$.

The log-likelihood function

First, let us introduce the following notation. Define $\bm {\mathcal X}:=\text{vec}(\bm X^\prime)=(\mathbf x_{1}^\prime\cdots \mathbf x_{T}^\prime)^\prime$ and $\bm{\mathcal Z}:=\text{vec}(\bm \Xi^\prime)=(\bm\xi_{1}^\prime\cdots \bm\xi_{T}^\prime)^\prime$ as the $NT$-dimensional vectors of observations and idiosyncratic components, $\bm{\mathcal A}:=\text{vec}(\bm\alpha\bm\iota_T^\prime)$ as the $NT$-dimensional vector of constants, $\bm {\mathcal F}:=\text{vec}(\bm F^\prime)=(\mathbf F_{1}^\prime\cdots \mathbf F_T^\prime)^\prime$ as the $rT$-dimensional vector of factors, and $\bm {\mathfrak L}:=\mathbf I_T\otimes \bm\Lambda$ as the $NT\times rT$ block-diagonal matrix containing the loadings matrix replicated $T$ times.

Then, by taking the vectorized version of the transposed of (ref), we can write model (ref) as:\footnote{Recall that $\text{vec}(\bm A\bm B\bm C)= (\bm C^\prime \otimes \bm A) \text{vec}(\bm B)$.}

Let $\bm\Omega^x:=\mathbb{E}_{}[(\bm {\mathcal X}- \bm{\mathcal A})(\bm {\mathcal X}- \bm{\mathcal A})^\prime]$ and $\bm\Omega^\xi:=\mathbb{E}_{}[\bm{\mathcal Z}\bm{\mathcal Z}^\prime]$, be the $NT\times NT$ covariance matrices of the $NT$-dimensional vectors of data and idiosyncratic components, respectively. Let also $\bm\Omega^F:= \mathbb{E}_{}[\bm {\mathcal F}\bm {\mathcal F}^{\prime}]$ be the $rT\times rT$ covariance matrix of the $rT$-dimensional factor vector. Then, because of Assumption (ref), \[ \bm\Omega^x = \bm {\mathfrak L}\,\bm\Omega^F\bm {\mathfrak L}^\prime + \bm\Omega^\xi. \] The parameters of the model that need to be estimated are then $\bm\alpha$ and $$ \bm\varphi:=(\mathrm{vec}(\bm\Lambda)^\prime, \mathrm{vech}(\bm\Omega^{\xi})^\prime,\mathrm{vech}(\bm\Omega^{F})^\prime)^\prime. $$

Exact log-likelihood

The gaussian quasi-log-likelihood computed in a generic value of the parameters, denoted as $\underline{\bm{\alpha}}$ and $\underline{\bm{\varphi}}$, is given by

Let $\bar{\bm{\mathcal X}}=\text{vec}(\bar {\mathbf x}\bm\iota_T^\prime)$ with $\bar {\mathbf x}=T^{-1}\sum_{t=1}^T\mathbf x_t$. Then, it holds that \[ \left(\bm {\mathcal X}-\underline{\bm {\mathcal A}}\right)^\prime \left(\underline{\bm\Omega}^x\right)^{-1}\left(\bm {\mathcal X}-\underline{\bm {\mathcal A}}\right)= \left(\bm {\mathcal X}-\bar{\bm{\mathcal X}}\right)^\prime \left(\underline{\bm\Omega}^x\right)^{-1}\left(\bm {\mathcal X}-\bar{\bm{\mathcal X}}\right)+\left(\bar{\bm{\mathcal X}}-\underline{\bm {\mathcal A}}\right)^\prime \left(\underline{\bm\Omega}^x\right)^{-1}\left(\bar{\bm{\mathcal X}}-\underline{\bm {\mathcal A}}\right), \] which shows that (ref) is maximized when $\underline{\bm {\mathcal A}}=\bar{\bm{\mathcal X}}$, i.e., the QML estimator of $\bm\alpha$ is $\widehat{\bm\alpha}^{\text{\tiny QML}}:=\bar{\mathbf x}$.

The concentrated quasi-log-likelihood computed in a generic value of the parameters $\underline{\bm{\varphi}}$ is then given by:

In general, it might be very hard, if not impossible, to find the QML estimator of $\bm\varphi$, i.e., the maximizer of (ref). This is due to the very large number of parameters that need to be estimated. The QML estimator might not even exist unless other restrictions are imposed on the model. The main problem here is that, under our assumptions, we should estimate both ${\bm\Omega}^\xi$ and ${\bm\Omega}^F$, which are in general full block Toeplitz matrices containing autocovariance matrices up to lag $(T-1)$. Moreover, each of the $T$ blocks of ${\bm\Omega}^\xi$ is in general an $N\times N$ full matrix, so, without imposing further restrictions, ${\bm\Omega}^\xi$ alone contains already $NT(NT+1)/2$ parameters to be estimated, which is more than the $NT$ available data points.

Mis-specified log-likelihoods

Henceforth, we will be working with a mis-specified log-likelihood where we treat the idiosyncratic components as serially and cross-sectionally uncorrelated. Thus, we consider a mis-specified model for the idiosyncratic 2nd order structure described by: ${\bm\Omega}^\xi=\mathbf I_T\otimes {\bm\Sigma}^\xi$, which reduces the parameters in $\bm\Omega^\xi$ to be estimated to just the $N$ elements of the diagonal matrix $\bm\Sigma^\xi$. The log-likelihood (ref) is then replaced by:

In general, maximization of (ref) is still unfeasible because it requires estimating $\bm\Omega^F$ which contains other $rT(rT+1)/2$ parameters. In this paper we consider two approaches. First, we consider a further mis-specification by imposing zero autocorrelation in the factors, i.e., by setting ${\bm\Omega}^F=\mathbf I_T\otimes{\bm\Gamma}^F=\mathbf I_{rT}$, where we used also Assumption (ref)(b). This is the approach reviewed in Section (ref), which defines the QML estimator $\widehat{\bm\varphi}^{\text{\tiny QML,S}}$. Second, we consider a parametric linear model describing the factor dynamics, so that $\bm\Omega^F$ depends on the parameters characterizing such model. In particular, in the case of the VAR(1) specification given in (ref), the $(t,s)$-block of $\bm\Omega^F$ is the $r\times r$ lag-$(t-s)$ autocovariance matrix: $\bm\Gamma_{t-s}^F:=\mathbb{E}[\mathbf F_t\mathbf F_s^\prime]$ such that $\text{vec}(\bm\Gamma_{t-s}^F)=(\mathbf I_r\otimes \mathbf A)^{t-s} (\mathbf I_{r^2}- \mathbf A\otimes\mathbf A)^{-1}\text{vec}(\mathbf H\mathbf H^\prime)$. This reduces the number of parameters in $\bm\Omega^F$ to be estimated to $r^2+r(r+1)/2$. This is the approach reviewed in Section (ref), which defines the QML estimator $\widehat{\bm\varphi}^{\text{\tiny QML,D}}$.

Identification of the QML estimator

Hereafter, we define the parameter space as the set $\mathcal O\in\mathbb R^Q$ where $Q$ is the number of parameters to be estimated. It is intended that any generic parameter vector $\underline{\bm\varphi}\in\mathcal O$ has elements satisfying all the assumptions given in Section (ref). In particular, notice that, under those assumptions, the dimension of $\mathcal O$ grows with $N$.

Usually in QML estimation it is required to first prove the existence of the QML estimator. However, to this end we cannot rely on compactness of $\mathcal O$ since its dimension increases to infinity. Nevertheless, as shown in the next sections, in the present setting the parameter estimates maximizing the log-likelihood have an asymptotic linear representation, so we can easily bound the estimation errors directly by the distance between the estimated and the true parameter vector without the need of using any Taylor expansion for asymptotic analysis. This allows us to bypass the need to work with a compact parameter space. The only condition we need is for each element of the parameter vector $\bm\varphi$ to belong to a bounded set, and this is guaranteed by Assumptions (ref)(a), (ref)(a), (ref)(b), (ref)(a), and (ref)(b).

Finally, in order to have a well defined log-likelihood in correspondence of the QML estimator of $\bm\varphi$, it is common to make the following assumption.

This assumption is identical to the requirement imposed by baili16. In fact, it is a redundant condition which we assume only for the sake of simplicity. Indeed, it can be easily proved that in the present context Assumption (ref) is always satisfied GGM21.

QML estimation of an approximate static factor model

In this section we review the case in which the dynamics of the factors is not explicitly modeled, i.e., when we set $\bm\Omega^F=\mathbf I_{rT}$. So we consider an approximate static factor model described by the equation: \[ x_{it}=\alpha_i+\bm\lambda_i^\prime\mathbf F_t+\xi_{it}, \quad i=1,\ldots, N, \quad t=1,\ldots, T, \] together with Assumptions (ref) through (ref). Then, either $\{\mathbf F_t\}$ is effectively an uncorrelated process, so that it can also be considered as a deterministic sequence, or we are introducing a further mis-specification on top of the mis-specifications related to the idiosyncratic 2nd order structure introduced in the previous section.

For a static factor model the log-likelihood (ref) simplifies to

Notice that, up to a scaling by $T$, this is also known in the literature as the Wishart log-likelihood for the sample covariance $T^{-1}(\mathbf x_{t}-\bar{\mathbf x})(\mathbf x_{t}-\bar{\mathbf x})^\prime$ AR56. The parameters to be estimated are: ${\bm\varphi}=(\mathrm{vec}(\bm\Lambda)^\prime,\sigma^2_1,\ldots,\sigma^2_N )^\prime$, which is a $Q=(r+1)N$-dimensional vector. We shall denote the maximizer of $\ell_{0,\text{\tiny S}}(\bm {\mathcal X};\underline{\bm\varphi})$ as $\widehat{\bm\varphi}^{\text{\tiny QML,S}}$.

Loadings

Despite all introduced simplifications, no closed form solution exists for the elements of the vector of QML estimators, but there exist few numerical ways to compute it. Classical iterative approaches, dealing with cross-sections of small dimension $N$, were proposed by, e.g, joreskog69 and lawleymaxwell71, and extended to the large $N$ case by sundbergfeldmann16. Similarly, ng15 adapted the Newton-Raphson maximization algorithm to the large $N$ case by exploiting the properties of tri-diagonal blocked matrices. However, all these algorithms are heuristic and their relation to proper QML estimation is unclear. In this respect, an EM algorithm, thus in principle accomplishing QML estimation, was introduced by RT82, and considered in the large $N$ setting by baili12,baili16.

Under our assumption of an approximate factor model, baili16 prove that

Although this result is derived under the alternative identifying condition $N^{-1}\bm\Lambda^\prime(\bm\Sigma^\xi)^{-1}\bm\Lambda$ being diagonal (IC3 in the original papers), it is possible to show that it is applicable also under our identifying conditions. In particular, in MBQMLPCA it is shown that, under the assumptions of this paper, the following holds:

Therefore, (ref), together with (ref), implies (ref), and we immediately see that the QML, PC, and the unfeasible OLS estimators of the loadings are asymptotically equivalent, as long as $N\to\infty$. This is formalized in the following result, of which part (a) is an immediate consequence of (ref) and Proposition (ref), and part (b), which holds regardless of the imposed identification conditions, is proved by baili16.

There are three mis-specifications that we introduced in the log-likelihood (ref). First, we treat the factors as serially uncorrelated. This is, however, an innocuous mis-specification that has no effect on the asymptotic properties of the loadings. Indeed, the OLS estimator, which is asymptotically equivalent to the QML estimator, is not affected by the autocorrelation of the regressors, i.e., of the factors, as long as their covariance matrix is well defined as required by Assumptions (ref)(b) and (ref)(c).

Second, the QML and PC estimators of the loadings are not the most efficient estimators, since they both neglect the possible idiosyncratic serial correlations, hence the sandwich form of the asymptotic covariance matrix, which for the QML estimator is due to the use of a mis-specified log-likelihood, while for the PC estimator is due to its non-parametric nature. Clearly, if the idiosyncratic components were truly serially uncorrelated, then, we would have $\mathbb{E}[\bm\zeta_i\bm\zeta_i^\prime]=\sigma^2_i\mathbf I_T$, and the asymptotic covariance matrix would reduce to $\sigma_i^2 (\bm\Gamma^F)^{-1}=\sigma_i^2\mathbf I_r$ (because of Assumption (ref)(b)), which is the Gauss-Markov lower bound under serial homoskedasticity. This is the case studied by baili12, and we refer to Section (ref) for more details.

Third, as a consequence of the fact that we are estimating an approximate factor model using the mis-specified log-likelihood (ref) of an exact factor model, we would not get consistent estimators of the loadings if $N$ were fixed. Indeed, only in the limit $N\to\infty$ we can disentangle the common components, i.e., the loadings and the factors, from the idiosyncratic components (see Lemma (ref)). This is reflected in the asymptotic expansion (ref), where it is clear that if $N$ is fixed then the estimation error is non-vanishing. In other words, the mis-specification error, which we introduce by using a mis-specified log-likelihood, vanishes asymptotically only if $N\to\infty$. So, despite being fully parametric, the QML estimator does not suffer of the curse of dimensionality, but, in fact, it produces consistent estimates only in a high-dimensional setting, i.e., it enjoys a blessing of dimensionality.

Factors

In this section we consider estimation of the factors when their dynamics is not specified. If the factors are treated as as a sequence of constant parameters, then it is possible to rewrite the mis-specified log-likelihood (ref) as (see, e.g., AR56 and baili12):

It is immediate to see that for known parameters $\bm\varphi$, the estimator of the factors maximizing (ref), is the unfeasible Weighted Least Squares (WLS):

Once we have also the QML estimator of the parameters $\widehat{\bm\varphi}^{\text{\tiny QML,S}}$, maximizing the log-likelihood (ref), we can compute (ref), and obtain the feasible WLS:

which is nothing else but the classical “least-squares” estimator of the factors originally proposed by bartlett37 (see also thomson36, bartlett38, and lawleymaxwell71). In the large $N$ case, baili12 and baili16 prove that:

The following result follows.

If instead we treat the factors as random variables, but we do not model their dynamics, then their optimal (in mean-squared sense) linear estimator is the linear projection of the true factors onto the observed data. Given that we are considering a mis-specified exact static factor model, such linear projection can be approximated by:

Once we have also the QML estimator of the parameters $\widehat{\bm\varphi}^{\text{\tiny QML,S}}$ maximizing the log-likelihood (ref) of a static factor model, we obtain the feasible Linear Projection (LP) estimator

which is nothing else but the classical “regression” estimator of the factors proposed by thomson51 (see also scott66, and lawleymaxwell71). Clearly, the unfeasible LP estimator (ref) has always a smaller MSE than the unfeasible WLS estimator (ref). Nevertheless, since because of Proposition (ref), and Assumptions (ref)(a) and (ref)(a), $\Vert\widehat{\bm\Lambda}^{\text{\tiny QML,S}}\Vert = O_{\mathrm P}(\sqrt N)$ and $\Vert (\widehat{\bm\Sigma}^{\xi\text{\tiny QML,S}})^{-1}\Vert = O_{\mathrm P}(1)$, then, the feasible WLS and LP estimators coincide asymptotically baili12:

The following result follows.

The asymptotic covariance matrix of the estimated factors is the same for both the WLS and LP estimator. It has a sandwich form reflecting the mis-specification of the cross-sectional correlation structure of the idiosyncratic components. It differs from the asymptotic covariance of the PC estimator given in Proposition (ref), which is an OLS-type estimator. Notice that, in general, under our assumptions of an approximate factor model, we cannot say if the WLS/LP estimators are more or less efficient than the PC estimator, i.e., the matrix $\bm{\mathcal W}_{t}^{\text{\tiny \upshape OLS}}-\bm{\mathcal W}_{t}^{\text{\tiny \upshape WLS}}$ is neither positive nor negative definite. The only thing we can say is that neither the WLS/LP nor the PC estimator are the most efficient, the efficiency loss of the former being due to the use of a mis-specified log-likelihood, and the efficiency loss of the latter being due to its non-parametric nature.

Clearly, if the true model were an exact factor model then we would have $\mathbb{E}[\bm\xi_t\bm\xi_t^\prime]=\bm\Sigma^\xi$, and $\bm{\mathcal W}_{t}^{\text{\tiny \upshape WLS}}$ would reduce to $(\bm\Sigma_{\Lambda\xi\Lambda})^{-1}$, which is the Gauss-Markov lower bound under cross-sectional heteroskedasticity. In this case the WLS/LP estimators are always more efficient than the PC estimator, since the latter, being non-parametric does not even account explicitly for idiosyncratic cross-sectional heteroskedasticity.

It is clear from Propositions (ref) and (ref), that in the fixed $N$ case no consistency can be proved for the estimated factors. Indeed, for fixed $N$ there is no hope of estimating the whole $rT$-dimensional vector of factors $\bm F$ consistently, simply due to a lack of degrees of freedom. Moreover, for any given $t=1,\ldots,T$, the WLS/LP estimators of the factors in (ref) are the result of the aggregation of data along the cross-sectional dimension, i.e., they are the result of a cross-sectional regression. Hence, for consistency we must require $N\to\infty$.\footnote{Notice that although the PC estimator of the factors proposed by bai03 is just an eigenvector so no cross-sectional aggregation is required, we still require $N\to\infty$ to prove its consistency. Indeed, in this case we need to compute eigenvectors of a $T\times T$ covariance matrix, which are consistently estimating the true eigenvectors only if $N\to\infty$.} This is another manifestation of the blessing of dimensionality.

In fact, from (ref) it is clear that we must have also $T\to\infty$ in order to consistently estimate the factors. This is because the feasible WLS/LP estimators in (ref) and (ref) require the use of the QML estimator of the loadings and the idiosyncratic variances, which are consistent if and only if $T\to\infty$ (see Proposition (ref)).

QML estimation of an approximate dynamic factor model

In this section we review the case in which the dynamics of the factors is explicitly modeled, so we consider an approximate dynamic factor model described by the linear system:

together with Assumptions (ref) through (ref). In this case, $\{\mathbf F_t\}$ is an autocorrelated stochastic process and $\bm\Omega^F$ is a function of $\mathbf A$ and $\mathbf H$.

For a dynamic factor model the log-likelihood (ref) can be decomposed as:

where (recall that $\mathbf F_0=\mathbf 0_r$ because of Assumption (ref)(g))

with $\mathbb C\text{ov}_{\underline{\bm\varphi}}(\mathbf F_t|\bm{\mathcal X})=\mathbb{E}_{\underline{\bm\varphi}}[ (\mathbf F_t-\mathbb{E}_{\underline{\bm\varphi}}[\mathbf F_t|\bm{\mathcal X}])(\mathbf F_t-\mathbb{E}_{\underline{\bm\varphi}}[\mathbf F_t|\bm{\mathcal X}])^\prime|\bm{\mathcal X}]$. The parameters to be estimated are $\bm{\varphi}=(\bm\phi^\prime,\bm\theta^\prime)^\prime$ with $\bm\phi=(\mathrm{vec}(\bm\Lambda)^\prime, \sigma_1^2,\ldots,\sigma_N^2)^\prime$ and $\bm\theta=(\mathrm{vec}(\mathbf A)^\prime, \mathrm{vec}(\mathbf H)^\prime)^\prime$, so ${\bm\varphi}$ is a $Q=(r+1)N+2r^2$-dimensional vector. We shall denote the maximizer of $\ell_{0,\text{\tiny D}}(\bm {\mathcal X};\underline{\bm\varphi})$ as $\widehat{\bm\varphi}^{\text{\tiny QML,D}}$.

We study QML estimation of the approximate dynamic factor model under the following simplifying assumption.

Gaussianity is a standard assumption in this setting, see, e.g., shumwaystoffer82, watsonengle83, and DGRqml. It strengthens the moment conditions made in Assumptions (ref) and (ref). The non-Gaussian case is discussed in quahsargent93 and BLqml. In particular, from this assumption it follows that, for all $N,T\in\mathbb N$ (see Appendix (ref) for a proof):

Loadings

The classical approach to estimate a dynamic factor model by QML is based on the maximization of the prediction error log-likelihood which depends on the linear prediction of the factors, usually computed via the Kalman filter, which is defined in the next section. This approach is equivalent to maximizing the log-likelihood (ref) (hannan2012statistical). Although there is no analytic solution for the QML estimator of the parameters obtained in this way, numerical maximization is standard practice at least for low dimensional exact dynamic factor models, see, e.g., stockwatson89,stockwatson91 and harvey90. In high-dimensions JK15 propose to maximize the prediction error log-likelihood subject to a preliminary step in which the data are projected onto a lower dimensional space. It is not clear, however, what are the effects of this first step on the asymptotic properties of the final estimators.

In high-dimensions maximization of the log-likelihood (ref) is usually achieved by means of EM algorithm, which is an iterative procedure proposed by DLR77 as a general way to achieve QML estimation when dealing with incomplete information (see also sundberg74,sundberg76, and sundberg2019statistical, for EM and incomplete data in exponential families). The use of the EM algorithm for estimating small dimensional state space models dates back to shumwaystoffer82, watsonengle83, harveypeters90, and GZ96, while its use to estimate a high-dimensional approximate dynamic factor model was first suggested by quahsargent93. The asymptotic properties of the estimated factors were firstly derived by DGRqml, and then completed by BLqml who derived also the asymptotic properties of the estimated loadings.

In the present context, the EM algorithm has two main features: (i) it gives closed form solutions for the estimated parameters, and (ii) it runs in few iterations. In a nutshell, it works as follows. Consider a given iteration $k\ge 0$ of the EM algorithm and let us assume to have an estimate of the parameters $\widehat{\bm\varphi}^{(k)}$. Then, for each $k$ repeat two steps.

The algorithm is intialized and terminated as follows.

The rationale for the EM algorithm is the following. By taking the conditional expectations of the log-likelihood (ref) and computed using $\widehat{\bm\varphi}^{(k)}$, we get:

The QML estimator of $\bm{\varphi}$ is then a maximum of the right hand side of (ref). In fact, since, by definition of Kullback-Leibler divergence, $\mathcal H(\underline{\bm\varphi},\widehat{\bm\varphi}^{(k)})\le \mathcal H(\widehat{\bm\varphi}^{(k)},\widehat{\bm\varphi}^{(k)})$ for any pair $(\underline{\bm\varphi},\widehat{\bm\varphi}^{(k)})$ and any $k\ge0$ (see DLR77, and wu83), it is enough to maximize $\mathcal Q(\underline{\bm\varphi},\widehat{\bm\varphi}^{(k)})$. It follows that the EM algorithm defines a continuous path in the parameter space from the starting point to the stopping point along which the log-likelihood monotonically increases without leaping over valleys.

In general, the implementation of the EM steps might be non-trivial. However, things simplify considerably in the present setup. First, the full-information expected log-likelihood $\mathcal Q(\underline{\bm\varphi},\widehat{\bm\varphi}^{(k)})$ can be further decomposed as:

Then, it is clear from (ref) that, in order to compute the expected full information log-likelihood in the E-step, we just need to compute the conditional moments of the factors: $\mathbb{E}_{\widehat{\bm\varphi}^{(k)}}[\mathbf F_t| \bm{\mathcal X}]$, $\mathbb{E}_{\widehat{\bm\varphi}^{(k)}}[\mathbf F_t\mathbf F_t^\prime| \bm{\mathcal X}]$, and $\mathbb{E}_{\widehat{\bm\varphi}^{(k)}}[\mathbf F_t\mathbf F_{t-1}^\prime| \bm{\mathcal X}]$. And, in turn, these moments can be obtained by means of the Kalman smoother implemented using the estimated parameters $\widehat{\bm\varphi}^{k}$, as detailed in the next section.

The M-step has then a closed form solution. Indeed, from (ref) we immediately see that the final EM estimator of the parameters, $\bm\phi$, in the measurement equation (ref) is:

such that $\widehat{\bm\Lambda}^{\text{\tiny EM}}:=(\widehat{\bm\lambda}_1^{\text{\tiny EM}}\cdots \widehat{\bm\lambda}_N^{\text{\tiny EM}})^\prime$ is the estimated matrix of all loadings and $\widehat{\bm\Sigma}^{\xi\text{\tiny EM}}:=\text{diag}(\widehat{\sigma}_1^{2\,\text{\tiny EM}},\ldots, \widehat{\sigma}_N^{2\,\text{\tiny EM}})$ is the estimated diagonal matrix of idiosyncratic variances.

Notice that, since to compute the estimators in (ref) and (ref) we use the mis-specified expected conditional log-likelihood (ref) of an exact factor model, we reduce a potentially hard maximizations problem of estimating a high-dimensional parameter vector $\bm\phi$ into the maximizations of $N$ expected log-likelihoods, each depending on a finite $r+1$-dimensional parameter vector $(\text{vec}(\bm\lambda_i)^\prime,\sigma_i^2)^\prime$. This can be seen as an extreme form of regularization of the idiosyncratic covariance matrix which makes estimation of the large dimensional parameter vector $\bm\phi$ straightforward.

Closed form expressions are immediately derived also for the final estimators of the parameters, $\bm\theta$, in the state equation (ref):

For a log-likelihood belonging to the exponential family, as the Gaussian one, it is possible to prove that the EM algorithm produces a sequence of estimators $\{\widehat{\bm\lambda}_i^{(k)}\}$ which converges to a local maximum of the log-likelihood (ref) as $k\to\infty$, say $\widehat{\bm\lambda}_i^{(\infty)}$ wu83.

It is then common practice to consider the output of the EM algorithm $\widehat{\bm\lambda}_i^{\text{\tiny EM}}$ as equivalent to the QML estimator $\widehat{\bm\lambda}_i^{\text{\tiny QML,D}}$. This, however, is not always true, because there are two additional sources of errors: (i) the local maximum might not be the global maximum, and (ii) the EM algorithm runs for a finite number of iterations $k^*$.

The two aforementioned sources of error are addressed by BLqml who prove that $\widehat{\bm\lambda}_i^{\text{\tiny EM}}$ is asymptotically equivalent to the QML estimator $\widehat{\bm{\lambda}}_i^{\text{\tiny QML,D}}$ maximizing the full log-likelihood (ref), which, in turn, is asymptotically equivalent to the QML estimator $\widehat{\bm{\lambda}}_i^{\text{\tiny QML,S}}$ obtained by maximizing the log-likelihood of a static factor model (ref). Namely,

and (up to logarithmic terms)

In particular, (ref) is a consequence of the fact that we initialize the EM algorithm with the PC estimator which is fully identified under our Assumption (ref)(b), so the sequence $\{\widehat{\bm\lambda}_i^{(k)}\}$ actually converges to the global identified maximum $\widehat{\bm\lambda}_i^{\text{\tiny QML,D}}$ (see also ruud91, for a similar result in the case of one-to-one mapping from the factors to the data, corresponding to the case of no idiosyncratic component). Moreover, (ref) follows from the fact that the numerical error is proportional to the information loss between using the expected joint log-likelihood in (ref) instead of the full log-likelihood (ref), i.e., it depends on the fraction of missing information due to the fact that the factors are not observed (sundberg76, MR94, and MLT07). Last, (ref) follows from the fact that the exact factor model log-likelihood (ref) and the approximate factor model log-likelihood (ref) are both asymptotically equivalent to the log-likelihood of $\bm{\mathcal X}$ conditional on the factors, which is the same in both cases, i.e., $\ell_{0,\text{\tiny S}}(\bm {\mathcal X}|\bm{\mathcal F};\underline{\bm\varphi})=\ell_{0,\text{\tiny D}}(\bm {\mathcal X}|\bm{\mathcal F};\underline{\bm\varphi})$, so asymptotically their maxima should coincide and they both coincide with the unfeasible OLS estimator (see also Remark (ref)).

From (ref), (ref), and (ref), jointly with Proposition (ref), it follows that.

Up to logarithmic terms, the EM estimator of the loadings has the same asymptotic properties of the QML estimator of a static factor model. Therefore, the same comments made after Proposition (ref) about efficiency apply also in this case. Moreover, the EM estimator of the loadings has also the same asymptotic properties of the PC estimator. This is not surprising since the EM algorithm is initialized with the PC estimator, and, therefore, $\widehat{\bm\lambda}_i^{\text{\tiny EM}}$ can be seen as a one-step estimator, which, by construction, must also be consistent and as efficient as the initial estimator $\widehat{\bm\lambda}_i^{\text{\tiny PC}}$ (see, e.g., VDV2000, and LC06).

Factors

In a dynamic factor model, the factors are explicitly treated as autocorrelated random variables, thus, their optimal linear predictor is the linear projection of the true factors onto the all available data which is collected into the $NT$-dimensional vector $\bm {\mathcal X}$. For generic values of the parameters and any given $t=1,\ldots,T$, such prediction is given by the linear projection: $\text{Proj}_{\underline{\bm\varphi}}\left({\mathbf F}_t|\bm{\mathcal X}\right) =: \bm{\mathcal P}_{\underline{\bm\varphi},t}^\prime\left(\bm{\mathcal X}-\underline{\bm{\mathcal A}}\right)$, where $\bm{\mathcal P}_{\underline{\bm\varphi},t}$ is an $NT\times r$ matrix, which depends on $\underline{\bm\varphi}$, and such that it minimizes the mean-squared error (MSE):

By solving (ref) for all $t=1,\ldots, T$ and for known parameters $\bm\varphi$, we obtain:

Now, under Assumption (ref), it is clear that $\text{Proj}_{\underline{\bm\varphi}}\left(\bm{\mathcal F}|\bm{\mathcal X}\right)=\mathbb{E}_{\underline{\bm\varphi}} [\bm{\mathcal F}|\bm{\mathcal X}]$, so (ref) coincides with the optimal predictor of the factors. In this case, (ref) can also be derived by noticing that given the joint distribution of $\bm{\mathcal F}$ and $\bm{\mathcal X}$ in (ref), we have (the conditional covariance is the Schur complement of the unconditional covariance):

In practice, we know that we can always estimate $\bm{\mathcal A}$ with its QML estimator $\bar{\bm{\mathcal X}}$ and, as in the previous sections, we treat the idiosyncratic components as if they were serially and cross-sectionally uncorrelated. In this case, from (ref), we get the following estimator of the factors:

This is nothing else but the unfeasible estimator obtained by taking the inverse Fourier transform of the smoother originally proposed by wiener49 and kolm41. Now since $\{\mathbf F_t\}$ is an autocorrelated process, $\bm\Omega^F$ is not diagonal and at each point in time the estimator of the factors (ref) is a weighted average of all $T$ present, past, and future values of all $N$ time series. Hence, differently from the PC, WLS, and LP estimators, this estimator is obtained not only by cross-sectional aggregation, but also by temporal aggregation of the data, thus effectively taking into account all its 2nd order dependencies. Given that (ref) uses all available history of $\{\mathbf x_t\}$, it is also called a smoother estimator, rather than a filter which, at a given point in time $t$, would use only past and present information.

Direct implementation of (ref) is still unfeasible due to the presence of the $rT\times rT$ full matrix, ${\bm\Omega}^{F}$ which needs to be estimated and inverted.\footnote{Although in this paper we do not consider such possibility, it is also worth noticing that in the case of macroeconomic data we often have factors with a singular spectral density, i.e., of rank $q<r$ dagostinogiannone12. This implies $\mbox{rk}(\bm\Omega^F)<rT$, so that, in such case, the smoother estimator (ref) is not even defined.} There are two main solutions to this problem. First, we could ignore the autocorrelation of the factors and impose $\bm\Omega^F=\mathbf I_{rT}$, because of Assumption (ref)(b), and then it is easily seen that (ref) would become the LP estimator defined in (ref).

Second, in the spirit of this section on estimation of a dynamic factor model, we can use a parametric model for the dynamics of the factors, so that $\bm\Omega^F$ is parametrized accordingly. In this case, we can estimate the factors by means of the well-known iterative procedures proposed by kalman60, producing the Kalman filter and smoother estimators, together with their conditional second moments. For known parameters $\bm\varphi$, these are nothing else but the linear projections: ${\mathbf F}_{t|t}:= \mbox{Proj}_{\bm\varphi}(\mathbf F_t|\mathbf x_t,\ldots, \mathbf x_1)$ and ${\mathbf F}_{t|T}:= \mbox{Proj}_{\bm\varphi}(\mathbf F_t|\bm{\mathcal X})$, respectively (see also rauch63, and dejong89). Hence, under the chosen parametrization of $\bm\Omega^F$, the Kalman smoother coincides with the estimator given in (ref).

It is now sensible to compute the Kalman smoother by using the EM estimator of the parameters $\widehat{\bm\varphi}^{\text{\tiny EM}}$, which approximate the QML estimators, defined in (ref) through (ref). However, those estimators of the parameters depend again on the factors. To break this mutual dependence, we can iterate between the Kalman smoother and the EM algorithm to obtain joint estimates of the parameters and the factors. In particular, at each iteration $k\ge 0$ of the EM algorithm we can implement the Kalman smoother using the parameters $\bm\varphi^{(k)}$ estimated at the previous iteration, thus giving estimated factors, $\mathbf F_{t|T}^{(k)}$, their conditional MSE, $\mathbf P_{t|T}^{(k)}$, and lag-1 conditional MSE, $\mathbf C_{t,t-1|T}^{(k)}$ (see below for the explicit definitions), necessary to compute the expected log-likelihood in the E-step. By means of these moments we can compute new estimates of the parameters $\bm\varphi^{(k+1)}$ in the M-step, as detailed in (ref) through (ref).

In practice, since, because of Assumption (ref), the conditional distribution of the factors given the data is Gaussian, see (ref), then, the conditional mean coincides with the linear projection. Therefore, the sufficient statistics (moments of 1st and 2nd order) needed to compute the final estimators of the parameters given in (ref) through (ref), are given by:

Notice, however, that, since we treat the idiosyncratic components as cross-sectionally uncorrelated, then $\mathbf P_{t|T}^{(k^*)}$ is just an approximation of the true MSE of the Kalman smoother. In particular, since the Kalman smoother is unbiased (duncan1972, and harvey90), the true MSE coincides with the conditional covariance as given in (ref). Nevertheless, since, given $\widehat{\bm\varphi}^*$, both matrices are $O_{\mathrm P}(N^{-1})$ (BLqml) this approximation error is negligible.

The final Kalman filter and Kalman smoother estimators are obtained at the end of the EM algorithm, i.e., they are given by $\widehat{\mathbf F}_t^{\text{\tiny KF}}:= \mathbf F_{t|t}^{(k^*+1)}$ and $\widehat{\mathbf F}_t^{\text{\tiny KS}}:= \mathbf F_{t|T}^{(k^*+1)}$. In particular, under the VAR(1) model in (ref), these are defined by the following iterations:

where $\widehat{\mathbf P}_{t|t-1} = \widehat{\mathbf A}^{\text{\tiny EM}}\widehat{\mathbf P}_{t-1|t-1} \widehat{\mathbf A}^{\text{\tiny EM}\prime} + \widehat{\mathbf H}^{\text{\tiny EM}} \widehat{\mathbf H}^{\text{\tiny EM}\prime}$, with $\widehat{\mathbf P}_{t|t}$, being the filtered conditional MSEs, which evolves according to the iterations given in Appendix (ref). These iterations are initialized by $\widehat{\mathbf F}_{0}^{\text{\tiny KF}}=\mathbf 0_r$, $\widehat{\mathbf F}_{T+1}^{\text{\tiny KS}}=\widehat{\mathbf A}^{\text{\tiny EM}} \widehat{\mathbf F}_{T}^{\text{\tiny KF}}$, and $\widehat{\mathbf P}_{0|0} = c \mathbf I_r$, for some finite positive real $c$.

Under the present setting, BLqml show that the Kalman smoother and the Kalman filter are asymptotically equivalent

and also (up to logarithmic terms):

that is, the Kalman filter and the feasible WLS/LP estimators are asymptotically equivalent (see also baili16 where however it is also required $T/N^3\to 0$). From (ref), (ref), and Propositions (ref) or (ref), we get the following result.

Up to logarithmic terms, the Kalman smoother estimator, computed using the EM estimator of the parameters, has the same asymptotic properties of the feasible WLS/LP estimators. Therefore, the same comments made after Propositions (ref) and (ref) about efficiency apply also in this case.

Classical ML estimation

In this section, we review classical factor analysis for cross-sectional data, which is characterized by: (i) fixed $N$, (ii) uncorrelated Gaussian, i.e., independent, idiosyncratic components $\bm\Omega^\xi=\mathbf I_T\otimes \bm\Sigma^\xi$, and (iii) non-stochastic factors, i.e., $\bm\Omega^F=\mathbf I_{T}\otimes \bm\Gamma^F$. This is the setting considered mainly in psychometric applications (see lawleymaxwell71 and references therein).

Classical Maximum Likelihood (ML) estimation of the factor model can be carried out in two ways.

In principle, approach A is considered to be the appropriate one in the case stochastic factors, while approach B is considered to be the appropriate one in the case non-stochastic factors anderson2003. However, approach A can be used also with non-stochastic factors. Indeed, the log-likelihood (ref) does not depend on the factors, nor on their second moments if we also assume $\bm\Gamma^F=\mathbf I_r$. Hence, approach A is commonly considered the most sensible one.

Another reason to prefer approach A is that the log-likelihood (ref), used in approach B, might diverge to infinity under certain choices of the parameters. To be more precise, this issue is related to just one specific configuration, namely the case in which at least one idiosyncratic component has zero variance (AR56, anderson2003, BT11). This implies that the ML estimator of the idiosyncratic variances might not be defined. To solve this problem, lawley_1942 suggested to look for solutions satisfying just the first-order conditions, but, as shown by solari69, this does not lead to discovering a maximum, but just a stationary point of the log-likelihood. Notice that, if we instead follow approach A, then the log-likelihood (ref) can easily accommodate some idiosyncratic components with zero variance, see, e.g., bartholomew2011latent.

Nevertheless, approach B has the nice feature of allowing us to exploit the mutual dependence of factors and loadings, and, moreover, the log-likelihood (ref) has a more tractable form than the log-likelihood (ref).

Regardless of the approach considered, when $N$ is fixed ML estimation of factor models suffers of two main problems: (I) inference for the loadings is computationally challenging, and (II) the factors cannot be estimated consistently and their estimation might even pose an incidental parameter problem for the estimation of the other parameters of the model. We now consider both issues in detail and show how letting $N\to\infty$ helps solving both problems and would allow us to use the estimation approach in B.

Consider first, the estimation of the loadings. Under the assumption of an exact factor model, baili12 prove that the ML estimator of the loadings, maximizing (ref), is such that:

From this result it is clear that since now the factor model is truly exact, then, the ML estimator of the loadings is $\sqrt T$-consistent, regardless of $N$. However, while if $N\to\infty$, the asymptotic covariance is the one given in Proposition (ref), when $N$ is fixed the expression is much more complex since in this case the error coming from (ref) is also $O_p(T^{-1/2})$ thus contributes to the asymptotic distribution. Indeed, from the asymptotic expansion (ref) it is clear that if $N$ is fixed we still have $\sqrt T$ consistency, but an additional term on top of the OLS estimation error, must be included in the asymptotic distribution. The specific expression of the asymptotic covariance in the fixed $N$ case is given, for example, in AR56, AFP87, and AA88 (see also MBPCA). Since for fixed $N$ defining a consistent estimator of such covariance matrix is not easy, inference might become a computationally hard problem. Notice that the same problem exists even if we assumed the simplest possible case of an exact factor model with the unrealistic assumption of cross-sectionally homoskedastic idiosyncratic components, as in young40 and whittle52.

Second, the factors represent an additional $rT$-dimensional vector of parameters that need to be estimated. They can either be estimated after we have an ML estimator of the loadings and the idiosyncratic variances (approach A) or they can be estimated jointly with the other parameters (approach B). Now, when $N$ is fixed none of the two approaches allows to retrieve the factors consistently. The impossibility to retrieve the factors consistently when $N$ is fixed is precisely the reason why in classical factor analysis the factors are usually identified only indirectly through their associated loadings which, as we have seen above are $\sqrt T$-consistent. This is the spirit of confirmatory analysis joreskog69, an example of which is in mardia1979multivariate, where the factors are identified through a VARIMAX orthogonal rotation of the loadings, providing rotated loadings with a few large entries and as many near-zero entries as possible. Another example is in terada2014strong, where the loadings are identified by the reduced $K$-means clustering method.

Moreover, besides the aforementioned problem of zero idiosyncratic variances, approach B, although appealing, poses also an incidental parameter problem. Indeed, we need to simultaneously estimate $rT+(r+1)N$ parameters using $NT$ observations, and we might easily run out of degrees-of-freedom. This is especially true if $N$ is small and $T$ is large, while the larger $N$ gets, the less binding the issue is (see Remark (ref) for details).

Now, if we let $N\to\infty$, then we can prove that under approach A the factors are consistently estimated (see Propositions (ref) and (ref)). Furthermore, as shown in Section (ref) below, when $N\to\infty$ also the incidental parameter problem vanishes asymptotically and we could use the estimation approach B by just iterating between the two sets of estimates. Say we start at iteration $k=0$ with an estimate of the loadings and the idiosyncratic variances, e.g., given by the PC estimator. Then, at any iteration $k\ge 0$ we would have the feasible WLS and OLS estimators:

Clearly since we initialize the algorithm with a consistent estimator of the parameters, then at each iteration $k\ge 0$ both estimators can be shown to be consistent too as $N,T\to\infty$ (see also pz23), but no proof of the convergence of $\widehat{\bm \lambda}_i^{(k+1)}$ to the ML estimator is available so far. Approach B would then become feasible and, possibly, more preferable than approach A. Indeed, approach B would allow us to explicitly take into account the mutual dependence of factors and loadings and we would also have closed form expressions for all estimators.

The incidental parameter problem

In this section we formally show that if $N$ is fixed, then the score of the log-likelihood (ref) with respect to the loadings does not satisfy the orthogonality condition by neyman79. Hence, any estimator the loadings using an estimator of the factors, is going to be biased. More precisely, for any specific value of the factors, $\widetilde{\mathbf F}_t$, of the loadings, $\widetilde{\bm\Lambda}$, and of the idiosyncratic variances, $\widetilde{\bm\Sigma}^\xi$, consider the score with respect to the loadings $\underline{\bm\lambda}_i$, $i=1,\ldots, N$,

Clearly, in the true values of factors and parameters the first-order conditions are satisfied:

since we treat the factors as constant parameters (but this would hold even for stochastic factors uncorrelated with the idiosyncratic components).

We now ask ourselves what happens if we replace the true factors with an estimator, and let us choose the best possible estimator, which is the unfeasible WLS in (ref). To this end, consider the curve $\bm\gamma:[0,1]\to \mathbb R^r$ connecting $\mathbf F_t$ and ${\mathbf F}_t^{\text{\tiny WLS}}$, i.e, such that $\bm\gamma(\eta)= {\mathbf F}_t+\eta({\mathbf F}_t^{\text{\tiny WLS}}-\mathbf F_t)$ for $\eta\in[0,1]$, and consider the $r$-dimensional function $\mathbf D_{\eta}: \mathbb R^r\to\mathbb R^r$ describing the variation of the score along the curve $\bm\gamma$, so that for any $\eta\in[0,1]$ we have the $r\times 1$ differential vector \[ \mathbf D_{\eta}\left({\mathbf F}_t^{\text{\tiny WLS}}-\mathbf F_t\right):=\partial_\eta \left\{ \mathbb{E} \left[ \bm s_i\left(\mathbf x_t;{\mathbf F}_t+\eta({\mathbf F}_t^{\text{\tiny WLS}}-\mathbf F_t),{\bm\Lambda},{\bm\Sigma}^\xi\right) \right] \right\}. \] In order, to have unbiased estimators of the loadings the first-order conditions should be locally insensitive to the value of the factors, which here play the role of nuisance parameters. A necessary condition for this to happen is that the Gateaux derivative along $\bm\gamma$ should be such that (see also chernozhukov18):

If condition (ref) is satisfied, then we could estimate the loadings by plugging into the log-likelihood (ref) noisy estimates of the factors, and then maximize the concentrated log-likelihood $\ell(\bm {\mathcal X};{\bm{\mathcal F}}^{\text{\tiny WLS}},\underline{\bm\varphi})$, without strongly violating the moment conditions (ref). However, it is easy to see that for fixed $N$ condition (ref) is not satisfied. Indeed, by using (ref) and the expression of the WLS estimator (ref) into (ref), we get:

which for fixed $N$ is never zero, not even if the model is exact, in which case the summation will just made of one term when $j=i$. So even if we use the ideal unfeasible WLS estimator of the factors, we would get biased estimates for the loadings. However, we clearly see that in the limit $N\to\infty$, the expression in (ref) is $O_{\mathrm P}(N^{-1})$, because of Assumptions (ref)(a), (ref)(a), and (ref)(b), and condition (ref) is satisfied. So, as expected, the problem of incidental parameters vanishes asymptotically.

On the mis-specification error and the search for efficient estimators

Consider the two equivalent ways of writing the factor model (ref) when stacking all $NT$ observations in one $NT$ dimensional vector, i.e., by vectorzing the matrix representation (ref) or its transposed:

Let $\bm\Theta^\xi:=\mathbb{E}[\text{vec}(\bm \Xi)\text{vec}(\bm \Xi)^\prime]=\mathbb{E}[\bm{\mathfrak E}\bm{\mathfrak E}^\prime]$ and $\bm\Omega^\xi:=\mathbb{E}[\text{vec}(\bm \Xi^\prime)\text{vec}(\bm \Xi^\prime)^\prime]=\mathbb{E}[\bm{\mathcal E}\bm{\mathcal E}^\prime]$, which are $NT\times NT$ matrices having the same elements just ordered differently.

If the factors were observed, the true, i.e., not mis-specified, log-likelihood (ref) could be written also as:

Then, for known $\bm\Theta^{\xi}$ the QML estimator of the loadings, maximizing (ref), would be the unfeasible GLS estimator:

This would be the QML estimator of the loadings for an approximate factor model with known factors.

Likewise, for known loadings the true, i.e., not mis-specified, log-likelihood (ref) could be written also as:

Then, for known $\bm\Omega^{\xi}$ the QML estimator of the factors, maximizing (ref), would be the unfeasible GLS estimator:

This would be the QML estimator of the factors for an approximate factor model with known loadings.

Now, even if the factors or the loadings were known, the GLS estimators (ref) and (ref) would still require the unfeasible task to estimate and invert $\bm\Theta^\xi$ or $\bm\Omega^\xi$. As done in Section (ref), we can then consider a mis-specified version of the log-likleihoods (ref) or (ref), where we do not account for idiosyncratic correlations and we set $\bm\Theta^\xi=\bm\Sigma^\xi\otimes \mathbf I_T$ and $\bm\Omega^\xi=\mathbf I_T\otimes\bm\Sigma^\xi$. Then, it is easy to see that the maximizers of such mis-specified log-likelihoods are the unfeasible OLS estimator of the loadings and the unfeasible WLS estimator of the factors. Indeed, in this case (ref) would simplify as follows:

with components $\bm\lambda_i^{\text{\tiny OLS}}$ given in (ref). And, similarly, (ref) would simplify as follows:

with components $\mathbf F_t^{\text{\tiny WLS}}$ given in (ref).

As shown in Section (ref), these unfeasible OLS and WLS estimators are asymptotically equivalent to the QML and feasible WLS estimators studied in Propositions (ref) and (ref), respectively. So to understand what is the effect of the mis-specifications introduced in the log-likelihood (ref), it is enough to study how different the unfeasible OLS and WLS estimators in (ref) and (ref) are from the most efficient estimators which are the unfeasible GLS estimators in (ref) and (ref), respectively.

On average the full and the mis-specified idiosyncratic covariances are asymptotically equivalent since, for all $N,T\in\mathbb N$,

because of Assumption (ref)(b) and since $M_\xi$ and $\rho$ do not depend on $i,j,t,s$.\footnote{Recall the inequalities: $\left\Vert\bm\Theta^\xi\right\Vert_F\le \sqrt{NT}\left\Vert\bm\Theta^\xi\right\Vert\le \sqrt{NT}\left\Vert\bm\Theta^\xi\right\Vert_1$ and $\left\Vert\bm\Omega^\xi\right\Vert_F\le \sqrt{NT}\left\Vert\bm\Omega^\xi\right\Vert\le \sqrt{NT}\left\Vert\bm\Omega^\xi\right\Vert_1$.} Nevertheless, the difference between the two sets of estimators depends on the inverse idiosyncratic covariances, i.e., on

which, even by assuming that $\bm\Theta^\xi$ and $\bm\Omega^\xi$ are positive definite, are still finite for all $N$ and $T$, but do not vanish. Therefore, as expected, the GLS estimators are in principle better than the unfeasible QML and WLS estimators described in the previous section.

In light of this discussion, it is natural to consider extensions of QML estimation of the static approximate factor model where we model explicitly some of the serial and cross-sectional idiosyncratic correlations. This would allow us to compute new estimators of the loadings and the factors which are closer to the unfeasible GLS estimators and thus are more efficient than the OLS and WLS.

In particular, the second best that we can hope for are the following pseudo-GLS estimators. Consider, first, estimation of the loadings when we do not model either the idiosyncratic cross- or cross-autocorrelations, i.e., when imposing $\mathbb{E}[\bm\zeta_i\bm\zeta_j^\prime]=\mathbf 0_{T\times T}$ for $i\ne j$. So we maximize the log-likelihood (ref) when $\bm\Theta^\xi$ is block-diagonal with $N$ blocks $\bm\Delta_i^\xi:=\mathbb{E}[\bm\zeta_i\bm\zeta_i^\prime]$, $i=1,\ldots, N$, each being a $T\times T$ matrix, with entries $\mathbb{E}[\xi_{it}\xi_{is}]$, $t,s=1,\ldots, T$, hence, under serial homoskedasticity, $\bm\Delta_i^\xi$ is a Toeplitz matrix. The QML estimator of the loadings is then the pseudo-GLS estimator:

This estimator has an asymptotic covariance matrix which attains the Gauss-Markov lower bound: $\lim_{T\to\infty} T(\mathbb{E}[\bm F^\prime (\bm \Delta_i^\xi)^{-1}\bm F])^{-1}$.

However, even for known factors, computing (ref) is in general unfeasible, since it requires estimation and inversion of a large matrix $\bm\Delta_i^\xi$. A common approach is to assume a parametric expression for $\bm\Delta_i^\xi$. For example, let $\xi_{it}=\alpha_i\xi_{i,t-1}+\nu_{it}$ with $\vert \alpha_i\vert<1$ and $\nu_{it}\sim (0,\omega_i^2)$ and uncorrelated across $i$ and $t$. Then, both $\bm\Delta_i^\xi$ and its inverse depend only on $\alpha_i$ and $\omega_i^2$.\footnote{ Indeed, we have that the $(t,s)$ entry of $\bm\Delta_i^\xi$ is $\frac{\omega_i^2\alpha_i^{|t-s|}}{1-\alpha_i^2}$ and \[ \left(\bm\Delta_i^\xi\right)^{-1}\!\!=\frac{\omega_i^2}{(1-\alpha_i^2)^2} \left(

\right). \] } In practice, we could first estimate the loadings by QML, then the factors by WLS, and last the idiosyncratic components, from which, by fitting an AR(1), we could get estimates of $\alpha_i$ and $\omega_i^2$ and thus of $\bm\Delta_i^\xi$ and its inverse. By using all these estimates we can compute a feasible version of the pseudo-GLS estimator of the loadings in (ref). This is similar to the approach used in PC estimation by BT11. and, for known factors, it coincides with the classical feasible GLS estimator of the loadings, which is asymptotically equivalent to their QML estimator when considering the mis-specified log-likelihood (ref), but for the filtered data $(1-\alpha_iL) x_{it}$, $i=1,\ldots, N$, and, thus, with $\bm\Sigma^\xi$ replaced by a diagonal matrix with entries $\omega_i^2$, $i=1,\ldots, N$ CU49.

Similarly, in order to explicitly model the autocorrelation of the idiosyncratic components, the EM algorithm presented in Section (ref) can also be modified. Namely, in the M-step of any iteration $k\ge 0$, we could iterate between estimates of the loadings conditional on $\alpha_i$ and $\omega_i$, and estimates of $\alpha_i$ and $\omega_i$ conditional on the loadings, until convergence, and then move to iteration $k+1$ of the EM algorithm. This is the procedure adopted by reiswatson10 and it is an instance of an Expectation Conditional Maximization (ECM) algorithm MR93.

Second, consider estimation of the factors when we do not model either idiosyncratic auto- and cross-autocorrelations, i.e., when we impose $\mathbb{E}[\bm\xi_t\bm\xi_s^\prime]=\mathbf 0_{N\times N}$ for $t\ne s$. So we maximize the log-likelihood (ref) when $\bm\Omega^\xi=\mathbf I_T\otimes \bm\Gamma^\xi$, where $\bm\Gamma^\xi:=\mathbb{E}[\bm\xi_t\bm\xi_t^\prime]$ is a $N\times N$ matrix, with entries $\mathbb{E}[\xi_{it}\xi_{jt}]$, $i,j=1,\ldots, N$, hence, under serial homoskedasticity, it is independent of $t$. The QML estimator of the factors is then the pseudo-GLS estimator:

This estimator has an asymptotic covariance matrix which attains the Gauss-Markov lower bound: $\lim_{N\to\infty} N(\bm \Lambda^\prime (\bm \Gamma^\xi)^{-1}\bm \Lambda)^{-1}$.

However, even for known loadings computing (ref) is in general unfeasible, since it requires estimation and inversion of a large matrix $\bm\Gamma^\xi$. To this end, bailiao16 propose an EM algorithm to maximize the log-likelihood (ref) when replacing $\bm\Sigma^\xi$ with the full matrix $\bm\Gamma^\xi$ but subject to an $\ell_1$ penalty imposed on its off-diagonal entries. They then use the QML estimators of the loadings and the idiosyncratic covariance obtained in this way to compute a feasible version of the pseudo-GLS estimator of the factors in (ref). Similar approaches are also in wang2019penalized, and poignard2020statistical. Generalizations of the EM algorithm, e.g., by considering penalized M-steps in order to account for cross-sectional idiosyncratic correlations, are in principle possible too. This is the spirit of the approach by LM02 who assume a sparse VAR model for the idiosyncratic components, thus accounting also for cross-autocorrelation, and they embed into the EM algorithm a penalized M-step.

Nevertheless, in general, the pseudo-GLS estimators (ref) and (ref) are not necessarily more efficient than the QML and WLS estimators studied in Proposition (ref) and (ref), since, they do not address possible cross-autocorrelations, so they are still maximizers of mis-specified log-likelihoods. We might have gains in efficiency only if the cross-autocorrelations are small enough, i.e., if the $T\times T$ off-diagonal blocks of $\bm\Theta^\xi$ or the $N\times N$ off-diagonal blocks of $\bm\Omega^\xi$ are sparse enough. Moreover, the feasible pseudo-GLS estimators discussed above my suffer if we do not correctly specify the autoregressive order of the idiosyncratic components or we do not select correctly the degree of penalization when thresholding the idiosyncratic covariance matrix. For these reasons such estimators are hardly considered in empirical econometric applications.

Simulations

Throughout, we let $N\in\{20, 50, 100,200\}$, $T=100$, and $r=2$, and, for all $i=1,\ldots, N$ and $t=1,\ldots, T$, we simulate the data according to

where $\bm\ell_i$ and $\bm f_t$ are $r$-dimensional vectors. Specifically,

The parameters $\tau$ and $\delta$ control the degrees of cross-sectional and serial idiosyncratic correlation in the idiosyncratic components. The the noise-to-signal ratio for series $i$ is given by $\theta_i$.

Finally, in order for the simulated loadings and factors to satisfy Assumptions (ref)(a) and (ref)(b) we proceed as follows. Given the common components is generated as $\chi_{it}={\bm\ell}_{i}^\prime {\bm f}_t$, let $\bm\chi_t=(\chi_{1t}\cdots \chi_{Nt})^\prime$ and compute $\widetilde{\bm\Gamma}^\chi=\frac 1T\sum_{t=1}^T \bm\chi_t\bm\chi_t^\prime$. Collect its $r$ non-zero eigenvalues into the $r\times r$ diagonal matrix $\widetilde{\mathbf M}^{\chi}$ and the corresponding normalized eigenvectors as the columns of the $N\times r$ matrix $\widetilde{\mathbf V}^{\chi}$, with rows $\widetilde{\mathbf v}_i^{\chi\prime}:=(\widetilde{v}_{i1}^{\chi}\cdots \widetilde{v}_{ir}^{\chi})$, $i=1,\ldots,N$. The loadings and factors are then simulated as \[ \bm \lambda_i:=(\widetilde{\mathbf M}^{\chi})^{1/2}\widetilde{\bm{\mathcal S}}\widetilde{\mathbf v}_i^{\chi}\; \text{ and }\; \mathbf F_t := (\widetilde{\mathbf M}^{\chi})^{-1/2} \widetilde{\bm{\mathcal S}}\widetilde{\mathbf V}^{\chi\prime}\bm\chi_t, \] with $\widetilde{\bm{\mathcal S}}$ being an $r\times r$ diagonal matrix having entries $\pm 1$ according to $\widetilde{\bm{\mathcal S}}_{jj} :=\mathbb I(\widetilde{v}^{\chi}_{1j}> 0)-\mathbb I(\widetilde{ v}^{\chi}_{1j}\le 0)$, $j=1,\ldots, r$, so that Assumption (ref)(c) is also satisfied. It is shown in Appendix (ref) that such transformation is equivalent to applying a linear invertible transformation to $\bm \ell_i$ and $\mathbf F_t$.

We simulate the model described above $B=500$ times, and at each replication, $b=1,\ldots, B$, we estimate the loadings by means of: (i) unfeasible OLS using the true factors, (ii) PC as in MBPCA, (iii) QML implemented via the EM algorithm in baili12,baili16, and (iv) EM algorithm as in DGRqml and BLqml. Both the QML and EM estimators are obtained by using the PC estimator to initialize the iterations. We estimate the factors in the following ways: (i) unfeasible OLS using the true loadings, (ii) PC using the PC estimator of the loadings as in MBPCA, (iii) LP and WLS using the QML estimators of the loadings and idiosyncratic variances as in baili12,baili16, and (iv) Kalman smoother using the EM estimators of the loadings and idiosyncratic variances as in DGRqml and BLqml.

In Table (ref) we report the Mean-Squared-Error (MSE) for each column, $j=1,\ldots, r$, of the considered loadings estimators, averaged over the $B$ replications (with standard deviations in parenthesis).

An application on a large euro area dataset

We analyze a new macroeconomic dataset of the euro area (EA) of $N=116$ quarterly series observed in the period 2000:Q1-2019:Q4. Missing values are imputed using the EM algorithm by stockwatson02JASA. This data is available from BLEA.

After transforming data to stationarity, the criterion by baing02 suggests the presence of $r=4$ common factors, explaining on averaged 63% of total variance in the data.

Figure (ref) shows the common component of: GDP growth rate, unemployment rate, 3-months interest rate, and inflation rates (measured as the growth rate of the GDP deflator and of the Harmonized Index of Consumer Prices). For each variable we consider the estimates obtained by means of PC analysis as in bai03, QML plus WLS as in baili16, and EM algorithm plus Kalman smoother as in DGRqml. Given Propositions (ref), (ref), (ref), (ref), (ref), and (ref), each estimated common components satisfy the following CLTs, as $N,T\to\infty$,

The 95% confidence bands are computed accordingly, and by estimating $\bm{\mathcal V}_i$ using the classical HAC estimator andrews91 and $\bm{\mathcal W}_t$ using a similar Cross-Sectional-HAC estimator baing06.

Concluding remarks

In this paper, we have shown two main results. First, when estimating an approximate static or dynamic factor model by QML, possibly via the EM algorithm, and using a mis-specified log-likelihood for an exact factor model, the estimated loadings are asymptotically equivalent to the loadings estimated via PC analysis. Why should then we be using the QML estimator rather than the PC estimator of the loadings? The main reason put forward by the literature seems to be the possibility of easily imposing constraints on the loadings, e.g., incorporating prior beliefs.

Second, we have also shown that the factors estimated via WLS or via the Kalman filter/smoother, are all asymptotically equivalent, and these estimators should be always preferred over the PC estimators as they address the heteroskedasticity of the idiosyncratic components so they might even be more efficient (depending on the amount of neglected cross-sectional correlations). Furthermore, in a time series setting the Kalman filter/smoother should be preferred as it allows to easily deal with missing values and given its dynamic aggregation property it can handle very persistent or even non-stationary data without any modification.

Summing up, although the QML plus WLS estimation approach, proposed by baili16, is asymptotically equivalent to the EM plus Kalman smoother approach, proposed by DGRqml, the latter is likely to produce estimated factors which enjoy better final sample properties and is a more flexible approach for dealing with complex economic datasets. Relevant examples of the use of this estimation technique are: