Recent Developments on Factor Models and its Applications in Econometric Learning

Introduction

The recent decade has witnessed a blossom of developments on statistical learning theories and practice, embraced with the exciting progresses on large-scale optimizations and dimension reduction techniques. Factor models, as one of the central machinery on summarizing and extracting information from large scale datasets, have received much attention in this revolutionary era of data science, and many breakthrough methodologies and applications have been developed in this exciting area.

This paper makes a selective overview on the recent developments of the factor model and its applications on econometric learning. Our review focuses on the perspective of the low-rank structure of factor models, and draws particular attentions to estimating the model from the low-rank recovery point of view. A central focus in the progress of this literature is the understanding and recovering low-rank structures of high-dimensional models. Many new learning theories and methods have been developed, which have revolutionized the modern understanding of econometric modeling. Meanwhile, the low-rank structure is one of the key properties of factor models. While this structure has long been aware of by researchers, studying the factor model from the perspective of low-rank matrix recovery is relatively new, and has led to many exciting new discoveries and understanding.

The survey mainly consists of three parts: the first part is a review on new factor estimation based on modern techniques on recovering low-rank structures of high-dimensional models. The second part discusses statistical inferences of several factor-augmented models and applications in statistical learning models. The final part summarizes new developments dealing with unbalanced panels from the matrix completion perspective.

We concentrate on recent developments on methodologies and applications in econometric learning. For a more comprehensive account on this topic, see Chapters 9-11 of the book by Fan et al. (2020c). Meanwhile, several important topics are not covered in this survey, but have also generated extensive researches in the literature. Those include selecting the number of factors, weak factors, identification, continuous-time and time-varying models, nonstationarity and structural breaks, Bayesian methods, bootstrap factors, as well as more sophisticated panel data models. Several excellent reviews have been written with emphasis on these topics. For those reviews, we refer to Stock and Watson (2016) for dynamic factor models with applications on macroeconomics, to Bai and Wang (2016) for time series and panel data models, and to Gagliardini et al. (2019) for a recent review on conditional factor models with applications to finance. Another class of estimation is a hybrid of PCA-method and the state space approach, see Giannone et al. (2008) and Doz et al. (2011) for more discussions. In addition, the generalized dynamic factor model is another important strand of literature, where factors are often estimated using the dynamic principal components, the frequency domain analog of principal components, developed by Brillinger (1964). Forni et al. (2000, 2005) provided rates of convergence of the common component estimated by dynamic principal components. Finally, we refer to the following papers for more detailed developments, among others: Bai and Ng (2002); Ahn and Horenstein (2013); Onatski (2010); Li et al. (2017), Bai and Li (2012, 2016), Onatski (2012); Chudik et al. (2011), Cheng et al. (2016); Massacci (2017); Gagliardini et al. (2016); Goncalves and Perron (2018); Baltagi et al. (2017); Barigozzi et al. (2018), Aït-Sahalia and Xiu (2017); Chen et al. (2019a); Liao and Yang (2018); Li et al. (2019); Pelger (2019), Su and Wang (2017).

Spiked Incoherent Low-Rank Models

Modern high-dimensional factor models can be viewed as a type of spiked incoherent low-rank model, a broad class of models that have drawn active research in the recent decade. A spiked incoherent low-rank model typically refers to a large matrix $\boldsymbol{\Sigma}$ (either observable or not), having the following decomposition:

Such decomposition satisfies the following three properties:

The rank of $\mathbf{L}$ is either bounded or grows very slowly compared to its dimensions.

The nonzero singular values of $\mathbf{L}$ grow fast, while the largest singular value of $\mathbf{S}$ is either bounded or grows much slower.

(also known as “pervasiveness”) The left and right singular vectors of $\mathbf{L}$, corresponding to the nonzero singular values, should have diversified elements, which means, elements of the rescaled singular vectors should be uniformly bounded.

The low-rank structure achieves dimension reductions: suppose the matrix $\boldsymbol{\Sigma}$ is of $N\times N_{1}$ dimensions, while the rank of $\mathbf{L}$ is $r$. Then the low-rank structure reduces the dimension from $O(NN_{1})$ to $O(N+N_{1})r$; the latter is the magnitude of the number of parameters in $\mathbf{L}$. Meanwhile, the spikedness helps seperate $\mathbf{L}$ from $\mathbf{S}$ approximately, and ensures that the large “signals” concentrate on $\mathbf{L}$, the low rank component. Finally, the incoherence, a condition that excludes matrices being low-rank and sparse simultaneously, enables us to estimate well the singular eigenvectors.

We explain these three properties using the matrix form of factor models. Consider

where $\mathbf{f}_{t}$ is a $r$-dimensional vector of factors; $\mathbf{b}_{i}$ is the loading vector and $u_{it}$ is the idiosyncratic noise. Specifically, (2.1) applies to two decompositions of this model.

Factor Decomposition. The matrix form of the factor model gives

where $\mathbf{Y}$ and $\mathbf{U}$ are $N\times T$ matrices of $y_{it}$ and $u_{it}$; $\mathbf{B}$ is the $N\times r$ matrix of $\mathbf{b}_{i}$ while $\mathbf{F}$ is the $T\times r$ matrix of $\mathbf{f}_{t}.$ Then corresponding to the notation (2.1), $\boldsymbol{\Sigma}=\mathbf{Y}$, $\mathbf{L}=\mathbf{M}$ and $\mathbf{S}=\mathbf{U}.$ In this decomposition, $\boldsymbol{\Sigma}$ is observable. Apparently, $\mathbf{M}$ is a low-rank matrix with rank $r$. The nonzero singular values of $\mathbf{M}$, under the strong factor assumption, grows much faster than those of $\mathbf{U}$, which gives rise to the spikedness property. Now let $\boldsymbol{\xi}$ be the $N\times r$ matrix whose columns are the left singular vectors of $\mathbf{M}$, and let $\boldsymbol{\xi}_{i}^{\prime}$ denote its $i$ th row. Then under the assumption that the nonzero eigenvalues of $\mathbf{B}^{\prime}\mathbf{B}$ grow fast with $N$, for some constant $C>0$,

which gives rise to the incoherent singular vectors. The right singular vectors can be bounded similarly.

Covariance Decomposition. It is also well known from the factor model (2.2) that the covariance matrix of $\mathbf{y}_{t}=(y_{1t},\cdots,y_{Nt})^{\prime}$, denoted by $\boldsymbol{\Sigma}_{y}$, can be decomposed as follows:

where $\boldsymbol{\Sigma}_{u}$ denotes the covariance matrix of $\mathbf{u}_{t}$. The above decomposition is well known for portfolio allocations and risk managements, where the total volatility is decomposed into the systematic risk $\mathbf{L}$, plus the (sparse) idiosyncratic risk $\boldsymbol{\Sigma}_{u}$. It also leads to the spiked incoherent low-rank model, but $\boldsymbol{\Sigma}_{y}$ is unknown and needs to be estimated.

2 Estimation

There are two general approaches to estimating model (2.1): (i) Principal Components Analysis (PCA), and (ii) low-rank regularization. Here we present a general PCA estimation setting, and defer the discussion of low-rank regularization to Section 3.2. We shall assume rank$(\mathbf{L})=r$ to be known.

For any matrix $\mathbf{A}$, let $\mathbf{A}=\mathbf{U}_{A}\mathbf{D}_{A}\mathbf{V}^{\prime}_{A}$ denote the singular value decomposition (SVD) of $\mathbf{A}$. Define the singular value hard thresholding operator as

where $\bar{}\mathbf{D}_{R}$ is a diagonal matrix that keeps the top $R$ diagonal elements of $\mathbf{D}_{A}$ and replaces the remaining elements by zeros. So $H_{R}(\mathbf{A})$ is the best rank $R$ matrix approximation to $\mathbf{A}$.

Suppose an estimator of $\boldsymbol{\Sigma}$, denoted by $\widehat{}\boldsymbol{\Sigma}$, is available, satisfying

for some sequences $\eta_{N}$ and $c_{N}$. We use $\widehat{}\boldsymbol{\Sigma}$ as the input matrix, which can be the sample covariance matrix or its robustfied versions (Fan et al., 2019c). The goal is to estimate $\mathbf{L}$ in (2.1) and its $N\times r$ matrix of the left singular vectors, denoted by $\boldsymbol{\xi}$ (also let $\boldsymbol{\zeta}$ denote its right singular vectors). We use respectively $\widehat{}\mathbf{L}:=H_{R}(\widehat{}\boldsymbol{\Sigma})$ with $R=r$, which is the rank $r$ projection of $\widehat{\boldsymbol{\Sigma}}$, and the $N\times r$ matrix $\widehat{}\boldsymbol{\xi}$ whose columns are the left singular vectors of $\widehat{}\boldsymbol{\Sigma}$. The following theorem, adapted from Fan et al. (2018), provides deviation bounds of the estimators. To make the paper self-contained, we also provide a simpler proof with slightly different conditions.

Consider the general model (2.1) with bounded $r:=\text{rank}(\mathbf{L})$. Suppose that $\min_{2\leq i\leq r+1}|\lambda_{i-1}(\mathbf{L})-\lambda_{i}(\mathbf{L})|\asymp\max_{2\leq i\leq r+1}|\lambda_{i-1}(\mathbf{L})-\lambda_{i}(\mathbf{L})|:=g_{N}$ and $\eta_{N}+\|\mathbf{S}\|=o_{P}(g_{N})$. Then, under condition (2.6), we have (i)

(ii) If additionally, $\|\mathbf{S}\|_{\infty}+\|\mathbf{L}\|_{\infty}=O_{P}(1)$, $N_{1}c_{N}=o_{P}(g_{N})$, then

This theorem is relatively general, and is applicable to low-rank models that are not necessarily consequences from factor models. The proof relies on perturbation bounds for singular vectors/values, and the achieved rates are sharp. Result (i) is simple and gives asymptotic bounds under the operator norm. Result (ii) gives element-wise deviation bound for the singular vectors, which requires more dedicated technical arguments.

Estimation under Factor Models

We observe an $N\times T$ data matrix $\mathbf{Y}$, which can be decomposed as

where $\mathbf{B}$ is $N\times r$ factor loadings matrix, $\mathbf{F}$ is $T\times r$ factors matrix and $\mathbf{U}$ is $N\times T$ idiosyncratic errors, which are uncorrelated with $\mathbf{M}:=\mathbf{B}\mathbf{F}^{\prime}$. All the three parts $\mathbf{B}$, $\mathbf{F}$ and $\mathbf{U}$ are unobserved. The $t$ th column of this expression can be written as

Under the model’s specification, we have the covariance structure (2.4). One of the most widely used estimation methods for the factor model is principal components analysis (PCA). Define the sample covariance $\mathbf{S}_{y}=\frac{1}{T}\sum_{t=1}^{T}\mathbf{y}_{t}\mathbf{y}_{t}^{\prime}=\frac{1}{T}\mathbf{Y}\mathbf{Y}^{\prime}$. Let $\widehat{}\boldsymbol{\xi}_{j}$ be the $j$th eigenvector corresponding to the largest $j$ th eigenvalues of $\mathbf{S}_{y}$. The PCA estimates $\mathbf{B}$ by taking $\widehat{}\mathbf{B}=\sqrt{N}(\widehat{}\boldsymbol{\xi}_{1},\cdots,\widehat{}\boldsymbol{\xi}_{R})$, which estimates $\mathbf{B}$ up to a diagonal transformation. Given $\widehat{}\mathbf{B}$, the factors can be estimated via the least squares:

This also leads to the estimated low-rank component $\frac{1}{T}\widehat{}\mathbf{B}\widehat{}\mathbf{F}^{\prime}\widehat{}\mathbf{F}\widehat{}\mathbf{B}^{\prime}$ for $\mathbf{B}\operatorname{cov}(\mathbf{f}_{t})\mathbf{B}^{\prime}$.

PCA is equivalent to the singular value hard thresholding by taking the input matrix $\widehat{}\boldsymbol{\Sigma}=\mathbf{S}_{y}$. Then $\frac{1}{T}\widehat{}\mathbf{B}\widehat{}\mathbf{F}^{\prime}\widehat{}\mathbf{F}\widehat{}\mathbf{B}^{\prime}=H_{R}(\mathbf{S}_{y})$. One can then apply Theorem 2.1 to infer the rates of convergence of the PCA estimators, which were obtained by Stock and Watson (2002a). Bai (2003) proved the asymptotic normality of PCA estimators for the factors and loadings. Results with general input $\widehat{\boldsymbol{\Sigma}}$ can be found in Chapter 10 of Fan et al. (2020c).

1.2 Maximum Likelihood Estimations

Another popular method to estimate a factor model is the maximum likelihood (ML) method (see, e.g., Lawley and Maxwell (1971), Bai and Li (2012), Doz et al. (2012)). Under the independence and normality assumptions, the log-likelihood function based on $\mathbf{y}_{t}$ is, for some constant $C$,

The log-likelihood function is then maximized with respect to the matrix parameters $(\mathbf{B},\operatorname{cov}(\mathbf{f}_{t}),\operatorname{diag}(\boldsymbol{\Sigma}_{u}))$ under additional restrictions that $\boldsymbol{\Sigma}_{u}$ is diagonal (Bai and Li, 2012, 2016) or sparse with regularizations (Bai and Liao, 2016; Wang et al., 2019b). Recently Barigozzi and Luciani (2019) explicitly accounted for autocorrelations of the factors in the likelihood function.

The factors can be estimated by two methods, one of which is the projection method. Under the joint normality assumptions of $\mathbf{f}_{t}$ and $\mathbf{u}_{t}$, we have

This provides the basis of estimating factors. The other approach is the generalized least squares: for given $\mathbf{B}$ and $\boldsymbol{\Sigma}_{u}^{-1}$, the GLS estimator for $\mathbf{f}_{t}$ is

Replacing the unknown parameters with their ML estimators, one obtains two estimators for the latent factors. Under large-$N$ setup, the difference of the two methods (PCA and MLE) for estimating factors are asymptotically negligible.

2 Low rank estimation

Given the low-rank structure of $\mathbf{M}$ (sparsity in singular value of $\mathbf{M}$), we can estimate the model via solving the following penalized optimization:

for some tuning parameter $\nu>0$. The solution is $\widehat{}\mathbf{M}=S_{\nu}(\mathbf{Y})$, where $S_{\nu}(\cdot)$ is the singular value thresholding operator (Ma et al., 2011), defined as follows. Let $\mathbf{Y}=\mathbf{U}_{y}\mathbf{D}\mathbf{V}_{y}^{\prime}$ be its SVD. Then $S_{\nu}(\mathbf{Y}):=\mathbf{U}_{y}\mathbf{D}_{\nu}\mathbf{V}_{y}^{\prime},$ where $\mathbf{D}_{\nu}=\operatorname{diag}(\{D_{ii}-\nu\}_{+})$ with $D_{ii}$ being the diagonal entries of $\mathbf{D}$. So $S_{\nu}(\mathbf{Y})$ applies “soft-thresholding” on the singular values of $\mathbf{Y}$. One can additionally estimate the factors and loadings using the singular vectors.

We note that this method is closely related to the PC-estimator, except the soft-thresholding is replaced by hard-threshoding. Let $R$ denote the “working number of factors”, which is the number of principal components one takes when applying the PC-method. We note that the PC-estimator for $\mathbf{M}$ with $R$ factors is given by (see Section 2.2):

This estimator is the solution to the penalized least squares problem (3.2) except that the nuclear norm is replaced by $\sum_{i=1}^{\min\{N,T\}}p_{\nu}(\psi_{i}(\mathbf{M}))$, where $p_{\nu}(\theta)=\nu^{2}-(\nu-|\theta|)_{+}^{2}$ is the harding thresholding penalty and $\psi_{i}(\mathbf{M})$ is the $i^{th}$ singular value of $\mathbf{M}$.

Therefore the difference between (3.2) and PCA is more fundamentally about that of hard- and soft- thresholding. Despite of many good properties, the soft-thresholding estimator possesses shrinkage bias, while the hard-thresholding reduces the bias. As a matter of fact, the shrinkage bias is on the singular values, rather than on the singular vectors. Indeed, the singular vectors of the two estimators are the same, and equal to the top $R$ singular vectors of $\mathbf{Y}.$ An important implication is that the factor estimator building on $\widehat{}\mathbf{M}$ is numerically equivalent to the PC-estimators for the factors, which does not suffer from any shrinkage bias. A formal statement and proof of the unbiasedness of eigenvectors can be found in Fan et al. (2019b).

2.2 Low-rank plus sparse decomposition

Recall that $\boldsymbol{\Sigma}_{y}$ and $\boldsymbol{\Sigma}_{u}$ denote the $N\times N$ covariance matrices of $\mathbf{y}_{t}$ and $\mathbf{u}_{t}$ in model (3.1), and that we have the following decomposition

We now demonstrate that this decomposition also provides a nice structure for estimating the covariance components. A key assumption is conditionally sparsity, namely, $\boldsymbol{\Sigma}_{u}$ is sparse. While the definition of sparsity may differ in different contexts, here we mean

should not grow too fast as $N\to\infty.$ This requirement can be weakened to approximate sparsity. In addition, $\mathbf{L}$ is a low-rank matrix. Thus we can directly estimate the above covariance decomposition via solving the following penalized optimization:

The above optimization can be solved by alternating the estimation of $\mathbf{L}$ and $\boldsymbol{\Sigma}_{u}$, and closed form solutions are available in both iterations. Given $\boldsymbol{\Sigma}_{u}$, solving for $\mathbf{L}$ leads to the singular value soft-thresholding: $\widehat{\mathbf{L}}=S_{\nu_{1}}(\mathbf{S}_{y}-\boldsymbol{\Sigma}_{u})$, and given $\mathbf{L}$, solving for $\boldsymbol{\Sigma}_{u}$ leads to the element-wise soft-thresholding: $\widehat{\boldsymbol{\Sigma}}_{u}=\widetilde{S}_{\nu_{2}}(\mathbf{S}_{y}-\mathbf{L})$. While both iterations solve convex problems, standard convergence analysis can be applied to show that the iterative algorithm converges in polynomial time.

A key quantity is the restricted strong convexity (RSC) constant, which is defined as follows:

We then have the following theorem, adapted from Agarwal et al. (2012). To make the paper self-contained, we also provide a proof with slightly different conditions. See the online supplement.

Conditioning on events $4\|\mathbf{S}_{y}-\boldsymbol{\Sigma}_{y}\|\leq\nu_{1}$ and $4\|\mathbf{S}_{y}-\boldsymbol{\Sigma}_{y}\|_{\infty}\leq\nu_{2}$, there is $C>0$ that only depends on $\text{rank}(\mathbf{L})$, so that

The optimal tuning parameters can be set to satisfy $\nu_{1}\asymp\frac{N}{\sqrt{T}}$ and $\nu_{2}\asymp\sqrt{\frac{\log N}{T}}$, respectively, accounting for estimating errors under two matrix norms:

both can be shown to hold with high probability under weak serial dependence and sub-Gaussian conditions. In additionally, if $\kappa(\nu_{1},\nu_{2})$ is bounded away from zero, with the choice of tunings, the convergence rate in Theorem 3.1 is $O_{P}(1+\frac{J\log N}{N^{2}})\frac{1}{T}$, which is sufficient to guarantee the convergence of the estimated factors and loadings. We refer to Lemma 2 of Agarwal et al. (2012) for more refined lower bound of $\kappa(\nu_{1},\nu_{2})$.

The above problem is also called “robust PCA” (Candès et al., 2011). For recent advance and references, see Chen et al. (2020c) where factorization methods are also discussed.

3 Covariance estimation

Fan et al. (2013) proposed a nonparametric estimator of $\boldsymbol{\Sigma}_{y}$, named POET (Principal Orthogonal complEment Thresholding), when the factors are unobservable. It is basically an one-step solution to optimization (3.4) with initialization $\boldsymbol{\Sigma}_{u}=0$. To motivate the estimator, suppose $r=R$. Then, heuristically

Thus, one estimates $\mathbf{L}$ by $H_{R}(\mathbf{S}_{y})$ and sets $\mathbf{S}_{u}:=\mathbf{S}_{y}-H_{R}(\mathbf{S}_{y})$. To account for the sparsity assumption on $\boldsymbol{\Sigma}_{u}$, Fan et al. (2013) estimates $\boldsymbol{\Sigma}_{y}$ and $\boldsymbol{\Sigma}_{u}$ as

where $h(x,\lambda_{ij})$ denotes the element-wise thresholding operator with thresholding value $\lambda_{ij}$. Here, we emphasize element-dependent thresholding $\lambda_{ij}$ to adapt to varying scales of covariance. For correlation thresholding at level $\lambda$, we take $\lambda_{ij}=\lambda\sqrt{s_{u,ii}s_{u,jj}}$ with $s_{u,ii}$ a diagnonal element of $\mathbf{S}_{u}$(Fan et al., 2013); we can also take other form such as the adaptive thresholding in Cai and Liu (2011). In general, the thresholding function should satisfy: (i) $h(x,\lambda)=0$ if $|x|<\lambda$, (ii) $|h(x,\lambda)-x|\leq\lambda$. (iii) there are constants $a>0$ and $b>1$ such that $|h(x,\lambda)-x|\leq a\lambda^{2}$ if $|x|>b\lambda$.

Note that condition (iii) requires that the thresholding bias should be of higher order. It is not necessary for consistent estimations, but we recommend using nearly unbiased thresholding (Antoniadis and Fan, 2001) for inference applications. One such example is known as SCAD. As noted in Fan et al. (2015), the unbiased thresholding is required to avoid size distortions in a large class of high-dimensional testing problems involving a “plug-in” estimator of $\boldsymbol{\Sigma}_{u}$. In particular, this rules out the popular soft-thresholding function, which does not satisfy (iii) due to its first-order shrinkage bias.

4 Projected PCA

In empirical asset pricing, factor loadings are known to depend on individual-specific observables $\mathbf{X}_{i}$, which represent a set of time-invariant characteristics such as individual stocks’ size, momentum, and values. To incorporate the information carried by the observed characteristics, Connor and Linton (2007) and Connor et al. (2012) model explicitly the loading matrix as a function of covariates $\mathbf{X}_{i}$. Fan et al. (2016) extended the model to allowing components in factor loadings that are not explainable by characteristics:

Here $\mathbf{g}(\cdot)$ is a vector of nonparametric functions. With this model, they introduced an improved factor estimator, known as projected PCA.

The basic idea of projected PCA is to smooth the observations $\{y_{it}\}_{i=1}^{N}$ for each given $t$ against their associated covariates $\{\mathbf{X}_{i}\}_{i=1}^{N}$ (cross-sectional smoothing), and apply PCA to the smoothed data (fitted values). Let $\{\phi_{j}(\mathbf{x})\}_{j=1}^{J}$ be a set of basis functions. This can be either unstructured, such as kernel machines, or structured such as a basis for additive models (Fan et al., 2020c). Set $\phi(\mathbf{X}_{i})^{\prime}=(\phi_{1}(\mathbf{X}_{i}),\cdots.,\phi_{J}(\mathbf{X}_{i}))$ and $\Phi(\mathbf{X})=(\phi(\mathbf{X}_{1}),\cdots,\phi(\mathbf{X}_{N}))^{\prime}$, an $N\times J$ matrix. Then the projection matrix on characteristics can be taken as $\mathbf{P}=\Phi(\mathbf{X})(\Phi(\mathbf{X})^{\prime}\Phi(\mathbf{X}))^{-1}\Phi(\mathbf{X})^{\prime}.$ The projected data $\mathbf{P}\mathbf{Y}$ is the fitted value of regressing $\mathbf{Y}$ on to the basis functions.

With probability approaching one, all the eigenvalues of $\frac{1}{N}({\mathbf{P}}\mathbf{B})^{\prime}{\mathbf{P}}\mathbf{B}$ are bounded away from both zero and infinity as $N\to\infty$.

The above conditions require that the strengths of the loading matrix should remain strong after the projection. Condition (ii) implies that if we apply $\mathbf{P}$ to both sides of $\mathbf{Y}=\mathbf{B}\mathbf{F}^{\prime}+\mathbf{U}$, then

we infer that the columns of $\mathbf{F}$ are approximately the eigenvectors of the $\mathbf{Y}^{\prime}\mathbf{P}\mathbf{Y}$, scaled by a factor $\sqrt{T}$. This motivates estimating factors by using the top $R$ eigenvectors of $\mathbf{Y}^{\prime}\mathbf{P}\mathbf{Y}$.

Fan et al. (2016) derived the rates of convergence of the projected PCA method. A nice feature is that the consistency of latent factors is achieved even when the sample size $T$ is finite so long as $N$ goes to infinity. Intuitively, the idiosyncratic noise is removed from cross-sectional projections, which does not require a long time series.

Similarly, in many applications, while we do not know the latent factors $\mathbf{f}_{t}$, we do know that factors are related to some proxy variables $\mathbf{W}_{t}$. For example, the latent factors are unknown for equity markets, but they are related to Fama-French factors (Fama and French, 2015); latent factors for disaggregated macroeconomics time series are unknown, but they are related to aggregated ones (McCracken and Ng, 2016). Switching the roles of rows and columns, longitudinal regression of each series $\{y_{it}\}_{t=1}^{T}$ on $\{\mathbf{W}_{t}\}_{t=1}^{T}$ yields the projected data matrix, from which latent factors and loadings can be extracted similarly. See Fan et al. (2020a) for details on how latent factor learning is augmented by instruments $\{\mathbf{W}_{t}\}_{t=1}^{T}$.

5 Diversified projection

In this section, we continue denoting by $R$ as the number of factors we use, and by $r$ as the true number of factors. Fan and Liao (2020) proposed a simpler factor estimator that does not rely on eigenvectors, by using cross-sectional diversified projections (DP). Let $\mathbf{W}=(\mathbf{w}_{1},\cdots,\mathbf{w}_{R})$ be a given exogenous (or deterministic) $N\times R$ matrix, where each of its $R$ columns $\mathbf{w}_{k}$ is an $N\times 1$ vector of “diversified weights”, whose definition is to be clear below. We estimate $\mathbf{f}_{t}$ by simply taking

By substituting $\mathbf{y}_{t}=\mathbf{B}\mathbf{f}_{t}+\mathbf{u}_{t}$ into the definition, immediately we have

Thus $\widehat{}\mathbf{f}_{t}$ (consistently) estimates $\mathbf{f}_{t}$ up to an $R\times r$ affine transform $\mathbf{H}$, with the estimation error $\mathbf{e}_{t}:=\frac{1}{N}\mathbf{W}^{\prime}\mathbf{u}_{t}$. The assumption that $\mathbf{W}$ should be diversified ensures that as $N\to\infty$, $\mathbf{e}_{t}$ is “diversified away” (converging to zero in probability). More specifically, we impose the following assumption.

There is a constant $c>0$, so that as $N\to\infty$, (i) The $R\times R$ matrix $\frac{1}{N}\mathbf{W}^{\prime}\mathbf{W}$ satisfies $\lambda_{\min}(\frac{1}{N}\mathbf{W}^{\prime}\mathbf{W})>c.$ (ii) $\mathbf{W}$ is independent of $\{\mathbf{u}_{t}:t\leq T\}$. (iii) Suppose $R\geq r$, $\operatorname{rank}(\mathbf{H})=r$ and $\psi^{2}_{\min}(\mathbf{H})\gg\frac{1}{N}$, where $\psi_{\min}(\mathbf{H})$ denotes the minimum nonzero singular value of $\mathbf{H}=\frac{1}{N}\mathbf{W}^{\prime}\mathbf{B}$.

Conditions (i) and (ii) define the “diversified weights” $\mathbf{W}$. When $(u_{1t},\cdots,u_{Nt})$ are cross-sectionally weakly dependent, they ensure that $\mathbf{e}_{t}$ is diversified away. Condition (iii) of Assumption 3.2 is a key condition, which requires that $\mathbf{W}$ should not diversify away the factor components in the time series. Several choices of $\mathbf{W}$ can be recommended to satisfy this condition. For instance, if factor loadings satisfy (3.6), then fix $R$ components of sieve basis functions: $(\phi_{1}(\cdot),\cdots,\phi_{R}(\cdot))$, we can define

Alternatively, we can also use transformations of the initial observation $\mathbf{x}_{t}$ for $t=0$, which was considered by Juodis and Sarafidis (2020). If $\mathbf{y}_{0}$ is independent of $\{\mathbf{u}_{t}:t\geq 1\}$, we can apply $w_{i,k}=\phi_{k}(y_{i,0})$ . These weights are correlated with $\mathbf{B}$ through $\mathbf{y}_{0}=\mathbf{B}\mathbf{f}_{0}+\mathbf{u}_{0}$.

An important benefit of the DP is that it is robust to over-estimating the number of factors. Theoretical studies of factor models have been crucially depending on the assumption that the number of factors, $r$, should be consistently estimated. This usually requires strong conditions on the strength of factors and serial conditions. Recently, Barigozzi and Cho (2018) proposed a PCA-based method to estimate factors that are robust to over-estimated $r$. They provided rates of convergence of the estimated common components when $R\geq r$.

Fan and Liao (2020) applied DP to several inference problems in factor-augmented models, including the post-selection inference, high-dimensional covariance estimation, and factor specification tests. They formally justified the robustness to over-estimating the number of factors in these applications. In particular, DP admits $r=0$ but $R\geq 1$ as a special case. That is, the inference is still valid even if there are no common factors present, but factors are nevertheless estimated for insurance. In addition, Karabiyik et al. (2019) applied DP to the context of panel data models in the presence of common factors.

6 Factor estimators robust to heavy tails

To apply either the PCA or the MLE to estimate the model, we need an initial covariance estimator $\mathbf{S}_{y}$, whose application requires elements of $\mathbf{y}_{t}$ have sufficient moments. Some technical results of factor estimations even require sub-Gaussian conditions on data’s tail distributions. However, heavy tailed data are not uncommon in economic applications. For instance, about thirty percent of 131 disaggregated macroeconomic variables of Ludvigson and Ng (2016) have excess kurtosis greater than six, so their distributions are fatter than the t-distribution with degrees of freedom five. Indeed, heavy tails are a stylized feature of high-dimensional data, as it is unlikely that all variables have sub-Gaussian tails.

Because the presence of heavy-tailed data invalidates many conditions required for estimating factor models, the recent literature has proposed several methods that are robust to the tail distributions. Here we describe two of them: truncation and robust M-estimation.

In the high-dimensional setting, consider estimating multivariate means from an independent triangular array variables $y_{i1},...,y_{iT}$ with $\max_{i\leq N}\operatorname{Var}(y_{it})\leq\sigma^{2}$. Truncate the data

Catoni (2012) constructed a robust M-estimator that shares the same Gaussian concentration. Fan et al. (2017a, 2019c) used the adaptive Huber’s loss to define the mean estimator:

where $\tau_{i}$ is a growing sequence, and

(ii) The robust M-estimation approach: Suppose $\log N=o(T)$, and the truncation parameter is set to satisfy $\tau_{i}\asymp\sqrt{\frac{T}{\log N}}.$ Then there is $c>0$ which does not depend on any moments of $y_{it}$, or $(N,T)$, with probability at least $1-4N^{-3}$,

Based on the above robust covariance inputs, we can create factor estimators and derive their theoretical properties following the guidance of Section 2.2. See Chapter 10 of Fan et al. (2020c) for further generalizations.

7 Use of cross-covariance

contains valuable information about $\mathbf{B}$. This motivates to estimate loadings by applying PCA to aggregated $\{\boldsymbol{\Sigma}_{h}:h=1,\cdots\}$, and we studied by Lam and Yao (2012). A related idea has been extended to matrix-variate PCA (Wang et al., 2019a; Chen et al., 2020a). Fan and Zhong (2018) also provided a procedure to efficiently aggregate the cross-covariance information with the covariance information when $h=0$.

8 Which method to use?

Many references have documented the comparisons among various estimation methods. Westerlund and Urbain (2013) made a comparison between PCA and cross-sectional averages in the panel data setting. Meanwhile, the PCA and low-rank penalized regressions are practically very similar. So we do not distinguish their use in practice. In general, because of the simplicity for implementations and relatively weak required conditions, the PCA still seems to be the most widely used method in applied research. Meanwhile, robust covariance inputs can also be integrated with the surveyed low-rank recovery methods.

In addition, when either factors or loadings can be partially explained by observed characteristics, the projected PCA is recommended. This is particularly useful in asset pricing applications where the explanatory power of asset characteristics has been well documented in the literature.

Factor-Augmented Inference and Econometric Learning

Forecasting in a data-rich environment has been an important research topic in economics and finance. Typical examples include forecasts of the aggregate output or inflation rate using a large number of the categorized macroeconomic variables.

Stock and Watson (2002a); Bai and Ng (2006) considered factor-augmented regression model for $h$-step ahead forecast:

Here $\mathbf{w}_{t}$ in (4.1) is the observed predictors, which may include lagged dependent variables. Equation (4.2) is a high-dimensional factor model that includes a vector of latent factors $\mathbf{f}_{t}$. The forecast can be implemented by regressing $y_{t+h}$ onto $\mathbf{w}_{t}$ and estimated factors. The factor model (4.2) serves as an important dimension reduction tool.

Fan et al. (2017b) generalized (4.1) to the nonlinear model with multi-indices. Consider the following forecasting model:

where $h(\cdot)$ is an unknown link function, and $\boldsymbol{\varepsilon}_{t+1}$ is the error independent of $\mathbf{f}_{t}$ and $\mathbf{u}_{t}$. Vectors $\boldsymbol{\phi}_{1},\dots,\boldsymbol{\phi}_{R}$ are $r$-dimensional linear-indepencent prediction indices. In contrast with linear forecasting, the above model specifies that the predicting function is nonlinear and depends on multiple indices of extracted factors. If we specify $R<r$, further dimension reductions are achieved.

A prominent result related to model (4.3) is given by Li (1991), which shows that under some regularity conditions such as $\mathbf{f}_{t}$ is elliptically symmetric, we have

which is a nonparametric covariance estimator. The above sliced covariance estimator is based on the observable factors. If the factors are unknown, they are replaced by their estimators, which leads to the following sufficient forecasting algorithm based on the factor models.

Sufficient forecasting algorithm based on the factor models.

Estimate factors in model (4.2) for $t=1,\dots,T$;

Construct the covariance estimator as in (4.5) with $\widehat{}\mathbf{f}_{t}$ in place of $\mathbf{f}_{t}$;

Obtain $\widehat{}\boldsymbol{\phi}_{1},\widehat{}\boldsymbol{\phi}_{2},\dots,\widehat{}\boldsymbol{\phi}_{R}$ by the top $R$ eigenvectors of the covariance in Step 2;

Construct the predictive indices $\widehat{}\boldsymbol{\phi}_{1}^{\prime}\widehat{}\mathbf{f}_{t},\dots,\widehat{}\boldsymbol{\phi}_{R}^{\prime}\widehat{}\mathbf{f}_{t}$;

Nonparametrically estimate $h(\cdot)$ with indices from Step 4, and forecast $y_{t+1}$.

2 Factor-adjusted regularized model selection

Consider a high-dimensional regression model

where $\mathbf{g}_{t}$ is a treatment variable whose effect $\boldsymbol{\beta}$ is of the main interest. The model contains high-dimensional exogenous control variables $\mathbf{x}_{t}=(x_{1t},\cdots,x_{Nt})$ that determine both the outcome and treatment variables. Having many control variables creates challenges for statistical inferences, as such, we assume that $(\boldsymbol{\nu},\boldsymbol{\theta})$ are sparse vectors.

Control variables are often strongly correlated due to the presence of confounding factors

This invalidates conditions of using penalized regressions to directly select among $\mathbf{x}_{t}$. Instead, if we substitute (4.8) to (4.6), we reach a factor-adjusted regression model:

where $\boldsymbol{\alpha}_{g}^{\prime}=\boldsymbol{\theta}^{\prime}\mathbf{B}$, $\boldsymbol{\alpha}_{y}^{\prime}=\boldsymbol{\beta}\boldsymbol{\alpha}_{g}^{\prime}+\boldsymbol{\nu}^{\prime}\mathbf{B}$, and $\boldsymbol{\gamma}^{\prime}=\boldsymbol{\beta}\boldsymbol{\theta}^{\prime}+\boldsymbol{\nu}^{\prime}$. Here $(\boldsymbol{\alpha}_{y},\boldsymbol{\alpha}_{g},\boldsymbol{\beta})$ are low -dimensional coefficient vectors while $(\boldsymbol{\gamma},\boldsymbol{\theta})$ are high-dimensional sparse vectors. Importantly, the model contains high-dimensional latent controls $\mathbf{u}_{t}$, which are weakly dependent due to the nature of idiosyncratic noises. The use of $\mathbf{u}_{t}$ instead of $\mathbf{x}_{t}$ validates conditions for many high-dimensional variable selection methods.

Fan et al. (2020b) and Hansen and Liao (2018) showed that the penalized regression can be successfully applied to (4.9) to select components in $\mathbf{u}_{t}$, which are cross-sectionally weakly correlated. Motivated by Belloni et al. (2014), the algorithm can be summarized as follows. For notational simplicity, we focus on the univariate case $\dim(\boldsymbol{\beta})=1$.

Estimate $\boldsymbol{\beta}$ as follows.

Estimate $\{(\mathbf{f}_{t},\mathbf{u}_{t}):t\leq T\}$ from (4.8) to obtain $\{(\widehat{}\mathbf{f}_{t},\widehat{}\mathbf{u}_{t}):t\leq T\}$.

Run penalized variable selections on $\widehat{}\mathbf{u}_{t}$:

Obtain residuals: $\widehat{}\boldsymbol{\varepsilon}_{y,t}=y_{t}-(\widehat{}\boldsymbol{\alpha}_{y}^{\prime}\widehat{}\mathbf{f}_{t}+\widehat{}\boldsymbol{\gamma}^{\prime}\widehat{}\mathbf{u}_{t}),$ and $\widehat{}\boldsymbol{\varepsilon}_{g,t}=\mathbf{g}_{t}-(\widehat{}\boldsymbol{\alpha}_{g}^{\prime}\widehat{}\mathbf{f}_{t}+\widehat{}\boldsymbol{\theta}^{\prime}\widehat{}\mathbf{u}_{t}).$

Estimate $\boldsymbol{\beta}$ by residual-regression: $\widehat{}\boldsymbol{\beta}=(\sum_{t=1}^{T}\widehat{}\boldsymbol{\varepsilon}_{g,t}^{2})^{-1}\sum_{t=1}^{T}\widehat{}\boldsymbol{\varepsilon}_{g,t}\widehat{}\boldsymbol{\varepsilon}_{y,t}.$

Note that $\boldsymbol{\gamma}:\to P_{\tau}(\boldsymbol{\gamma})$ is a sparse-induced penalty function with a tuning parameter $\tau$. When $\boldsymbol{\theta}$ and $\boldsymbol{\gamma}$ are sufficiently sparse, and the PC-estimator is used in step 1 with the correct selection of the number of factors, the above procedure is asymptotically valid:

where $\sigma_{g}^{2}$ and $\sigma_{\eta,g}^{2}$ are the asymptotic variances of $\boldsymbol{\varepsilon}_{g,t}$ and $\eta_{t}\boldsymbol{\varepsilon}_{g,t}$.

More recently, Fan and Liao (2020) showed that the assumption of correct selection of the number of factors can be relaxed if we use the diversified projection in step 1 instead, and (4.12) is still valid as long as we select $R\geq r$ factors (over selection). Importantly, this admits $r=0$, and $R\geq 1$ as a special case, i.e., there are no factors so that $\mathbf{x}_{t}=\mathbf{u}_{t}$ itself is cross-sectionally weakly dependent, but nevertheless we estimate $R\geq 1$ number of factors to run post-selection inference to alleviate the dependence among $\mathbf{x}_{t}$. This setting is empirically relevant as it allows to avoid pre-testing the presence of common factors for inference.

Figure 1, taken from Fan and Liao (2020), plots the histograms of the t-statistics based on estimated $\boldsymbol{\beta}$ over 200 simulations, superimposed with the standard normal density, where $R$ diversified projections are used to estimate factors in step 1. Here the weights are the initial transformations ($t=0$) so that the $i^{th}$ row of $\mathbf{W}$ is $(x_{it},x_{it}^{2},\cdots,x_{it}^{R})$ at $t=0$. The “double selection” is the algorithm used in Belloni et al. (2014) that directly selecting among $\mathbf{x}_{t}$, corresponding to the case $R=0$. The factor-augmented algorithm works well even if $r=0$; but when $r\geq 1$ factors are present, “double selection” leads to severely biased estimations.

Therefore as a practical guidance, we recommend that one should always run factor-augmented post-selection inference, with $R\geq 1$, to guard against confounding factors among the control variables.

3 Factor-adjusted robust multiple testing

Controlling the false discovery proportion (FDP) in large-scale hypothesis testing based on strongly dependent tests has been an important problem in many scientific discoveries across disciplines. See Fan et al. (2019a) and references therein, and Barras et al. (2010); Harvey et al. (2015); Harvey and Liu (2018); Giglio et al. (2020) for applications in empirical asset pricing.

Suppose we observe realizations of a random vector $\{\mathbf{y}_{t}=(y_{1t},\cdots,y_{Nt})^{\prime}\}_{t=1}^{T}$. Let $\boldsymbol{\alpha}=(\alpha_{1},\cdots,\alpha_{N})^{\prime}$ denote its mean vector. We are interested in testing individual hypotheses:

Let $p_{i}$ denote the $p$-value for testing $H_{0}^{i}$ based on a test statistic such as $t$-test, which rejects if $p_{i}<x$ given some critical value $x$. Define the number of false discoveries (rejections) and the total number of rejections as follows:

In large-scale multiple testing problems, researchers often aim to control the false discovery proportion (FDP) and the false discovery rate (FDR) defined by

The goal is to find the critical value $x$ so that FDR$(x)\leq\tau$ for a desired level $\tau$ (e.g., 0.10) or more relevantly FDP$(x)\leq\tau$ with high confidence. While $\mathcal{V}(x)$ is known, $\mathcal{F}(x)$ is not in practice. A general principle of finding $x$ proceeds as the following two steps.

Find $\bar{\mathcal{F}}(x)$ such that either it upper bounds $\mathcal{F}(x)$ for all $x\in(0,1)$, or it estimates $\mathcal{F}(x)$ uniformly well.

Set the critical value to $x^{*}=\sup\{x\in(0,1):\bar{\mathcal{F}}(x)\leq\tau\max\{\mathcal{V}(x),1\}\}$.

One of the most popular procedures, proposed by Benjamini and Hochberg (1995), proceeds as follows. Denote $p_{(1)}\leq\cdots\leq p_{(N)}$ as the sorted p-values for the individual tests. Then the critical value is set to

This method fits into Algorithm 4.3 with $\bar{\mathcal{F}}(x)=Nx$, which is an asymptotic upper bound for $\mathcal{F}(x)$ when the individual p-values are independent. One of the limitations of this upper bound is that it is too conservative if the number of true negatives is small compared to $N$. More fundamentally, it requires the test statistics be weakly dependent, a topic we shall discuss in more detail next. Other methods, such as Storey (2002); Fan et al. (2012), etc., aim to directly estimate $\mathcal{F}(x)$ in step 1 in the presence of strong dependence among test statistics, and are also adaptive to the unknown number of true negatives.

In addition, instead of Algorithm 4.3, Romano and Wolf (2007); Romano et al. (2008) provided alternative procedures for FDR control.

3.2 Removing dependence by factor adjustments

The key to the success of FDR control is that the individual test statistic should be either weakly dependent or independent. This makes the FDR and FDP approximately the same and easier to control. On the other hand, suppose the cross-sectional dependence of $\mathbf{y}_{t}$ is generated from a latent factor model:

To illustrate consequences of omitting adjusting latent factors as well as the effectiveness of the use of the factor-adjusted method (to be detailed below), let us consider a numerical example of a single factor model, where elements of $\mathbf{u}_{t}$, $\mathbf{f}_{t}$ and $\mathbf{B}_{t}$ are generated from the standard normal distribution. We take the true means to be $\alpha_{i}=0.6$ for $1\leq i\leq N/4$ and 0 otherwise, and compare two estimated $\alpha_{i}$: 1) the sample means of $\mathbf{y}_{t}$, without using factor adjustments; 2) the factor-adjusted estimator based on PCA. We apply the method of Benjamini and Hochberg (1995) for multiple testing, setting $\tau=0.05$.

The top panels of Figure 2 plot the histograms, from a single simulation, of the estimators for $\alpha_{i}$, corresponding to those that satisfy the null hypotheses $\alpha_{i}=0$ and those that satisfy the alternatives $\alpha_{i}=0.6$. Clearly, there is a large overlap (on the upper left panel) between sample means from the null and alternative, making tests based on sample means difficult to distinguish the alternatives from the nulls. In contrast, the PCA-based estimator can easily separate the nulls and alternatives, as shown on the upper right panel in Figure 2.

The middle two panels of Figure 2 plot the histograms of the true FDP over 1000 simulations based on the two estimators. It is evident that the distribution of the FDP corresponding to the factor-adjusted estimator concentrates around the nominal level. In contrast, the one based on the sample mean has a noticeable long tail as well as a larger mean and variance, which demonstrate the challenge to control FPD in presence of common factors, as explained above.

Finally, omitting confounding factors would lead to larger standard errors and conservative inference. The bottom two panels in Figure 2 plot the standard errors of individual estimated alphas and the sorted p values for the two estimation methods. The sample-mean estimator has much fewer sorted p-values below the B-H threshold line (i.e., fewer rejections), compared to the factor-adjusted estimator.

Hence it is recommended to estimate and remove the latent factors before applying standard FDR control algorithms.

3.3 Identifying skilled hedge funds

Giglio et al. (2020) studied the problem of identifying hedge funds that are able to produce positive alphas (i.e., have “skill”), among thousands of existing funds. They considered a linear pricing model, where hedge fund returns are:

In the model $\mathbf{f}_{t}$ contains both observable and latent factors. The model allows nontradable observable factors and $\boldsymbol{\lambda}$ is the vector of factor risk premia.

At a broad level, their methodology proceeds as the Fama-MacBeth regression integrated with the PCA to extract latent factors:

Estimating alphas in the presence of latent and nontradable factors.

Run fund-by-fund time series regressions to estimate fund exposures (betas) to observable factors.

Apply PCA to the residuals to recover the latent factors and betas.

Implement cross-sectional regressions like Fama-MacBeth to estimate the risk premia of the factors (including both observable and latent factors) and the alphas.

Because of many negative alphas from unskilled hund managers, the multiple testing problem should be properly formulated as one-sided hypotheses:

Hence rejecting $H_{0}^{i}$ indicates skilled fund manger $i$. On the other hand, the existence of potentially a very large number of negative alphas gives rise to the issue of power loss, only to add noises to the model. The loss of power associated with testing inequalities is well known as the problem of “deep in the null”, and is often seen in the econometric literature. To address this issue, Giglio et al. (2020) proposed to first screen off very bad funds, identified as:

where $c_{NT}>0$ is a slowly growing sequence to ensure sure screening (Fan and Lv, 2008): $P(\mathcal{I}\subseteq\mathcal{H}_{0})\to 1$. They recommended to apply FDR control algorithms on funds outside $\mathcal{I}$. Therefore, there are two ingredients that are recommended for identifying skilled fund managers via multiple testing: (1) adjust the effect of latent factors, and (2) remove the estimated alphas that are deep in the null. Both are playing essential roles of gaining good testing power.

4 Instrumental variable regression

The issue of endogeneity is often encountered in real data applications. Consider the following instrumental variable (IV) regression model

where $\mathbf{w}_{1t}$ is a $k_{1}$-dimensional vector of exogenous regressors and $\mathbf{w}_{2t}$ is a $k_{2}$-dimensional vector of endogenous regressors. Meanwhile, we have an $N$-dimensional IV $\mathbf{x}_{t}$ which admit a factor structure:

Below we introduce four estimators for $\boldsymbol{\beta}^{0}$, which differ on their choices of the instruments.

Use $\mathbf{f}_{t}$ as the instruments. Project $\mathbf{w}_{2t}$ on $\mathbf{f}_{t}$:

where $\boldsymbol{\phi}$ is a $k_{2}\times r$ matrix. We need $r\geq k_{2}$ for identification. Let $\mathbf{z}_{t}=(\mathbf{w}_{1t}^{\prime},\mathbf{f}_{t}^{\prime})^{\prime}$ be the set of instruments. As $\mathbf{f}_{t}$ is unobservable, we replace it with some factor estimator and apply the two stage least squares estimator $\widehat{}\boldsymbol{\beta}_{\mathbf{f}}$ with the feasible instruments.

Use $\mathbf{x}_{t}$ as the instruments. Project $\mathbf{w}_{2t}$ on $\mathbf{x}_{t}$:

where $\boldsymbol{\theta}$ is a $k_{2}\times N$ coefficient matrix. This projection motivates the use of $\mathbf{x}_{t}$ directly as a set of high-dimensional IV. Suppose that $\varepsilon_{t}$ is an i.i.d process, then the two-stage least squares estimator is efficient, and is given by

where $\mathbf{X}$ is $T\times N$ matrix of $\mathbf{x}_{t}$; $\mathbf{W}$ and $\mathbf{Y}$ are matrices of $\mathbf{w}_{t}$ and $y_{t}$. Note that $\widehat{}\boldsymbol{\Sigma}_{x}$ is the estimated covariance of $\mathbf{x}_{t}$, which can be constructed using factor-based covariance estimators as described in Section 3. It is interesting to compare the asymptotic behaviors of $\widehat{}\boldsymbol{\beta}_{\mathbf{f}}$ with $\widehat{}\boldsymbol{\beta}_{\mathbf{x}}$. Bai and Ng (2010) showed when $\mathbf{u}_{t}$ and $\mathbf{w}_{t}$ are uncorrelated, they have the same asymptotic variance, but $\widehat{}\boldsymbol{\beta}_{\mathbf{x}}$ has a $O(\frac{N}{T})$ bias term. So $\widehat{}\boldsymbol{\beta}_{\mathbf{x}}$ is consistent only if $N=o(T)$.

Use selected $\mathbf{x}_{t}$ as the instruments. We still consider the projection (4.14), but assume that rows of $\boldsymbol{\theta}$ are sparse vectors so that we can apply penalized regression to select among the components of $\mathbf{x}_{t}$:

where $P_{\tau}(\boldsymbol{\theta}_{j})$ is a sparse-induced penalty with tuning $\tau$. Let $\mathbf{x}_{t,\text{selec}}$ be the vector of selected components corresponding to nonzero components of $\{\widehat{}\boldsymbol{\theta}_{j}:j\leq\dim(\mathbf{w}_{2t})\}$. Belloni et al. (2012) used $(\mathbf{w}_{1t},\mathbf{x}_{t,\text{selec}})$ as the instruments to compute $\widehat{}\boldsymbol{\beta}_{\mathbf{x},\text{selec}}$, the two stage least squares estimator. This method however, would not work well in the presence of common factors. The strong dependence in $\mathbf{x}_{t}$ invalidates the variable selection procedure (4.15).

Use $\mathbf{f}_{t}$ and selected $\mathbf{u}_{t}$ as the instruments. We are not aware of any applications of this method in the IV literature, but it is still well motivated. Substitute the factor structure to (4.14), we obtain

where $\boldsymbol{\delta}=\boldsymbol{\theta}\mathbf{B}$. Hence we can carry out variable selections among $\mathbf{u}_{t}$:

Let $\widehat{}\mathbf{u}_{t,\text{selec}}$ be the vector of selected components corresponding to nonzero components of $\{\widehat{}\boldsymbol{\theta}_{j}:j\leq\dim(\mathbf{w}_{2t})\}$. We then use $(\mathbf{w}_{1t},\widehat{}\mathbf{f}_{t},\widehat{}\mathbf{u}_{t,\text{selec}})$ as the instruments to compute $\widehat{}\boldsymbol{\beta}_{\mathbf{f},\mathbf{u}}$, the two stage least squares estimator. This method is expected to work well because it marginalizes out the strong factors in $\mathbf{x}_{t}$, leaving remaining components $\mathbf{u}_{t}$ being weakly dependent.

Let us conduct a simple simulation to study the finite sample behaviors of the aforementioned four estimators. We consider a model $y_{t}=\mathbf{w}_{2t}^{\prime}\boldsymbol{\beta}^{0}+\varepsilon_{t}$, with a single endogenous regressor $\mathbf{w}_{2t}$ generated from (4.14) with $\mathbf{e}_{t}=\varepsilon_{t}/2$. Here $\boldsymbol{\theta}=(2,1,-1,0...,0)$ and $\mathbf{x}_{t}$ admits a two-factor structure. Variables $(\varepsilon_{t},\mathbf{f}_{t},\mathbf{B},\mathbf{u}_{t})$ are independent standard normal. Finally, variable selections are based on lasso with the oracle tuning parameter that controls the score of the least squares function. For instance, for problem (4.15) we set $P_{\tau}(\boldsymbol{\theta})=\tau\|\boldsymbol{\theta}\|_{1}$ with $\tau=2.2\|\frac{1}{T}\sum_{t}\mathbf{x}_{t}\mathbf{e}_{t}\|_{\infty}$.

Table 1 reports the bias and standard deviation of each estimator calculated from 1000 replications. First, using only estimated factors as the instruments ($\widehat{}\boldsymbol{\beta}_{\mathbf{f}}$) leads to the largest standard error. This is not surprising because it excludes the relevant information from $\mathbf{u}_{t}$ while the latter is correlated with $\mathbf{w}_{2t}$, so this method is less efficient. Secondly, using $\mathbf{x}_{t}$ as instruments without variable selection ($\widehat{}\boldsymbol{\beta}_{\mathbf{x}}$) has the smallest standard deviation, but is severely biased. Finally, the two instrumental selection based estimators ($\widehat{}\boldsymbol{\beta}_{\mathbf{x},\text{selec}}$ and $\widehat{}\boldsymbol{\beta}_{\mathbf{f},\mathbf{u}}$) perform favorably and similarly. But $\widehat{}\boldsymbol{\beta}_{\mathbf{x},\text{selec}}$ is not as stable, as it occasionally selects none of the instruments in our numerical experiments.

5 Boosting

Consider the following factor-augmented regression

where $\boldsymbol{\alpha}(\mathbf{L})=\boldsymbol{\alpha}_{0}+\boldsymbol{\alpha}_{1}\mathbf{L}+\dots+\boldsymbol{\alpha}_{p}\mathbf{L}^{p}$, $\boldsymbol{\gamma}(\mathbf{L})=\boldsymbol{\gamma}_{0}+\boldsymbol{\gamma}_{1}\mathbf{L}+\dots+\boldsymbol{\gamma}_{q}\mathbf{L}^{q}$ and $\boldsymbol{\beta}(\mathbf{L})=\boldsymbol{\beta}_{0}+\boldsymbol{\beta}_{1}\mathbf{L}+\dots+\boldsymbol{\beta}_{l}\mathbf{L}^{l}$, all are lag operator polynomials. Suppose that $\mathbf{w}_{t}$ is a $k$-dimensional vector and $\mathbf{f}_{t}$ is an $r$-dimensional vector. The above predictive regression has $n=1+(p+1)k+(q+1)+(l+1)r$ parameters. It is likely that partial parameters are zero. So model selection devices can be conducted to choose a parsimonious model. Here we briefly describe a model selection method, known as boosting, which was proposed to use by Bai and Ng (2009) in this context.

Initialize $\widehat{f}^{}(\cdot)$ an offset value. The default value is $\widehat{f}^{}(\cdot)\equiv\bar{y}$. Set $m=0$.

Increase $m$ by 1. Compute the residuals $e_{t}=y_{t}-\widehat{f}^{[m-1]}(\mathbf{z}_{t})$ for $t=1,2,\dots,T$.

Fit the residual vector $e_{1},\dots,e_{T}$ to $\mathbf{z}_{1},\dots,\mathbf{z}_{T}$ by the real-valued base procedure (e.g., regression):

Update $\widehat{f}^{[m]}(\cdot)=\widehat{f}^{[m-1]}(\cdot)+\nu\cdot\widehat{g}^{[m]}(\cdot)$, where $0<\nu\leq 1$ is a step-length factor.

One can apply the above $L_{2}$-Boosting to the factor-augmented predictive regression (4.16). As seen in Algorithm 4.5, one needs to specify the base procedure in step 3. Bai and Ng (2009) suggest two methods depending on the way to deal with lags, which leads to the component-wise $L_{2}$-Boosting and block-wise $L_{2}$-Boosting. In component-wise $L_{2}$-Boosting, one treats each lag of each variable as an independent predictor and the base procedure is a simple linear regression. Therefore, step 3 is given as follows.

Let $\mathbf{z}_{t,j}$ denote a typical regressor in the regressors pool with $j=1,2,\dots,n$. Regress the current residual $e_{t}$ (the residual in the $m$-th repetition) on each $\mathbf{z}_{t,j}$ to obtain the coefficient $\widehat{}\mathbf{b}_{j}$. Compute the sum of squared residuals, denoted by SSR($j$).

$\widehat{g}^{[m]}(\mathbf{x}_{t})=\mathbf{z}_{t,j_{m}}\widehat{}\mathbf{b}_{j_{m}}$ if $\mathbf{x}_{t}=\mathbf{z}_{t,j_{m}}$, and 0 otherwise.

Another way is to only differentiate the predictors in the current period and treat the predictor and its multiple lags as a block. This gives rise to the block-wise $L_{2}$-Boosting. The base procedure now is a multivariate regression with the regressors being one predictor and its lags. See Bai and Ng (2009) for details.

6 Threshold regression with mixed integer optimization

Threshold regressions have been used in economic applications to capture potential structural changes on regression coefficients. The early literature models the threshold effect using some observable scalar variable $q_{t}$ as in:

where $\mathbf{w}_{t}$ and $q_{t}$ are adapted to the filtration $\mathcal{F}_{t-1}$; $(\boldsymbol{\beta},\boldsymbol{\delta},\gamma)$ is a vector of unknown parameters, and $\varepsilon_{t}$ satisfies the conditional mean restriction. Hence when $q_{t}>\gamma$, the regression function becomes $\mathbf{w}_{t}^{\prime}(\boldsymbol{\beta}+\boldsymbol{\delta})$; when $q_{t}\leq\gamma$, it reduces to $\mathbf{w}_{t}^{\prime}\boldsymbol{\beta}$ (Chan, 1993; Hansen, 2000). In practice, it might be controversial to choose which observed variable plays the role of $q_{t}$. For example, if the two different regimes represent the status of two environments of the population, arguably it is difficult to assume that the change of the environment is governed by just a single variable.

Seo and Linton (2007) and Lee et al. (2020) extended the model to multivariate threshold:

where $\mathbf{f}_{t}$ is a vector of “factors” and $\boldsymbol{\gamma}$ is the corresponding unknown coefficients. So the model introduces a regime change due to a single index of factors. Allowing multivariate thresholding is important, because it permits the structural change to be governed by a potentially much larger dataset: $\mathbf{x}_{t}=\mathbf{B}\mathbf{f}_{t}+\mathbf{u}_{t},$ where $\dim(\mathbf{x}_{t})=N\to\infty$. So $\mathbf{f}_{t}$ can be unobserved factors that can be learned from $\mathbf{x}_{t}.$ For the identification purpose, suppose $\frac{1}{T}\sum_{t}\mathbf{f}_{t}\mathbf{f}_{t}^{\prime}=\mathbf{I}$ and $\mathbf{B}^{\prime}\mathbf{B}$ is diagonal, then $\boldsymbol{\gamma}$ and $\mathbf{f}_{t}$ are separately identified. This gives rise to the factor-driven two-regime regression model.

A natural strategy to estimate the model is to rely on least squares:

where $\widehat{}\mathbf{f}_{t}$ is the plugged-in PC-estimator of factors. Because the least squares problem is neither convex nor smooth in $\boldsymbol{\gamma}$, the computational task is demanding. Lee et al. (2020) recommended using algorithms based on mixed integer optimization (MIO). Introduce integers $d_{t}:=1\{\boldsymbol{\gamma}^{\prime}\widehat{}\mathbf{f}_{t}>0\}\in\{0,1\}$. The goal is to introduce linear constraints with respect to variables of optimization. Suppose there are known upper and lower bounds for $\delta_{j}$: $L_{j}\leq\delta_{j}\leq U_{j}$, where $\delta_{j}$ denotes the $j$th element of $\boldsymbol{\delta}$. Define $M_{t}\equiv\max_{\boldsymbol{\gamma}\in\Gamma}|\boldsymbol{\gamma}^{\prime}\widehat{}\mathbf{f}_{t}|$, where $\Gamma$ is the parameter space for $\boldsymbol{\gamma}$. Then it can be verified that the least squares problem is numerically equivalent to the following constraint MIO problem:

subject to (for any $\epsilon>0$), for each $t=1,\ldots,T$ and each $j=1,\ldots,\dim(\mathbf{w}_{t})$,

Then, we can apply modern MIO packages (e.g., Gurobi) to solve for the optimal $(\boldsymbol{\beta},\boldsymbol{\delta},\boldsymbol{\gamma})$.

Finally, Lee et al. (2020) also derived the asymptotic distribution of the estimated coefficients and proposed inferences based on bootstraps. Under the condition that $T=O(N)$, they showed that the effect estimating factors is negligible on the asymptotic distribution of the estimated $(\boldsymbol{\beta},\boldsymbol{\delta})$, but would affect both the rate of convergence and the limiting distribution of the estimated $\boldsymbol{\gamma}$.

7 Community detection

where $w_{k,l}$ is an unknown probability. We observe the matrix $\mathbf{A}$ and aim to recover the membership $\pi_{i}$ and the probabilities $w_{k,l}$ for all $k,l=1,\cdots,r$.

Therefore, $\mathbf{A}$ has the familiar decomposition (2.4), with $\mathbf{L}$ being similar to the systematic risk and $\mathbf{B}$ as a low-rank loading matrix. Since the elements in $\mathbf{S}$ are independent with mean-zero (Wigner matrix), the operator norm $\|\mathbf{S}\|$ does not grow too fast, compared to that of $\mathbf{L}$. We can then apply PCA on $\mathbf{A}$ to estimate $\mathbf{B}$. Suppose $r$ is known, then the estimator $\widehat{}\mathbf{B}$ is defined as $\sqrt{N}$ times the eigenvectors of $\mathbf{A}$, corresponding to the first $r$ eigenvalues.

Theorem 2.1 can be applied to obtain a deviation bound for the estimated loading matrix. If there is a sequence $g_{N}\to\infty$ and constants $c_{1},\cdots,c_{r}>0$ such that the eigenvalues $\lambda_{i}(\mathbf{W}^{1/2}\mathbf{B}^{\prime}\mathbf{B}\mathbf{W}^{1/2})=c_{i}g_{N}(1+o_{P}(1))$ for all $i\leq r$, then there is an $r\times r$ matrix $\mathbf{H}$, so that

Therefore, elements of a rotated $\mathbf{B}$ can be estimated uniformly well. Moreover, because each community has many nodes belong to, $\mathbf{B}\mathbf{H}$ has many identical rows, which makes the cluster analysis as a natural method for community detections. For instance, we can apply either the K-means cluster analysis, or the homogeneous pursuit of Ke et al. (2015) on the rows of $\widehat{}\mathbf{B}$ to consistently identify the communities.

8 Time varying models

So far we have been assuming that the factor loading and covariance matrices are time-invariant. Research on conditional factor models has also grown rapidly in recent years. Suppose

where $\mathbf{b}_{i,t}$ is a time-varying vector of loadings. There have been several approaches to addressing the issues of time-varying loadings. In this section we briefly review three of the most commonly used ones: (1) time-varying characteristics, (2) time-smoothing and (3) continuous-time models.

The first approach models $\mathbf{b}_{i,t}$ using a function of observed characteristics $\mathbf{z}_{i,t-1}$:

where $\mathbf{b}_{i}(\cdot)$ is either a linear function or an unknown nonparametric function of the characteristics. Therefore, the time-varyingness is mainly captured by the characteristics. An advantage of this approach, over the other two approaches to be reviewed below, is that if $\mathbf{z}_{i,t-1}$ is correctly specified and indeed can fully capture the degree of time-varyingness of the model, then $\mathbf{b}_{i,t}$ allows a large degree of varyingness, and potentially, structural breaks. On the other hand, the limitation of this approach is the potential misspecification of $\mathbf{z}_{i,t-1}$ and omitted variable problems. Above all, we refer to Gagliardini et al. (2019) for an excellent review on conditional factor models using this approach, and their applications in empirical asset pricing.

8.2 Time-smoothing

The second approach assumes that factor loadings change smoothly over time. Suppose $\mathbf{b}_{i}(\cdot)$ is an unknown smooth function, we assume

Then locally, $\mathbf{b}_{i,t}\approx\mathbf{b}_{i,r}$ for all $t\approx r$. So in a local window $\mathcal{B}(r)$ of each fixed $r$, the model is approximately time invariant:

Motivated by this assumption, Ang and Kristensen (2012) and Ma et al. (2020) tested the market mean-variance efficiency assumption in the case of known factor case. In the unknown factor case, Su and Wang (2017) first applied local smoothing on $y_{i,t}$ then employed PCA on the smoothed data to estimate the factors and loadings. While this approach does not require the specification of time-varying characteristics, it restricts to the smooth varying scenario and thus rules out structural breaks. In addition, slow rates of convergence appear near boundaries (that is, the beginning and the end of observing periods).

8.3 High-frequency factor models

where $\mathbf{y}_{t}$, $\mathbf{f}_{t},\mathbf{u}_{t}$ are vectors of asset prices, factors and idiosyncratic risks; $\boldsymbol{\alpha}_{t}$ is a drift term. The time-varying loading matrix $\mathbf{B}_{t}$ is an $N\times r$ matrix that is assumed to be continuous and locally bounded Itô semimartingale of the form:

where $\widetilde{}\boldsymbol{\alpha}_{s}$ and $\boldsymbol{\sigma}_{s}$ are optional processes and locally bounded; $\mathbf{W}_{s}$ is a Brownian motion. Roughly speaking, by the Burkholder-Davis-Grundy inequality (cf. chapter 2 of Jacod and Protter (2011)), $\mathbf{B}_{t}$ is also locally time-invariant, which is similar to the treatment of the time-smoothing approach. The major difference though, is that the use of high-frequency data has automatically “smoothed” the data. We refer to the following papers for recent developments on high-frequency factor models, among others: Aït-Sahalia and Xiu (2017); Chen et al. (2019a); Liao and Yang (2018); Li et al. (2019); Pelger (2019).

Unbalanced Panels

Missing data and unbalanced panels are not uncommon in economic and financial studies. Addressing the missing data issue in statistical modeling belongs to a larger category of problems, known as matrix completion. Low-rank matrix completion refers to the problem of recovering missing entries from low-rank matrices. It is particularly relevant to empirical asset pricing factor models, because many time series of returns have short histories or missing records. In this section we review several methods for matrix completions, which assume that the missing is at random, except for Cai et al. (2016); Bai and Ng (2019). Besides, the EM algorithm is also a classical approach to dealing with unbalanced panels. We refer to Stock and Watson (2002b); Su et al. (2019); Zhu et al. (2019) for detailed discussions on related issues.

Recall that the covariance matrix of $\mathbf{y}_{t}$, under the factor model (3.1), has the following decomposition, $\boldsymbol{\Sigma}_{y}=\mathbf{B}\operatorname{cov}(\mathbf{f}_{t})\mathbf{B}^{\prime}+\boldsymbol{\Sigma}_{u},$ where columns of $\mathbf{B}$ are approximately equal to the eigenvectors of $\boldsymbol{\Sigma}_{y}$ corresponding to the first $r$ eigenvalues. As such, let $\widehat{}\boldsymbol{\Sigma}_{y}$ be an input matrix, serving as an estimator for $\boldsymbol{\Sigma}_{y}$. Then as described in Section 2.2, we can estimate the space spanned by $\mathbf{B}$ using the leading eigenvectors of $\widehat{}\boldsymbol{\Sigma}_{y}$.

In the presence of missing data with exogenous missing, let $x_{it}=1\{y_{it}\text{ is observed}\}$ and we only observe $y_{it}x_{it}$ for all $(i,t)$, in which unobserved data is set to zero. Suppose for now $w_{i}:=P(x_{it}=1)$ is known. We can construct an unbiased estimator $\widehat{}\boldsymbol{\Sigma}_{y}=(\widehat{\sigma}_{ij})$ with

In the matrix form, let $\mathbf{Y}$ and $\mathbf{X}$ be the $N\times T$ matrices of $y_{it}$ and $x_{it}$. So we only observe $\mathbf{Y}\circ\mathbf{X}$, where $\circ$ represents the element-wise matrix product, the Hadamard product. Also let $\mathbf{W}$ be the diagonal matrix with $w_{i}$ being its $i$ th diagonal entry. Then

Therefore, columns of the loading matrix estimator $\widehat{}\mathbf{B}$ equal to $\sqrt{N}$ times the top right singular vectors of $\mathbf{Z}$. This method simply replaces the missing entries of $\mathbf{Y}$ by zero, and apply the inverse probability weighting (IPW) before applying PCA. The IPW has been popularly used in the causal inference literature (e.g., Imbens and Rubin (2015)). Here the same idea is applied to create an unbiased estimator for the covariance matrix.

In practice, we shall replace $w_{i}$ by its consistent estimators, such as $\widehat{w}_{i}:=\frac{1}{T}\sum_{t=1}^{T}x_{it}$. But in the case of homogeneous missing, that is, $w_{1}=\cdots=w_{N}$, the IPW is not needed, because $\mathbf{W}$ equals the identity matrix up to a constant, which does not affect the PCA on $\mathbf{Y}\circ\mathbf{X}$. In addition, factors can be further estimated using least squares by regressing $y_{it}x_{it}$ on the estimated loadings.

Theoretical properties were studied by Abbe et al. (2020); Su et al. (2019) under the assumption of homogenous missing. Su et al. (2019) used this estimator as their initial value for the EM algorithm. Xiong and Pelger (2019) allowed heterogenous missing and proved that the estimators are also asymptotically normal (they estimated $w_{i}w_{j}$ directly by $\frac{1}{T}\sum_{t=1}^{T}x_{it}x_{jt}$). We can also quickly derive the rate of convergence by applying Theorem 2.1. However, the IPW is the least efficient approach among all the methods to be discussed in this section. We shall verify this in a simulation study in Section 5.5.

2 Regularized matrix completion

Regularized matrix completion is a powerful technique to recover missing entries from low-rank matrices. This approach is also much faster than the EM algorithm in handling large panels. Due to these nice properties, it has also attracted much attention in the recent econometrics literature, e.g., Athey et al. (2018); Bai and Ng (2017); Moon and Weidner (2018); Giglio et al. (2020).

In the matrix form $\mathbf{Y}=\mathbf{M}+\mathbf{U}$, the goal is to recover the factor component $\mathbf{M}=\mathbf{B}\mathbf{F}^{\prime}$ when $\mathbf{Y}$ has missing elements. The nuclear-norm regularization is directly applicable:

with tuning parameter $\lambda$. The factors and loadings can be estimated by taking the singular vectors of $\widehat{}\mathbf{M}$. Negahban and Wainwright (2011) and Koltchinskii et al. (2011) derived the rate of convergence under the Frobenius norm. Under suitable conditions (e.g., missing at random, restricted strong convexity, sufficiently large noise) it can be proved that

Chen et al. (2020b) certifies further that the convex optimization (5.1) is optimal for all noise levels under Frobenius norm, operator norm, and elementwise-infinity norm. The proof is based on a novel technical device that bridges the convex optimization with a nonconvex optimization problem. However, this estimator is not asymptotically normal due to the presence of shrinkage bias, so is not suitable for statistical inferences.

3 Debiased estimators

Several recent progress in this literature focuses on debiasing the regularized regression in order to have valid confidence intervals, e.g., Chen et al. (2019b); Xia and Yuan (2019); Chernozhukov et al. (2019). When the missing is homogeneous, $P(x_{it}=1)=p$ for all $(i,t)$, Chen et al. (2019b) proposed the following simple debiased estimator

where $H_{R}(\cdot)$ is the best rank $R$ approximation in (2.5), $\widehat{\mathbf{M}}$ is given by (5.1), $\widehat{p}$ is the sample proportion of missing data. The idea is very intuitive. Ignoring the weak-dependence between $\widehat{\mathbf{M}}$ and $\mathbf{X}$ and estimating error in $\widehat{p}$, we have

which is approximately unbiased. However, the estimator $\widehat{\mathbf{M}}+\widehat{p}^{-1}(\mathbf{Y}-\widehat{\mathbf{M}}\circ\mathbf{X})$ is no longer of rank $R$, which increases the variances. This leads to use the projection as in (5.2), which is asymptotically efficient in terms of both rate and pre-constant.